李朋龍，董怡儲，譚攀，李曉龍

(1.重慶市地理信息中心，重慶 401121；2.重慶市遙感中心，重慶 401121；3.重慶市生產力發(fā)展中心，重慶 401121)

0 引言

數字正射影像(digital orthophoto map，DOM)是數字遙感影像利用數字高程模型(digital elevation model，DEM)，經過單張影像數字微分糾正，然后進行多片鑲嵌得到的[1-2]。影像鑲嵌是正射影像制作中最為關鍵的一步，而鑲嵌線的選擇更是影像鑲嵌中的重中之重，鑲嵌線的好壞直接影響著鑲嵌后正射影像的質量[3-4]。好的鑲嵌線應該避免穿過建筑物等區(qū)域形成幾何錯位，并且鑲嵌線兩側色彩過渡自然。目前國內外大多數DOM制作軟件仍需要人工對鑲嵌線進行編輯使其繞過建筑物等區(qū)域。因此研究鑲嵌線自動選擇方法對于提高DOM鑲嵌的自動化程度和質量都有重要的意義。目前，鑲嵌線網絡的自動選擇主要有三類方法：

第一類是基于廣義的Voronoi圖的鑲嵌線網絡。自Hsu等提出以常規(guī)Voronoi圖的方法生成鑲嵌線網絡[5]到潘俊等提出顧及重疊面的廣義Voronoi圖鑲嵌線網絡[6-8]，學者們解決了常規(guī)Voronoi圖在低重疊度下鑲嵌結果不完全覆蓋的問題，最大程度保留了鑲嵌結果投影差理論最小的特點。但此類方法中鑲嵌線沒有自動繞過建筑物等區(qū)域。

第二類是通過構建重疊區(qū)域異差影像進而自動選擇鑲嵌線。Kerschner等通過計算重疊區(qū)域對應像元的亮度和梯度差異，用Twin Snakes 算子動態(tài)滑動尋找最佳路徑。但是該算子通常會丟失全局最小值區(qū)域，不能確保全局最優(yōu)性[9-10]。王琳等對該算法進行了優(yōu)化，有效地減弱了這種影響，但是仍不能完全避免[11]。Davis提出了Dijkstra算法，但復雜度極高且容易丟失最佳路徑[12]。Jaechoon Chon和袁修孝等提出并改進最小化最大算法的搜索策略來優(yōu)化Dijkstra算法，很好地抑制了局部最大灰度，提高了鑲嵌線搜索的效率[13-14]。張劍清等提出的基于蟻群算法進行鑲嵌線自動選擇算法能避開房屋、樹木、水域等成像反差較大的區(qū)域[15-16]。陳繼溢等提出了采用最優(yōu)生成樹(MST)的鑲嵌線智能檢測算法[17-18]。

第三類方法是利用DSM、DLG等數據輔助鑲嵌線的自動選擇。孫杰等提出結合DSM與DEM，通過作差比較即可獲得測區(qū)內的建筑物等區(qū)域，從而使鑲嵌線避開建筑物區(qū)域[19-20]。左志權等提出了基于DSM輔助數據進行影像分割，然后采用貪吃蛇法來自動選擇鑲嵌線的方法[21]。王東亮、湯競煌等則利用歷史DLG數據輔助鑲嵌線網絡的自動選擇，能夠避開大多建筑物區(qū)域[22-23]。

綜上所述，第一類方法能夠實現影像局部投影變形最小，但無法保證鑲嵌線自動繞過建筑物等區(qū)域；而第二類和第三類方法繞過了建筑物，并沒有考慮到中心投影投影變形中央最小的特點。在已有工作的基礎上[24]，本文對提出的基于建筑物矢量數據對Voronoi圖鑲嵌線網絡進行優(yōu)化得到整個測區(qū)繞過建筑物的鑲嵌線網絡的方法，進行詳細的論述，并與Inpho、IPS等軟件進行了對比分析。通過對比分析本方法得到鑲嵌線網絡不僅繞過了建筑物成像區(qū)域，并且最大程度保留了Voronoi圖局部影像投影差最小的特點。

1 基于建筑物矢量獲取建筑物DOM成像區(qū)域

