建筑物矢量辅助的正射影像镶嵌线网络选择方法

李朋龙，董怡储，谭攀，李晓龙

(1.重庆市地理信息中心，重庆 401121；2.重庆市遥感中心，重庆 401121；3.重庆市生产力发展中心，重庆 401121)

0 引言

数字正射影像(digital orthophoto map，DOM)是数字遥感影像利用数字高程模型(digital elevation model，DEM)，经过单张影像数字微分纠正，然后进行多片镶嵌得到的[1-2]。影像镶嵌是正射影像制作中最为关键的一步，而镶嵌线的选择更是影像镶嵌中的重中之重，镶嵌线的好坏直接影响着镶嵌后正射影像的质量[3-4]。好的镶嵌线应该避免穿过建筑物等区域形成几何错位，并且镶嵌线两侧色彩过渡自然。目前国内外大多数DOM制作软件仍需要人工对镶嵌线进行编辑使其绕过建筑物等区域。因此研究镶嵌线自动选择方法对于提高DOM镶嵌的自动化程度和质量都有重要的意义。目前，镶嵌线网络的自动选择主要有三类方法：

第一类是基于广义的Voronoi图的镶嵌线网络。自Hsu等提出以常规Voronoi图的方法生成镶嵌线网络[5]到潘俊等提出顾及重叠面的广义Voronoi图镶嵌线网络[6-8]，学者们解决了常规Voronoi图在低重叠度下镶嵌结果不完全覆盖的问题，最大程度保留了镶嵌结果投影差理论最小的特点。但此类方法中镶嵌线没有自动绕过建筑物等区域。

第二类是通过构建重叠区域异差影像进而自动选择镶嵌线。Kerschner等通过计算重叠区域对应像元的亮度和梯度差异，用Twin Snakes 算子动态滑动寻找最佳路径。但是该算子通常会丢失全局最小值区域，不能确保全局最优性[9-10]。王琳等对该算法进行了优化，有效地减弱了这种影响，但是仍不能完全避免[11]。Davis提出了Dijkstra算法，但复杂度极高且容易丢失最佳路径[12]。Jaechoon Chon和袁修孝等提出并改进最小化最大算法的搜索策略来优化Dijkstra算法，很好地抑制了局部最大灰度，提高了镶嵌线搜索的效率[13-14]。张剑清等提出的基于蚁群算法进行镶嵌线自动选择算法能避开房屋、树木、水域等成像反差较大的区域[15-16]。陈继溢等提出了采用最优生成树(MST)的镶嵌线智能检测算法[17-18]。

第三类方法是利用DSM、DLG等数据辅助镶嵌线的自动选择。孙杰等提出结合DSM与DEM，通过作差比较即可获得测区内的建筑物等区域，从而使镶嵌线避开建筑物区域[19-20]。左志权等提出了基于DSM辅助数据进行影像分割，然后采用贪吃蛇法来自动选择镶嵌线的方法[21]。王东亮、汤竞煌等则利用历史DLG数据辅助镶嵌线网络的自动选择，能够避开大多建筑物区域[22-23]。

综上所述，第一类方法能够实现影像局部投影变形最小，但无法保证镶嵌线自动绕过建筑物等区域；而第二类和第三类方法绕过了建筑物，并没有考虑到中心投影投影变形中央最小的特点。在已有工作的基础上[24]，本文对提出的基于建筑物矢量数据对Voronoi图镶嵌线网络进行优化得到整个测区绕过建筑物的镶嵌线网络的方法，进行详细的论述，并与Inpho、IPS等软件进行了对比分析。通过对比分析本方法得到镶嵌线网络不仅绕过了建筑物成像区域，并且最大程度保留了Voronoi图局部影像投影差最小的特点。

1 基于建筑物矢量获取建筑物DOM成像区域

镶嵌线网络是相对于整个测区的，是对测区物方区域的划分，而建筑物矢量只是表达了建筑物在物方平面的区域范围。由于DEM不包含建筑物等具有高度的地物，因此DOM上建筑物依然存在投影变形，其成像区域远远大于建筑物矢量线所包围的区域。因此要想对镶嵌线进行优化使其绕过建筑物，必须计算出建筑物在每张DOM上的成像区域。

1.1 房顶矢量获取简易三维模型

规划管理部门每年都对城市建筑物进行信息采集工作，这些信息中包括了建筑物矢量线及高度等信息。可以根据建筑物矢量及建筑物高度和测区DEM得到建筑物的简易模型。如图1所示，ABCD为建筑物基底矢量多边形其高度信息可由DEM上内插得到，然后再加上建筑物高度H即可得到建筑物真实房顶abcd的位置，则ABCD-abcd即为该建筑物简易的三维模型(digital building model，DBM)。并将其以多面体结构模型的方式存储。

