曾 宇，苏晓斌

(湖南省第二测绘院，湖南 长沙 410019)

ADS100机载数字航空摄影测量系统是目前最先进的推扫式机载数字航空摄影测量系统[1]。其航飞一个架次可以同时获取前视、下视和后视3个视角的影像，具有分辨率高、光谱特性良好的特点。由于其独特的优势，在不动产测绘中得到广泛应用。

在处理ADS数据过程中，常用的处理方法是：首先用XPro软件进行空三加密及DSM匹配，其次将匹配的DSM导入第三方软件进行人机交互编辑处理，最后利用处理后的DEM生成正射影像。目前国内外软件在空三加密和DSM匹配环节上技术均较成熟，但从DSM到DEM的技术不够完善，因此自动提取的DSM仍需大量的人工编辑后才能满足用户实际需求[2]。

本文所介绍的方法是利用华浩超算平台对DSM进行编辑。该平台DSM编辑模块的原理为：首先通过SGBM匹配算法提取DSM点云数据；其次应用影像提取算法提取人工地物和水域边界；最后将提取的数据作为DSM点云处理矢量约束条件进行DSM处理，以达到有效去除人工地物的目的，生成满足需要的DEM数据。但在自动提取过程中难免存在错漏，故采用地理国情普查工作中的房屋矢量边界数据来替代和补充自动提取的人工地物矢量数据。通过该方法有效解决了DSM编辑难度大、效率低、精度难以达标的难题。

1 DEM和DOM

不动产测绘，目前最符合国情的手段是利用航飞影像生产数字线划图(DLG)和数字正射影像(DOM)。DLG可以直接在立体环境下采集得到，而DOM则需要以下步骤才能得到：①匹配DSM；②编辑DSM；③正射纠正。

DEM是通过有限的地形高程数据实现对地形曲面数字化模拟(即地形表面形态的数字化表达)，它是用一组有序数值阵列形式表示地面高程的一种实体地面模型。

DSM是指包含了地表建筑物、桥梁和树木等高度的地面高程模型。

数字地面模型(digital terrain model,DTM)是地形表面形态属性信息的数字表达，是带有空间位置特征和地形属性特征的数字描述。

DOM是根据像片的内外方位元素和DEM对数字化的航空像片/遥感影像进行数字微分纠正得到的影像图。但是DEM不能直接获取，因此需要对DSM或DTM进行编辑，使它接近真正的DEM。

由此可得，DOM质量取决于DEM的质量。而对DSM或DTM编辑又是决定DEM质量的重要步骤。

2 DSM的匹配与编辑

2.1 基于软件像素工厂的方法

像素工厂基于对多视角数据的自动多重相关匹配算法提取DSM[3]。像素工厂在平面环境下编辑DSM，程序实时纠正影像。

优点是：

(1) 能实时切换到正射影像界面，从而判断对DSM编辑是否有效。

(2) 能记录编辑的轨迹，实现多人同时编辑，不会产生重复编辑的现象。

缺点是：

(1) 需要频繁切换正射影像和DSM编辑界面。

(2) 编辑范围较大的DSM时，程序响应时间较长。

(3) 在平面环境下DSM基本上靠经验判断。

2.2 基于国产软件PixelGrid的方法

国产软件PixelGrid是基于物方几何约束的多影像相关匹配算法(geometrically constrained cross-correlation，GC3)提取DSM[1]。

该软件最为常用的编辑功能有：线编辑、区域平滑、曲面拟合和局部重构等。线编辑(选择特征线处理)主要用于沟谷、山脊处编辑;平滑处理主要用于处理平地的噪点;曲面拟合用于编辑倾斜的地形，也可用于山坡有植被区域的编辑;局部重构用于编辑小范围内的匹配错误，如阴影区域、雪覆盖区域。

优点是：

(1) 能够在立体环境下编辑DSM，实现高程点贴地。

(2) 能够叠加影像，编辑DSM，做到实时纠正影像。

缺点是：

(1) 需对数据进行转换，比较耗时。

(2) 立体环境下需要对人工地物进行单独处理，工序繁琐，效率低下。

2.3 基于华浩超算平台的DSM匹配编辑方法

华浩超算平台充分运用了SGBM匹配算法和图像区域增长算法匹配编辑DSM。

2.3.1 SGBM匹配算法

SGBM匹配算法是一种用于计算双目视觉中disparity的半全局匹配算法。具体思路：通过选取每个像素点的disparity，组成一个disparity map，设置一个与disparity map相关的全局能量函数，使这个能量函数最小化，以达到求解每个像素最优disparity的目的。

2.3.1.1 能量函数

能量函数形式如下

(1)

式中，D指disparity map；E(D)是该disparity map对应的能量函数；p、q代表影像中的某个像素；Np指像素p的相邻像素点(一般认为影像为8通道)；C(p，Dp)指当前像素点disparity为Dp时，该像素点的cost值；P1是一个惩罚系数，它用在像素p相邻像素的dsparity值与像素p的dsparity值相差为1的那些像素；P2是一个惩罚系数，它用在像素p相邻像素的dsparity值与像素p的dsparity值相差大于1的那些像素；I[]表示如果函数中的参数为真时返回1，否则返回0。