鑲嵌線網絡是相對于整個測區(qū)的，是對測區(qū)物方區(qū)域的劃分，而建筑物矢量只是表達了建筑物在物方平面的區(qū)域范圍。由于DEM不包含建筑物等具有高度的地物，因此DOM上建筑物依然存在投影變形，其成像區(qū)域遠遠大于建筑物矢量線所包圍的區(qū)域。因此要想對鑲嵌線進行優(yōu)化使其繞過建筑物，必須計算出建筑物在每張DOM上的成像區(qū)域。

1.1 房頂矢量獲取簡易三維模型

規(guī)劃管理部門每年都對城市建筑物進行信息采集工作，這些信息中包括了建筑物矢量線及高度等信息。可以根據建筑物矢量及建筑物高度和測區(qū)DEM得到建筑物的簡易模型。如圖1所示，ABCD為建筑物基底矢量多邊形其高度信息可由DEM上內插得到，然后再加上建筑物高度H即可得到建筑物真實房頂abcd的位置，則ABCD-abcd即為該建筑物簡易的三維模型(digital building model，DBM)。并將其以多面體結構模型的方式存儲。

圖1 建筑物矢量結合DEM獲取簡易三維模型

1.2 基于三維模型獲得建筑物成像區(qū)域

為準確得到建筑物在DOM上的完整成像區(qū)域，本文利用由DBM和DEM相結合來獲取建筑物在正射影像上成像區(qū)域矢量多邊形的方法[24-26]。如圖2所示，點O是攝影中心，長方體ABCD-abcd假設為簡易建筑物模型，將其所有的三角面片按照公式(1)投影到原始影像上，然后求并集即可得到該建筑物在原始影像上成像多邊形A′B′C′cba(多邊形端點為像方坐標表示)。

圖2 建筑物DOM成像區(qū)域獲取原理

(1)

然后再將原始影像上多邊形A′B′C′cba，按照公式(2)并基于DEM迭代逼近求解得到該多邊形在正射影像上對應區(qū)域，即該建筑物在DOM上成像區(qū)域多邊形A″B″C″D″c″b″a″(該多邊形端點由物方坐標表示)。

(2)

2 鑲嵌線網絡自動選擇選優(yōu)化

中心投影的成像方式和地形的起伏使航空影像存在投影變形，且越往影像邊緣投影變形越大。本文以每張影像像主點對應地面點坐標為散點集生成的Voronoi圖為初始鑲嵌線網絡，再利用建筑物在DOM上的成像區(qū)域對初始Voronoi圖鑲嵌線網絡進行優(yōu)化選擇，即繞開了建筑物的成像區(qū)域同時最大程度上保留了Voronoi圖投影變形理論最小的特點。

2.1 影像與建筑物從屬關系確定

一張影像上有多棟建筑物，一棟建筑物也會出現在多張影像上，為了對鑲嵌線進行選擇優(yōu)化，需要確定測區(qū)影像與建筑物之間的從屬關系。首先，利用將測區(qū)所有原始影像4個角點基于DEM迭代逼近的方法得到對應單片正射影像的成像區(qū)域多邊形S1，S2，…，Sn；然后，對影像Pi將測區(qū)所有建筑物按照1.2節(jié)所述方法計算其在Pi上成像區(qū)域多邊形B1，B2，…，Bm，進而與Si進行求交運算，得到影像Pi上所有成像建筑物Bi1，Bi2，…，Bim；最后，遍歷所有影像包含的建筑物，得到每棟建筑物如Bj對應的所有成像影像Pj1，Pj2，…，Pjm。至此，我們建立了測區(qū)所有影像與所有建筑物之間的雙向從屬關系。

2.2 初始鑲嵌線網絡的自動生成與簡化

以測區(qū)所有影像像主點對應的地面點為散點集P={p1，p2，…，pn}(n≥3)，自動生成測區(qū)初始的Voronoi圖鑲嵌線網絡。由于傳統Voronoi圖是純粹基于幾何位置關系生成的，因此網絡中會出現很多短邊。