图1 建筑物矢量结合DEM获取简易三维模型

1.2 基于三维模型获得建筑物成像区域

为准确得到建筑物在DOM上的完整成像区域，本文利用由DBM和DEM相结合来获取建筑物在正射影像上成像区域矢量多边形的方法[24-26]。如图2所示，点O是摄影中心，长方体ABCD-abcd假设为简易建筑物模型，将其所有的三角面片按照公式(1)投影到原始影像上，然后求并集即可得到该建筑物在原始影像上成像多边形A′B′C′cba(多边形端点为像方坐标表示)。

图2 建筑物DOM成像区域获取原理

(1)

然后再将原始影像上多边形A′B′C′cba，按照公式(2)并基于DEM迭代逼近求解得到该多边形在正射影像上对应区域，即该建筑物在DOM上成像区域多边形A″B″C″D″c″b″a″(该多边形端点由物方坐标表示)。

(2)

2 镶嵌线网络自动选择选优化

中心投影的成像方式和地形的起伏使航空影像存在投影变形，且越往影像边缘投影变形越大。本文以每张影像像主点对应地面点坐标为散点集生成的Voronoi图为初始镶嵌线网络，再利用建筑物在DOM上的成像区域对初始Voronoi图镶嵌线网络进行优化选择，即绕开了建筑物的成像区域同时最大程度上保留了Voronoi图投影变形理论最小的特点。

2.1 影像与建筑物从属关系确定

一张影像上有多栋建筑物，一栋建筑物也会出现在多张影像上，为了对镶嵌线进行选择优化，需要确定测区影像与建筑物之间的从属关系。首先，利用将测区所有原始影像4个角点基于DEM迭代逼近的方法得到对应单片正射影像的成像区域多边形S1，S2，…，Sn；然后，对影像Pi将测区所有建筑物按照1.2节所述方法计算其在Pi上成像区域多边形B1，B2，…，Bm，进而与Si进行求交运算，得到影像Pi上所有成像建筑物Bi1，Bi2，…，Bim；最后，遍历所有影像包含的建筑物，得到每栋建筑物如Bj对应的所有成像影像Pj1，Pj2，…，Pjm。至此，我们建立了测区所有影像与所有建筑物之间的双向从属关系。

2.2 初始镶嵌线网络的自动生成与简化

以测区所有影像像主点对应的地面点为散点集P={p1，p2，…，pn}(n≥3)，自动生成测区初始的Voronoi图镶嵌线网络。由于传统Voronoi图是纯粹基于几何位置关系生成的，因此网络中会出现很多短边。

如图3(a)所示，如果这些边的两端点全部落在建筑物成像的区域内则会加大对镶嵌线网络优化的难度。对于航带规则、重叠度合理的序列航空影像生成的Voronoi图，这些短边会规律地出现于中心影像与非四周影像(四周指上下左右)的边界处。因此可以删除这些短边对初始镶嵌线网络进行一定的简化，既可以最大程度地保证Voronoi图的特性又可以简化镶嵌线后续的优化难度。本文将Voronoi多边形中长度小于该多边形最长四条边长平均值的四分之一的短边删除，以该短边两端点的几何中点替代该短边，简化后的Voronoi图网络如图3(b)所示，每个节点最多归4个多边形共有，每条边最多归2个多边形共有。

图3 测区Voronoi图镶嵌线网络局部

2.3 镶嵌线网络中节点的选择优化

要想使镶嵌线自动绕过建筑物成像区域，必须先保证网络中节点不落在建筑物的成像区域内。由2.1节可知，每个节点最多归4个多边形共用，因此要保证节点不在4张影像中任一张影像上建筑物的成像区域内。

由于成像位置和姿态的不同，一栋建筑物在每张正射影像的成像区域也不相同，如图4所示，O1，O2是相邻两张影像的摄影中心，棱柱体ABCD-abcd为建筑物，则该建筑物在正射影像O2上成像区域为多边形A′B′C′cba，而在正射影像O1上成像区域为多边形B″C″D″dab。则该建筑物在2张影像的成像总区域为多边形A′B′EB″C″D″da[24]。

图4 建筑物在不同影像上成像示意图

以下给出了镶嵌线网络所有节点的优化过程：

①判断节点P是否归2张以上影像有效多边形共有(以下以4张共有为例)，如是进行步骤②，否则进行步骤⑦。

②以节点P为中心，搜索方圆一定距离范围的所有建筑物，记为B1，B2，…，Bn。

③对于每个建筑物，计算其在这4张DOM上的成像区域，并求并集作为该建筑物的整体城像区域，记为S1，S2，…，Sn。