利用式(1)在一个二维图像中寻找最优解，耗时巨大，但是该求解方式可归为NP-complete问题。因此该问题可以被近似分解为多个一维问题，即线性问题，而且每个一维问题都可以用动态规划方法来解决。又因为1个像素有8个相邻像素，因此可分解为8个一维问题。故每个像素的disparity只与其左边的像素相关，有如下公式

Lr(p,d)=C(p,d)+min(Lr(p-r,d),Lr(p-r,d-1)+

(2)

式中，r表示某个指向当前像素p的方向，在此可以理解为像素p左边的相邻像素；Lr(p，d) 表示在沿着当前方向(即从左至右)的情况下，当像素p的disparity取值为d时，其最小的cost值。

而这个最小值是从以下4种可能的候选值中选取的最小值：

(1) 前一个像素(左相邻像素)disparity取值为d时，其最小的cost值。

(2) 前一个像素(左相邻像素)disparity取值为d-1时，其最小的cost值+惩罚系数P1。

(3) 前一个像素(左相邻像素)disparity取值为d+1时，其最小的cost值+惩罚系数P1。

(4) 前一个像素(左相邻像素)disparity取值为其他值时，其最小的cost值+惩罚系数P2。

另外，当前像素p的cost值还需减去前一个像素取不同disparity值时最小的cost值。这是因为Lr(p，d)会随着当前像素的右移不断增长。为了防止数值溢出，需要其维持在一个较小的数值，该最小值为C(p，d)。C(p，d)的计算公式如下

(3)

式(3)的思路为：

首先设像素p的灰度/RGB值为I(p)，同时从I(p)、(I(p)+I(p-1))/2、(I(p)+I(p+1))/2三个值中选择出和I(q)差值最小的值(即d(p,p-d))。其次从I(q)、(I(q)+I(q-1))/2、(I(q)+I(q+1))/2三个值中选择出和I(p)差值最小的值(即d(p-d,p))；最后从上述2个值中选取最小值，即C(p，d)。

以上是从一个方向(从左至右)计算出像素在取值为某一disparity值时的最小cost值。但是，一个像素有8个邻域，因此从8个方向进行计算(左右，右左，上下，下上，左上右下，右下左上，右上左下，左下右上)8个cost值。

计算完成后，将8个方向上的cost值进行累加，然后选取累加值最小的disparity值作为该像素的最终disparity值。

2.3.1.2 disparity map

按照上述方法对每个像素进行操作后，便形成了整个图像的disparity map。公式表达如下

(4)

该方法首先创建一系列的矢量图层，包括房屋层、河湖层、湖泊层、山地层，该矢量图层作为后期编辑DSM点云的约束条件，用以滤除点云，生成DEM。

2.3.2 图像区域增长算法

DSM自动提取生成DEM利用图像区域增长算法。区域增长算法的实现需要3个条件：

(1) 种子点的选取：种子点选取的位置为地形点。

(2) 增长条件：如果种子点的相邻点为地形点，则归为新增长区域点。如果相邻点落在房顶或树木上，则该点不纳入新的增长区域点。

(3) 终止增长：若区域内不存在增长区域点，则终止增长。

这样便可以生成关于地形和地物的标志二值影像。为了显示生成的地物区域，必须获得二值影像中有序的边界坐标值。

二值影像的边界一般采用“虫随法”进行搜索，但是传统的“虫随法”存在以下缺陷：

(1) 可能会忽略边界上的小突出部分，但这个缺陷并不是致命的。因为几个像素被忽略，并不会对最终的DOM产生重大影响。

(2) 对于边界不封闭的区域，难以获得全部边界；“虫随法”的“小虫”有可能会掉入陷阱(即“小虫”围绕某一个小局部区域边界一直爬行，永远回不到起点)，这是一个致命的缺陷。

因此本算法在传统的“虫随法”基础上作出如下改进：

(1) 边界跟踪从二值影像的左上角点开始逐点扫描，当遇到边界点时便开始进行跟踪直到回到开始遇到的起始边界点。

(2) 如果回不到起始点，则将搜索到的点与边界点或角点一起组成封闭区域。

由上述方法可得：每一个像素至少产生一个边界点(直线型边界也是如此)，根据道格拉斯线综合法，可设置距离参数优化边界，并将优化后的边界转化为矢量文件。图1为自动生成的矢量边界。

图1 提取的房屋和河流层矢量边界

以矢量边界为约束条件，过滤房屋。过滤后的区域通过边界上的高程值进行内插，内插方法为三角剖分算法。效果如图2—图3所示。

实时预览正射影像的处理效果，如图4—图5所示。

图2 DSM被过滤前的效果

图3 DSM被矢量过滤后的效果

图4 某区域过滤房屋前的正射预览图

图5 某区域过滤房屋后的正射预览图

2.4 利用国情数据基于华浩超算平台进行DSM编辑试验

2.4.1 背 景

(1) 华浩超算平台匹配DSM和快速生成DOM初级数据。

(2) 地理国情房屋层数据，图6为某区域房屋层数据。

图6 某区域地理国情房屋层数据

2.4.2 基于华浩超算平台的处理方法

导入矢量数据，对矢量数据进行以下处理：

(1) 在平面环境下套合初级DOM并对地理国情线划进行修改，使尽可能多边界点为地形点，如图7所示。

图7 房屋层线划套合初级DOM

由图7可见：初级DOM中房屋存在较大变形；地理国情数据中房屋线划走向都是沿着屋顶。

(2) 对修改的线划形成的面域进行扩充，进一步补充边界地形点；修改线划图导入DSM，如图8所示。

图8 修改线划图导入DSM

(3) 线划图叠加DSM作为约束条件进行DSM处理，如图9—图10所示。

图9 房屋层叠加分析(Overlay)

图10 房屋层缓冲区分析(Buffer)

对比图9与图10可得，DSM中房屋被过滤，生成DEM。

3 DOM生成

加载处理后的DEM生成DOM，如图11所示。

图11 新旧DOM对比

对比新旧DOM可得，房屋变形已经修改。

4 检测精度

检测精度公式为

式中，f为相机的焦距；ai、bi、ci(i=1，2，3)为相机摄影瞬间的3个外方位角元素；(XS，YS，ZS)为相机摄影中心坐标；(X，Y，Z)是物方坐标；(x，y)为像方坐标。

从上述公式可以看出去除地物高程可以有效提高DOM的精度。

套合立体环境采集线划，查看精度，如图12所示。

图12 DLG叠加DOM

方法：将数字线划图叠加对应的1∶2000正射影像，进行精度评定。

规定：1∶500 1∶1000 1∶2000地形图航空摄影测量内业规范(GB/T 7930—2008)。

精度指标：平面位置中误差，平地、丘陵地区为实地±1.2 m，山地、高山地区为实地±1.6 m，最大误差应不超过中误差的2倍。

采集量取方式：分别采集量取低矮房屋的屋檐角、高层房屋基角、斑马线、坎角及围墙角等特征点。

比较结果如图13所示，可知，精度满足要求。

5 结 语

通过探索和实践，本文可以得出以下结论：充分利用地理国情和基础测绘线划图数据搭载华浩超算平台，对DEM进行快速匹配与编辑，并生成合格的DOM的方式方法，也同样适用于框幅式数码航空影像。因此该方式方法具有一定应用和推广价值。

图13 DLG叠加DOM精度检测表

[1] 张力,张继贤.基于多基线影像匹配的高分辨率遥感影像DEM自动生成[J].测绘科学,2008,33(S2)：35-39.

[2] 尤红建,苏林,李树楷.利用机载三维成像仪的DSM数据自动提取建筑物[J].武汉大学学报(信息科学版),2002,27(4):409-413.

[3] 丁勇，周利平，司玉琴，等.利用PixelGrid软件实现航空影像高保真高效率DSM的生产[J].遥感信息，2015,30(3)：85-88.

[4] 殷福忠，刘红军，张延波.基于地理信息的数字影像采集集成系统研究[J].测绘与空间地理信息，2010,33(4):42-45.

[5] 朱继文，李清华.基于影像匹配数据获取DEM方法探讨[J].黑龙江工程学院学报,2009,23(1)：36-38.

[6] 赵双明，李德仁.ADS 40机载数字传感器平差数学模型及其试验[J].测绘学报,2006,35(4)：342-346.

[7] 马旭燕，崔会娟，张富玲.基于 PixelGrid 系统的 DEM 自动匹配后处理关键技术研究[J].测绘与空间地理信息，2016,39(9)：147-149.

[8] 周晓敏,杨爱玲，孙丽梅，等.浅议基于像素工厂的 ADS80 影像处理技术[J].测绘与空间地理信息，2012,35(5)：146-148.

[9] 张雪萍.基于Pixel Factory 的 ADS80 影像快速处理方法[J].测绘与空间地理信息，2013,36(12)：12-15.

[10] 王海涛，武吉军，冯聪军，等.徕卡ADS40/ADS80数字航空摄影测量系统[J].测绘通报，2009(10)：78-79.

[11] 高立.ADS80航空摄影测量系统的特点与应用[J].测绘与地理空间信息,2011,34(6):212-214.

[12] 张永生,范大昭,纪松.用于ADS40传感器的多视觉立体匹配算法模型[J].测绘科学技术学报,2007,24(2):83-86.

[13] 中华人民共和国国家质量监督检验检疫总局 中国国家标准化管理委员会1∶500 1∶1000 1∶2000 地形图航空摄影测量内业规范：GB/T 7930—2008[S].北京: 中国标准出版社,2008.