如圖3(a)所示，如果這些邊的兩端點全部落在建筑物成像的區(qū)域內則會加大對鑲嵌線網絡優(yōu)化的難度。對于航帶規(guī)則、重疊度合理的序列航空影像生成的Voronoi圖，這些短邊會規(guī)律地出現于中心影像與非四周影像(四周指上下左右)的邊界處。因此可以刪除這些短邊對初始鑲嵌線網絡進行一定的簡化，既可以最大程度地保證Voronoi圖的特性又可以簡化鑲嵌線后續(xù)的優(yōu)化難度。本文將Voronoi多邊形中長度小于該多邊形最長四條邊長平均值的四分之一的短邊刪除，以該短邊兩端點的幾何中點替代該短邊，簡化后的Voronoi圖網絡如圖3(b)所示，每個節(jié)點最多歸4個多邊形共有，每條邊最多歸2個多邊形共有。

圖3 測區(qū)Voronoi圖鑲嵌線網絡局部

2.3 鑲嵌線網絡中節(jié)點的選擇優(yōu)化

要想使鑲嵌線自動繞過建筑物成像區(qū)域，必須先保證網絡中節(jié)點不落在建筑物的成像區(qū)域內。由2.1節(jié)可知，每個節(jié)點最多歸4個多邊形共用，因此要保證節(jié)點不在4張影像中任一張影像上建筑物的成像區(qū)域內。

由于成像位置和姿態(tài)的不同，一棟建筑物在每張正射影像的成像區(qū)域也不相同，如圖4所示，O1，O2是相鄰兩張影像的攝影中心，棱柱體ABCD-abcd為建筑物，則該建筑物在正射影像O2上成像區(qū)域為多邊形A′B′C′cba，而在正射影像O1上成像區(qū)域為多邊形B″C″D″dab。則該建筑物在2張影像的成像總區(qū)域為多邊形A′B′EB″C″D″da[24]。

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

以下給出了鑲嵌線網絡所有節(jié)點的優(yōu)化過程：

①判斷節(jié)點P是否歸2張以上影像有效多邊形共有(以下以4張共有為例)，如是進行步驟②，否則進行步驟⑦。

②以節(jié)點P為中心，搜索方圓一定距離范圍的所有建筑物，記為B1，B2，…，Bn。

③對于每個建筑物，計算其在這4張DOM上的成像區(qū)域，并求并集作為該建筑物的整體城像區(qū)域，記為S1，S2，…，Sn。

④對求解出的S1，S2，…，Sn進行求交判斷，如果有交集則合并同時記下Si與Bi的對應關系，直到剩下S1，S2，…，Sm(m≤n)的任意二者都不相交。

⑤逐一判斷節(jié)點P是否在S1，S2，…，Sm(m≤n)內部，如果在則記下該房屋的編號Bi，同時找到Bi對應的合并后Si，如果都不在，則該節(jié)點無需優(yōu)化進入步驟⑦。

⑥求出包含多邊形Si的最小外包矩形，然后找出距離節(jié)點P距離最小的一個矩形端點Q作為節(jié)點P的優(yōu)化后位置。然后將鑲嵌線網絡中所有的節(jié)點P更新為點Q。

⑦判斷P是否為最后一個節(jié)點，如果是則優(yōu)化完畢，如果不是則進行下一個節(jié)點的檢測更新，進入步驟①。

2.4 鑲嵌網絡中鑲嵌線的選擇優(yōu)化

每條鑲嵌線由2個影像有效多邊形共有，因此對于鑲嵌線自動繞過建筑物的優(yōu)化需要同時考慮2張影像共有的建筑物成像情況以及最優(yōu)化路徑選取。

以下給出了單條鑲嵌線優(yōu)化的方法步驟：

①判斷鑲嵌線L是否為2個影像有效多邊形共有，若僅為一個多邊形擁有，則該鑲嵌線無需優(yōu)化，直接進入步驟⑨。

②根據鑲嵌線L與影像的對應關系找到該鑲嵌線兩側的影像序號Pici，Picj，并根據建筑物與影像間的對應關系，將同時出現在該2張影像上的建筑物挑選出來記為B1，B2，…，Bn。

③對于每個建筑物，計算其在Pici和Picj2張正射影像上的成像區(qū)域，并求并集作為該建筑物的整體城像區(qū)域，記為S1，S2，…，Sn。