④对求解出的S1，S2，…，Sn进行求交判断，如果有交集则合并同时记下Si与Bi的对应关系，直到剩下S1，S2，…，Sm(m≤n)的任意二者都不相交。

⑤逐一判断节点P是否在S1，S2，…，Sm(m≤n)内部，如果在则记下该房屋的编号Bi，同时找到Bi对应的合并后Si，如果都不在，则该节点无需优化进入步骤⑦。

⑥求出包含多边形Si的最小外包矩形，然后找出距离节点P距离最小的一个矩形端点Q作为节点P的优化后位置。然后将镶嵌线网络中所有的节点P更新为点Q。

⑦判断P是否为最后一个节点，如果是则优化完毕，如果不是则进行下一个节点的检测更新，进入步骤①。

2.4 镶嵌网络中镶嵌线的选择优化

每条镶嵌线由2个影像有效多边形共有，因此对于镶嵌线自动绕过建筑物的优化需要同时考虑2张影像共有的建筑物成像情况以及最优化路径选取。

以下给出了单条镶嵌线优化的方法步骤：

①判断镶嵌线L是否为2个影像有效多边形共有，若仅为一个多边形拥有，则该镶嵌线无需优化，直接进入步骤⑨。

②根据镶嵌线L与影像的对应关系找到该镶嵌线两侧的影像序号Pici，Picj，并根据建筑物与影像间的对应关系，将同时出现在该2张影像上的建筑物挑选出来记为B1，B2，…，Bn。

③对于每个建筑物，计算其在Pici和Picj2张正射影像上的成像区域，并求并集作为该建筑物的整体城像区域，记为S1，S2，…，Sn。

④判断每个建筑物成像区域Si是否与镶嵌线L相交，将相交的多边形留下记为：S1，S2，…，Sm(m≤n)，若都不相交则该镶嵌线不需要优化处理，进入步骤⑨。

⑤对S1，S2，…，Sm进行两两求交判断，如果有交集则合并，直到剩下S1，S2，…，Sl(l≤m)的任意二者彼此之间不相交。

⑥将多边形S1，S2，…，Sl(l≤m)按照距离镶嵌线L起点LB的距离按照从小到大排序得到多边形序列SS1，SS2，…，SSl(l≤m)。

⑦从第一个多边形SS1起，求出镶嵌线L与多边形SSi(i≤l)的交点，记第一个交点为INSf，最后一个交点为INSe，则镶嵌线从点INSf到INSe这段沿着多边形SSi(i≤l)的轮廓走，选择总距离最短的路径代替原镶嵌线中从点INSf到INSe的部分，如图5所示。直到该条镶嵌线绕过所有的多边形SSi(i≤l)。

⑧更新2张影像多边形中的该条镶嵌线。

⑨进入下一条镶嵌线的优化，直到所有的镶嵌线优化完毕，得到整个测区自动绕过建筑物成像区域的镶嵌线网络。

图5 镶嵌线绕过建筑物优化示意图

3 实验与分析

为验证本文方法的有效性，以某地4个航带共32张由SWDC-5航摄仪拍摄的分辨率为0.1 m、航向重叠约为60%，旁向重叠约为45%、框幅为8 230×6 168的航摄影像及测区分辨率为2 m的DEM和656个建筑物房顶矢量进行实验。

3.1 建筑物成像区域计算结果

图6给出了利用本文方法从建筑物房顶矢量获得其在单张DOM上成像区域的过程和结果，其中图6(a)为测区一栋建筑物的房顶矢量，图6(b)为将其投影到DEM上得到的建筑物简易三维模型，图6(c)为由建筑物DBM投影到原始影像上的成像范围，图6(d)为基于DEM迭代逼近得到的建筑物在单片DOM上成像区域的矢量范围。从图中可以看出，利用本文方法可以由建筑物房顶矢量准确获得建筑物在单片DOM上的成像范围。

图6 获取建筑物成像区域物方矢量过程

3.2 镶嵌线网络自动选择分析

图7给出了镶嵌线网络中节点优化的结果，其中图7(a)～图7(d)是同一栋建筑物在节点相关4张单片DOM上的成像区域(黑色多边形)。图7(d)中外围多边形为该建筑在4张影像的总成像区域多边形，红圈内为节点优化前的位置，图7(e)为节点优化后的位置。从图中可以看出，利用本文的方法成功将落在建筑物成像区域内的网络节点检测并移出，为镶嵌线的优化奠定了基础。

