一種新的生成ADS正射影像方法

曾 宇，蘇曉斌

(湖南省第二測繪院，湖南 長沙 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.