④判斷每個建筑物成像區(qū)域Si是否與鑲嵌線L相交，將相交的多邊形留下記為：S1，S2，…，Sm(m≤n)，若都不相交則該鑲嵌線不需要優(yōu)化處理，進入步驟⑨。

⑤對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張影像多邊形中的該條鑲嵌線。

⑨進入下一條鑲嵌線的優(yōu)化，直到所有的鑲嵌線優(yōu)化完畢，得到整個測區(qū)自動繞過建筑物成像區(qū)域的鑲嵌線網絡。

圖5 鑲嵌線繞過建筑物優(yōu)化示意圖

3 實驗與分析

為驗證本文方法的有效性，以某地4個航帶共32張由SWDC-5航攝儀拍攝的分辨率為0.1 m、航向重疊約為60%，旁向重疊約為45%、框幅為8 230×6 168的航攝影像及測區(qū)分辨率為2 m的DEM和656個建筑物房頂矢量進行實驗。

3.1 建筑物成像區(qū)域計算結果

圖6給出了利用本文方法從建筑物房頂矢量獲得其在單張DOM上成像區(qū)域的過程和結果，其中圖6(a)為測區(qū)一棟建筑物的房頂矢量，圖6(b)為將其投影到DEM上得到的建筑物簡易三維模型，圖6(c)為由建筑物DBM投影到原始影像上的成像范圍，圖6(d)為基于DEM迭代逼近得到的建筑物在單片DOM上成像區(qū)域的矢量范圍。從圖中可以看出，利用本文方法可以由建筑物房頂矢量準確獲得建筑物在單片DOM上的成像范圍。

圖6 獲取建筑物成像區(qū)域物方矢量過程

3.2 鑲嵌線網絡自動選擇分析

圖7給出了鑲嵌線網絡中節(jié)點優(yōu)化的結果，其中圖7(a)～圖7(d)是同一棟建筑物在節(jié)點相關4張單片DOM上的成像區(qū)域(黑色多邊形)。圖7(d)中外圍多邊形為該建筑在4張影像的總成像區(qū)域多邊形，紅圈內為節(jié)點優(yōu)化前的位置，圖7(e)為節(jié)點優(yōu)化后的位置。從圖中可以看出，利用本文的方法成功將落在建筑物成像區(qū)域內的網絡節(jié)點檢測并移出，為鑲嵌線的優(yōu)化奠定了基礎。

經過節(jié)點優(yōu)化后的鑲嵌線網絡中每條鑲嵌線最多歸2個多邊形共有，因此鑲嵌線的優(yōu)化只需考慮鑲嵌線穿過的房屋在對應的2張影像的成像區(qū)域的并集即可。如圖8所示，一條鑲嵌線連續(xù)穿過了3個建筑物的成像區(qū)域，圖8(a)為鑲嵌線與鑲嵌線對應上方影像建筑物成像區(qū)域的位置關系，圖8(b)為鑲嵌線與對應下方影像上建筑物成像區(qū)域的位置關系，圖8(c)為原始鑲嵌線鑲嵌結果，圖8(d)是優(yōu)化后的鑲嵌線鑲嵌結果。比較圖8(c)和圖8(d)可知，利用本文方法優(yōu)化后的鑲嵌線成功繞過了建筑物成像區(qū)域，避免了鑲嵌線兩側紋理錯位的發(fā)生。

圖7 網絡節(jié)點優(yōu)化過程及結果

圖8 鑲嵌線優(yōu)化過程與結果

圖9給出了測區(qū)整體鑲嵌線網絡優(yōu)化前后的套合圖與局部對比圖，可以看出利用本文方法優(yōu)化選擇后的鑲嵌線不僅繞過了建筑物的成像區(qū)域保證了鑲嵌結果的紋理質量，同時也最大程度上保留了Voronoi圖鑲嵌線網絡局部投影變形理論最小的特點，保證了鑲嵌結果的幾何精度。

圖9 測區(qū)鑲嵌線網絡優(yōu)化前后對比

3.3 與其他方法對比分析