经过节点优化后的镶嵌线网络中每条镶嵌线最多归2个多边形共有，因此镶嵌线的优化只需考虑镶嵌线穿过的房屋在对应的2张影像的成像区域的并集即可。如图8所示，一条镶嵌线连续穿过了3个建筑物的成像区域，图8(a)为镶嵌线与镶嵌线对应上方影像建筑物成像区域的位置关系，图8(b)为镶嵌线与对应下方影像上建筑物成像区域的位置关系，图8(c)为原始镶嵌线镶嵌结果，图8(d)是优化后的镶嵌线镶嵌结果。比较图8(c)和图8(d)可知，利用本文方法优化后的镶嵌线成功绕过了建筑物成像区域，避免了镶嵌线两侧纹理错位的发生。

图7 网络节点优化过程及结果

图8 镶嵌线优化过程与结果

图9给出了测区整体镶嵌线网络优化前后的套合图与局部对比图，可以看出利用本文方法优化选择后的镶嵌线不仅绕过了建筑物的成像区域保证了镶嵌结果的纹理质量，同时也最大程度上保留了Voronoi图镶嵌线网络局部投影变形理论最小的特点，保证了镶嵌结果的几何精度。

图9 测区镶嵌线网络优化前后对比

3.3 与其他方法对比分析

本文利用Inpho OrthoVista 和IPS 4.2(icaros photogrammetric suite)的Stitch 模块分别进行了镶嵌线网络的自动选择，图10给出了3种方法镶嵌线网络的局部，表1给出了3种方法以该局部区域中心影像有效区域多边形(与周边影像镶嵌线组成的多边形)为对象的量化值统计情况。

图10 与其他镶嵌线网络自动选择方法结果对比

表1 与其他方法对比

首先，在镶嵌线网络选择效果上，从图10可以看出OrthoVista 自动选择的镶嵌线有5处穿过了建筑物，IPS 自动选择的镶嵌线有4处穿过建筑物(该片区共有6处)，造成了镶嵌线两侧纹理错位，而本文方法该片区镶嵌线网络绕过了所有的建筑物成像区域。特别是该片区右上方建筑物间距极小、密集程度高的情况，大多数自动选择方法无法绕过，而本文方法在镶嵌线选择优化过程中通过合并有交集的多个建筑物成像区域的方法，有效避开了建筑物密集区域。

其次，在测区镶嵌线网络的整体布局上，本文方法是对Voronoi图进行优化选择，最大程度保留了其镶嵌结果投影变形理论最小的优点。实验中通过测量建筑物同一房顶点与对应基底点的距离获得其投影变形大小，从表1可以看出，利用本文方法镶嵌结果中建筑物B1的投影变形最小为1.99 m，而OrthoVista镶嵌结果投影变形最大达到8.24 m。对于建筑物B2，本文方法和OrthoVista镶嵌结果中建筑物B2的纹理均取自同一影像，因此投影变形大小一致为1.86 m，而IPS Stich的镶嵌结果则取自于另一张影像，投影变形为3.47 m。

同时，本文算法基于矢量计算得到的镶嵌线网络，而另外2种算法是构建异差影像在像素级寻找最优路径的方法获得镶嵌线，因此本文方法算法更为简单，同时获得的镶嵌线网络为面结构，且节点数目更少。如表1所示，利用本文方法得到中心影像的有效范围矢量只有347个节点，远远低于前2种方法的1 370和735个节点。

最后，图11给出了3种方法自动选择生成整个测区镶嵌线网络的耗时，其中前2种方法是以单片正射影像为基础数据，对相邻影像进行全幅镶嵌线自动搜索，得到航向和旁向2个方向上相邻影像的镶嵌线，再进行镶嵌线相互切割得到最终镶嵌线网络。本文方法以建筑物矢量及高度、测区DEM、影像定向参数为基础数据，先生成Voronoi图初始镶嵌线网络并简化，再根据建筑物成像区域进行镶嵌线网络的选择优化。OrthoVista软件耗时1 925 s，IPS软件Stitch模块由于利用了GPU加速技术耗时为963 s，而本文方法由于计算复杂度远远低于前2种方法，且多为矢量运算，从数据读入到镶嵌线网络选择优化完成共耗时349 s。

图11 耗时对比

4 结束语

本文提出了一种建筑物矢量辅助的城市大比例尺正射影像镶嵌线网络自动选择方法。首先得到建筑物在单片正射影像上的成像区域，然后自动检测并优化初始Voronoi图镶嵌线网络中的所有节点和镶嵌线，得到整个测区绕开建筑物成像区域的最优镶嵌线网络。该方法不仅能够快速高效地得到绕过建筑物的镶嵌线网络，并且保留了Voronoi图网络镶嵌结果投影变形理论最小的特点，保证了镶嵌后DOM影像的质量，为城市大比例尺正射影像镶嵌线网络选择提供了一种方法。下一步将研究根据左右差值影像对镶嵌线网络进行优化处理，减少对建筑物房顶矢量的依赖。