本文利用Inpho OrthoVista 和IPS 4.2(icaros photogrammetric suite)的Stitch 模塊分別進行了鑲嵌線網絡的自動選擇，圖10給出了3種方法鑲嵌線網絡的局部，表1給出了3種方法以該局部區(qū)域中心影像有效區(qū)域多邊形(與周邊影像鑲嵌線組成的多邊形)為對象的量化值統計情況。

圖10 與其他鑲嵌線網絡自動選擇方法結果對比

表1 與其他方法對比

首先，在鑲嵌線網絡選擇效果上，從圖10可以看出OrthoVista 自動選擇的鑲嵌線有5處穿過了建筑物，IPS 自動選擇的鑲嵌線有4處穿過建筑物(該片區(qū)共有6處)，造成了鑲嵌線兩側紋理錯位，而本文方法該片區(qū)鑲嵌線網絡繞過了所有的建筑物成像區(qū)域。特別是該片區(qū)右上方建筑物間距極小、密集程度高的情況，大多數自動選擇方法無法繞過，而本文方法在鑲嵌線選擇優(yōu)化過程中通過合并有交集的多個建筑物成像區(qū)域的方法，有效避開了建筑物密集區(qū)域。

其次，在測區(qū)鑲嵌線網絡的整體布局上，本文方法是對Voronoi圖進行優(yōu)化選擇，最大程度保留了其鑲嵌結果投影變形理論最小的優(yōu)點。實驗中通過測量建筑物同一房頂點與對應基底點的距離獲得其投影變形大小，從表1可以看出，利用本文方法鑲嵌結果中建筑物B1的投影變形最小為1.99 m，而OrthoVista鑲嵌結果投影變形最大達到8.24 m。對于建筑物B2，本文方法和OrthoVista鑲嵌結果中建筑物B2的紋理均取自同一影像，因此投影變形大小一致為1.86 m，而IPS Stich的鑲嵌結果則取自于另一張影像，投影變形為3.47 m。

同時，本文算法基于矢量計算得到的鑲嵌線網絡，而另外2種算法是構建異差影像在像素級尋找最優(yōu)路徑的方法獲得鑲嵌線，因此本文方法算法更為簡單，同時獲得的鑲嵌線網絡為面結構，且節(jié)點數目更少。如表1所示，利用本文方法得到中心影像的有效范圍矢量只有347個節(jié)點，遠遠低于前2種方法的1 370和735個節(jié)點。

最后，圖11給出了3種方法自動選擇生成整個測區(qū)鑲嵌線網絡的耗時，其中前2種方法是以單片正射影像為基礎數據，對相鄰影像進行全幅鑲嵌線自動搜索，得到航向和旁向2個方向上相鄰影像的鑲嵌線，再進行鑲嵌線相互切割得到最終鑲嵌線網絡。本文方法以建筑物矢量及高度、測區(qū)DEM、影像定向參數為基礎數據，先生成Voronoi圖初始鑲嵌線網絡并簡化，再根據建筑物成像區(qū)域進行鑲嵌線網絡的選擇優(yōu)化。OrthoVista軟件耗時1 925 s，IPS軟件Stitch模塊由于利用了GPU加速技術耗時為963 s，而本文方法由于計算復雜度遠遠低于前2種方法，且多為矢量運算，從數據讀入到鑲嵌線網絡選擇優(yōu)化完成共耗時349 s。

圖11 耗時對比

4 結束語

本文提出了一種建筑物矢量輔助的城市大比例尺正射影像鑲嵌線網絡自動選擇方法。首先得到建筑物在單片正射影像上的成像區(qū)域，然后自動檢測并優(yōu)化初始Voronoi圖鑲嵌線網絡中的所有節(jié)點和鑲嵌線，得到整個測區(qū)繞開建筑物成像區(qū)域的最優(yōu)鑲嵌線網絡。該方法不僅能夠快速高效地得到繞過建筑物的鑲嵌線網絡，并且保留了Voronoi圖網絡鑲嵌結果投影變形理論最小的特點，保證了鑲嵌后DOM影像的質量，為城市大比例尺正射影像鑲嵌線網絡選擇提供了一種方法。下一步將研究根據左右差值影像對鑲嵌線網絡進行優(yōu)化處理，減少對建筑物房頂矢量的依賴。