基于幾何約束移動(dòng)最小二乘的TomoSAR山區(qū)點(diǎn)云高精度三維重建方法

李曉婉 梁興東* 張福博 劉云龍 李焱磊 郭其昌 萬(wàn)陽(yáng)良 卜祥璽

①(中國(guó)科學(xué)院空天信息創(chuàng)新研究院微波成像技術(shù)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室 北京 100190)

②(中國(guó)科學(xué)院大學(xué)電子電氣與通信工程學(xué)院 北京 100049)

1 引言

通過(guò)對(duì)同一觀(guān)測(cè)場(chǎng)景的多角度觀(guān)測(cè)，層析合成孔徑雷達(dá)(Tomographic Synthetic Aperture Radar,TomoSAR)可以將傳統(tǒng)的二維SAR圖像拓展到三維圖像，已然成為地形測(cè)繪、城市監(jiān)督、生物量估計(jì)等三維應(yīng)用領(lǐng)域中非常有效的數(shù)據(jù)獲取手段[1—3]。

作為典型的自然場(chǎng)景，山區(qū)三維重建在地形測(cè)繪、資源勘探、災(zāi)害評(píng)估等方面具有重要意義。然而，TomoSAR山區(qū)三維重建主要面臨著如下難題：(1)表面覆蓋有大量分布式目標(biāo)，對(duì)時(shí)間去相關(guān)非常敏感[4]；(2)高程跨度大，點(diǎn)云高程模糊問(wèn)題顯著[5]；(3)由于有限的基線(xiàn)長(zhǎng)度和通道數(shù)量，點(diǎn)云高程向定位精度非常差[6,7]。通過(guò)在跨航向布設(shè)天線(xiàn)陣列，機(jī)載陣列TomoSAR可以有效地提高圖像間的時(shí)間相關(guān)性[8,9]，為山區(qū)三維重建提供數(shù)據(jù)支持，解決了第1個(gè)難題?；诜指詈椭亟M的高程模糊求解方法可以恢復(fù)出無(wú)高程模糊的點(diǎn)云[10]，解決了第2個(gè)難題。目前針對(duì)第3個(gè)難題研究較少，迫切地需要提出一種可以提高點(diǎn)云高程向定位精度的TomoSAR山區(qū)三維重建方法。

現(xiàn)有的三維重建方法主要是針對(duì)TomoSAR和LiDAR提出的。前者多適用于符合平面幾何約束模型的建筑物，如RANSAC算法[11]、TomoSeed算法[12]，以及一系列基于建筑物足跡和機(jī)器學(xué)習(xí)的方法[13—15]。顯然，這些方法難以用于曲面結(jié)構(gòu)的山區(qū)。后者的處理對(duì)象較為豐富，常利用法向量和曲率等幾何特征、點(diǎn)云分布特征以及拉普拉斯正則化處理等方法進(jìn)行曲面重建[16—20]。TomoSAR點(diǎn)云密度和LiDAR點(diǎn)云密度相比擬[21]。然而，相比于LiDAR點(diǎn)云，有限基線(xiàn)數(shù)量、相參成像和電磁散射特性等因素使得TomoSAR點(diǎn)云的定位精度較差且可能存在離群點(diǎn)和虛假目標(biāo)等干擾；側(cè)視幾何使得TomoSAR點(diǎn)云存在固有的幾何畸變問(wèn)題。因此，LiDAR點(diǎn)云處理方法多難以直接用于TomoSAR山區(qū)三維重建[7]。最小二乘(Least-Squares,LS)法在曲線(xiàn)(曲面)散亂點(diǎn)云去噪和曲面擬合等方面得到了廣泛的應(yīng)用[22]。然而，LS追求全局最優(yōu)，對(duì)于復(fù)雜點(diǎn)云容易出現(xiàn)過(guò)平滑問(wèn)題。移動(dòng)最小二乘(Moving Least-Squares,MLS)法利用分段擬合的思想，可以在一定程度上改善過(guò)平滑問(wèn)題[23]。盡管如此，MLS難以實(shí)現(xiàn)高精度的TomoSAR山區(qū)三維重建。一個(gè)原因是坐標(biāo)系表征問(wèn)題。雷達(dá)坐標(biāo)系允許只進(jìn)行高程向定位誤差修正，但疊掩使其不滿(mǎn)足MLS函數(shù)逼近條件。大地坐標(biāo)系滿(mǎn)足此條件，但是山區(qū)地形起伏多變，這種對(duì)高度向定位誤差修正的方式卻可能引入更強(qiáng)烈的斜距誤差。此外，MLS缺乏對(duì)重建結(jié)果的幾何約束，難以對(duì)上述誤差進(jìn)行再次修正。

本文充分地利用了TomoSAR點(diǎn)云高程與地距間的單調(diào)遞增關(guān)系，提出了一種基于幾何約束MLS的TomoSAR山區(qū)點(diǎn)云三維重建方法。通過(guò)投影變化和幾何約束MLS兩個(gè)主要步驟，可以針對(duì)性地抑制高程向定位誤差。本文的主要?jiǎng)?chuàng)新工作如下：(1)提出了一種新的坐標(biāo)系投影變換方法，這是利用幾何約束進(jìn)行高程定位誤差修正的前提；(2)將TomoSAR點(diǎn)云的幾何約束關(guān)系融入MLS，保證了重建結(jié)果的全局約束性；(3)提出了一種幾何約束MLS迭代求解算法，保證了結(jié)果的準(zhǔn)確性。

2 TomoSAR點(diǎn)云幾何約束模型

TomoSAR通過(guò)多角度觀(guān)測(cè)實(shí)現(xiàn)了疊掩地形的高程分辨。圖1示意了TomoSAR山區(qū)觀(guān)測(cè)幾何，黑色原點(diǎn)示意陣列天線(xiàn)，紅色加粗曲線(xiàn)示意被照射山區(qū)地形，黑色加粗曲線(xiàn)示意被遮擋山區(qū)地形，H對(duì)應(yīng)載機(jī)高度，P1,P2和 P3疊掩在一起。假設(shè)基線(xiàn)均勻分布，最左側(cè)通道為主通道，對(duì)于觀(guān)測(cè)目標(biāo) P2，第n個(gè)通道與主通道間的相位差可以表示為[24]

圖1 TomoSAR山區(qū)觀(guān)測(cè)幾何示意圖(y和z分別表示地距向和高度向，r和s分別表示斜距向和高程向)Fig.1 Diagram of TomoSAR mountain observation geometry (y and z indicate the ground direction and the height direction,respectively,r and s indicate the range direction and the elevation direction,respectively)

其中，λ是波長(zhǎng)，Δb是相鄰?fù)ǖ篱g的基線(xiàn)長(zhǎng)度，θ是被觀(guān)測(cè)目標(biāo)相對(duì)于主通道的入射角，β是基線(xiàn)水平傾角。

根據(jù)壓縮感知高程分辨求解模型，第n個(gè)通道與主通道間的相位差也可以表示為[25]

其中，Nh是一個(gè)高程周期的離散點(diǎn)數(shù)，h是高程向譜估計(jì)值。

令ρs表示高程向超分辨率，h與真實(shí)高程值s間的關(guān)系為s=hρs。其中，真實(shí)高程值為目標(biāo)的高程向距離，即地距-高度平面原點(diǎn)o相距目標(biāo)所在方位-斜距平面的高程向距離。以觀(guān)測(cè)目標(biāo) P2為例，其斜距值和真實(shí)高程值分別為r2和s2。為了簡(jiǎn)便，下述內(nèi)容均將h稱(chēng)為高程值。

聯(lián)立式(1)和式(2)，可以得出

簡(jiǎn)化常數(shù)項(xiàng)，高程值可以進(jìn)一步表示為

其中，常數(shù)項(xiàng)k=λ/2ΔbNh。

通常來(lái)說(shuō)，0＜β ＜π/2，下視角0＜θ ＜π/2。在這種情況下，目標(biāo)高程值與下視角具有單調(diào)遞增關(guān)系。同時(shí)，照射幾何表明，只有下視角隨地距增加而增加的點(diǎn)才會(huì)被照射到。如圖1所示，P1,P2和 P3滿(mǎn)足照射幾何，均可以被照射到，其高程值隨地距增加而增加。盡管 P4和 P5對(duì)應(yīng)地距范圍內(nèi)黑色加粗地形不會(huì)被照射到，但其在TomoSAR點(diǎn)云中表現(xiàn)為兩端高程連續(xù)的空洞區(qū)域，不會(huì)破壞這種單調(diào)遞增關(guān)系。這是因?yàn)?P4和 P5的下視角相同，由式(4)可以推得其高程值也相同。

綜上，TomoSAR點(diǎn)云高程與地距具有單調(diào)遞增幾何約束關(guān)系。

3 幾何約束移動(dòng)最小二乘方法

本文方法的流程圖見(jiàn)圖2，主要包括投影變換和幾何約束MLS這兩個(gè)步驟。下面對(duì)其進(jìn)行詳細(xì)介紹。

圖2 幾何約束MLS方法流程圖Fig.2 Flowchart of the geometry constrained MLS method

3.1 投影變換

投影變換主要包括駐點(diǎn)定位和地距投影兩步。前者用于估計(jì)原始散射點(diǎn)的地距順序，后者用于按照已知地距順序?qū)υ忌⑸潼c(diǎn)進(jìn)行地距重排，從而避免地距倒置帶來(lái)的約束失效問(wèn)題。

3.1.1 駐點(diǎn)定位

圖3示意了常見(jiàn)的S形山區(qū)疊掩地形。A和B稱(chēng)為地形駐點(diǎn)，其具有位于左向開(kāi)口或右向開(kāi)口的U形曲線(xiàn)斜距邊緣的特性。以地形駐點(diǎn)高程為門(mén)限可以將地形分為I區(qū)，II區(qū)和III區(qū)3個(gè)部分。I區(qū)和III區(qū)內(nèi)點(diǎn)云呈現(xiàn)出正斜率走勢(shì)，真實(shí)地距隨斜距增加而增加；II區(qū)內(nèi)點(diǎn)云呈現(xiàn)出負(fù)斜率走勢(shì)，真實(shí)地距隨斜距減少而增加。由此可見(jiàn)，在正確定位駐點(diǎn)后，可以利用相鄰駐點(diǎn)高程區(qū)間內(nèi)散射點(diǎn)的斜率特性估計(jì)出原始散射點(diǎn)的相對(duì)地距關(guān)系。

圖3 斜距-高程平面S形山區(qū)疊掩地形結(jié)構(gòu)圖Fig.3 The S shaped layover mountainous terrain on the range-elevation plane

山區(qū)地形復(fù)雜多變，且存在地形遮擋等問(wèn)題。TomoSAR山區(qū)點(diǎn)云還可能呈現(xiàn)出如圖4所示的更復(fù)雜形狀，但不影響上述規(guī)律。為此，本文利用二分遞歸思想進(jìn)行分區(qū)駐點(diǎn)定位，具體步驟如下：

步驟1 利用同一方位-斜距單元內(nèi)疊掩目標(biāo)的高程跨度遠(yuǎn)大于非疊掩目標(biāo)高程跨度的特性執(zhí)行疊掩地形判斷。具體地說(shuō)，只要檢測(cè)地形中高程跨度大于閾值Th的點(diǎn)數(shù)大于閾值Tnum，就將其視為疊掩地形，并進(jìn)入步驟2。反之，將其視為非疊掩地形，進(jìn)入步驟3；

步驟2 利用前述駐點(diǎn)特性執(zhí)行駐點(diǎn)檢測(cè)。首先檢測(cè)斜距邊緣點(diǎn)。為了避免上下邊界點(diǎn)的影響，只有與輸入數(shù)據(jù)邊界點(diǎn)的斜距向距離和高程向距離大于閾值Trdis1和閾值Thdis1的點(diǎn)會(huì)被保留。最后，會(huì)進(jìn)一步判斷邊緣點(diǎn)是否符合U形曲線(xiàn)特征。對(duì)于右向開(kāi)口曲線(xiàn)，位于邊緣點(diǎn)左側(cè)Trdis2和Thdis2矩形框內(nèi)的鄰域點(diǎn)數(shù)需小于閾值Tnei，同時(shí)，還需要輸入數(shù)據(jù)中存在位于其右側(cè)高程兩側(cè)且與其斜距向距離大于Tr的點(diǎn)。左向開(kāi)口曲線(xiàn)判斷準(zhǔn)則與之相反。如圖4(a)，點(diǎn)A和點(diǎn)T為斜距邊緣點(diǎn)。但是，點(diǎn)T為邊界點(diǎn)會(huì)被舍棄。又因?yàn)辄c(diǎn)A符合右向開(kāi)口U形曲線(xiàn)特征，會(huì)被檢測(cè)為駐點(diǎn)。

步驟3 下述內(nèi)容均以dataij表示以pi和pj為上下邊界的待檢測(cè)地形，fact(dataij) 表示對(duì)待檢測(cè)地形dataij進(jìn)行駐點(diǎn)定位。如圖4(a)，data12則對(duì)應(yīng)以p1和p2為上下邊界的待檢測(cè)地形，fact(data12)對(duì)應(yīng)主函數(shù)。下面分3種情況執(zhí)行二分遞歸。

(1) 情況1為非疊掩地形，輸出結(jié)果為空且會(huì)結(jié)束駐點(diǎn)定位；

(2) 情況2為疊掩地形且已經(jīng)檢測(cè)到駐點(diǎn)的情況。這時(shí)會(huì)對(duì)兩個(gè)邊界點(diǎn)和駐點(diǎn)進(jìn)行高程排序，并對(duì)相鄰點(diǎn)內(nèi)的地形進(jìn)行遞歸。如圖4(a)所示，對(duì)駐點(diǎn)A (p3)和兩個(gè)邊界點(diǎn)排序后，相鄰點(diǎn)內(nèi)的地形為data13和data23，輸出結(jié)果為p3∪f(wàn)act(data13)∪f(wàn)act(data23)；

(3) 情況3為疊掩地形但沒(méi)有檢測(cè)到駐點(diǎn)的情況。如圖4(b)，p3為高程向區(qū)間二分后確定的分割點(diǎn)，使用二分法遞歸執(zhí)行主函數(shù)。不考慮p3恰好為駐點(diǎn)的情況時(shí)，輸出結(jié)果為 fact(data13)∪f(wàn)act(data23)。顯然，這可能會(huì)漏檢恰好同為駐點(diǎn)和分割點(diǎn)的點(diǎn)p3。為了避免這種情況的發(fā)生，會(huì)在位于p3上方的駐點(diǎn)和邊界點(diǎn)中找到與其高程距離最小的點(diǎn)p4，同理，會(huì)在p3下方找到p5，并再次執(zhí)行f act(data45)，輸出結(jié)果為fact (data13)∪f(wàn)act (data23)∪f(wàn)act(data45)。

圖4 兩種典型地形駐點(diǎn)定位示意Fig.4 Diagram of two typical terrain for stagnation point positioning

綜上，駐點(diǎn)定位主函數(shù)f act(data)為

3.1.2 地距投影

步驟1 方位向循環(huán)：

(1) 以駐點(diǎn)高程升序的方式估計(jì)所有高程向相鄰駐點(diǎn)間點(diǎn)的斜率；

(2) 對(duì)于正向斜率區(qū)域，以斜距升序方式調(diào)整點(diǎn)的順序；對(duì)于負(fù)向斜率區(qū)域，以斜距降序方式調(diào)整點(diǎn)的順序；

(4) 根據(jù)式(3)計(jì)算排序后點(diǎn)的地距：yn=rnsinθn；

(5) 保持每個(gè)散射點(diǎn)的方位和高程不變，投影坐標(biāo)系表征的點(diǎn)為

步驟2 結(jié)束循環(huán)；

步驟3 輸出投影變換后的點(diǎn)：P′={pi(ai,yi,hi)|ai ∈(A1,A2,...,AJ),i=1,2,...,Ntotal}及其斜距r={ri|i=1,2,...,Ntotal}。

圖5說(shuō)明了地距倒置導(dǎo)致幾何約束失效的問(wèn)題，這也是進(jìn)行地距投影的主要原因。圖5(a)示意斜距-高程平面的散射點(diǎn)，圖5(b)示意地距-高程平面內(nèi)的散射點(diǎn)。可以看出，受高程向定位誤差影響，點(diǎn)地距倒置，而且這使得幾何約束失效。通過(guò)駐點(diǎn)定位和斜率估計(jì)可以推斷出點(diǎn)位于II區(qū)，所以由可以推斷出這兩個(gè)點(diǎn)的真實(shí)地距順序。點(diǎn)示意了按照真實(shí)地距排序的點(diǎn)可以看出，點(diǎn)不再滿(mǎn)足單調(diào)遞增約束關(guān)系，幾何約束MLS可以對(duì)其進(jìn)行再次修正。

圖5 投影變換示意Fig.5 Diagram of projection transform

3.2 幾何約束MLS

3.2.1 函數(shù)模型

參考MLS局部加權(quán)的思想，任意點(diǎn)pn,n=1,2,...,N，局部子空間內(nèi)幾何約束MLS代價(jià)函數(shù)為

其中，第1項(xiàng)為MLS懲罰項(xiàng)，w(pi-pn)是點(diǎn)pn處具有緊支撐特性的光滑連續(xù)權(quán)函數(shù)，f(pi)是局部子空間的擬合函數(shù)；第2項(xiàng)為單調(diào)遞增懲罰項(xiàng)，單調(diào)遞增懲罰因子κ大 于0，κ越大，單調(diào)遞增約束越緊，vn是單調(diào)遞增懲罰局部子空間的權(quán)函數(shù)。

令影響域半徑r=‖pi-pn‖/dmI，高斯權(quán)函數(shù)w(pi-pn)[26]為

其中，dmI 決定局部子空間的大小，β決定權(quán)函數(shù)的衰減速度。

擬合函數(shù)可以表示為

其中，多項(xiàng)式函數(shù)矩陣D(p)的展開(kāi)式為

待求系數(shù)c=[c1c2...cM]T,d(x,y)=[d1d2...dM]T是最高次數(shù)為M的多項(xiàng)式基函數(shù)。

權(quán)函數(shù)v=[v1v2...vN]T需要保證只有在擬合結(jié)果違背單調(diào)遞增約束時(shí)，幾何約束懲罰項(xiàng)才起作用。任意點(diǎn)pn相 對(duì)其鄰域子空間p win的單調(diào)遞增權(quán)函數(shù)可以表示為

其中，zneiL和zneiR分別是與點(diǎn)pn方位向相同且位于其左右兩側(cè)的p win個(gè)鄰域點(diǎn)的第三維信息。

3.2.2 問(wèn)題求解

幾何約束MLS的解析解可以表示為

權(quán)函數(shù)矩陣W是對(duì)角元素為w(pi-pn)的對(duì)角陣，單調(diào)遞增權(quán)函數(shù)V是對(duì)角元素為vi的單位陣。

然而，下述幾點(diǎn)不足使得符合單調(diào)遞增約束的解析解難以得到：(1)矩陣V與解析解耦合。也就是說(shuō)，需要先給定多個(gè)解析解，才能求出對(duì)應(yīng)的矩陣V以及代價(jià)函數(shù)，最后再進(jìn)行最小代價(jià)函數(shù)查找才能確定最優(yōu)解，計(jì)算代價(jià)很大；(2)每個(gè)局部子空間內(nèi)點(diǎn)的分布不一致，對(duì)應(yīng)MLS懲罰項(xiàng)擬合誤差也不一致。所以，最好根據(jù)子空間數(shù)據(jù)分布設(shè)置N個(gè)單調(diào)遞增懲罰因子κ。顯然，這是比較困難的；(3)對(duì)于表現(xiàn)為下降趨勢(shì)的輸入數(shù)據(jù)，無(wú)法滿(mǎn)足單調(diào)遞增約束。

于是，本文給出了一種迭代求解方法。針對(duì)問(wèn)題(1)，該方法采取先執(zhí)行MLS后針對(duì)性地對(duì)非單調(diào)遞增點(diǎn)進(jìn)行高程修正的策略。這是因?yàn)楫?dāng)單調(diào)遞增懲罰項(xiàng)為0時(shí)，MLS結(jié)果就是幾何約束MLS結(jié)果；針對(duì)問(wèn)題(2)，引入不小于0的局部抖動(dòng)因子dh，這可以將擬合結(jié)果的高程抖動(dòng)限定在dh內(nèi)；針對(duì)問(wèn)題(3)，迭代修正可以重復(fù)修正不滿(mǎn)足單調(diào)遞增約束的點(diǎn)。整個(gè)幾何約束MLS求解算法見(jiàn)表1。對(duì)于MLS求解，本文采用不斷放大dmI的策略來(lái)自適應(yīng)確定影響域半徑，直到其滿(mǎn)足MLS求解的非病態(tài)條件。此外，與格網(wǎng)化的MLS求解流程[27]不同的是，為了保證斜距不變性，本文會(huì)先求所有點(diǎn)MLS擬合結(jié)果。最后，通過(guò)大地坐標(biāo)系投影變換和均勻格網(wǎng)化處理[27]得到重建結(jié)果。

表1 幾何約束MLS求解算法Tab.1 The solution algorithm of the geometry constrained MLS

4 實(shí)驗(yàn)結(jié)果分析

為了驗(yàn)證所提方法的有效性，給出了仿真數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)的處理結(jié)果。作為對(duì)比，給出了傳統(tǒng)的大地坐標(biāo)系下執(zhí)行的LS[22]和MLS[27]處理結(jié)果。

4.1 參數(shù)設(shè)置討論

投影變換參數(shù)對(duì)數(shù)據(jù)敏感度較高，需要根據(jù)實(shí)際數(shù)據(jù)進(jìn)行經(jīng)驗(yàn)選取，以保證駐點(diǎn)的U形分布特性。幾何約束MLS參數(shù)對(duì)數(shù)據(jù)敏感度較低。參數(shù)pwin值越大，單調(diào)遞增權(quán)函數(shù)vn為1的概率越高，會(huì)導(dǎo)致迭代求解的效率變低?？紤]到本文以地距升序局部單調(diào)遞增修正的方式來(lái)達(dá)到全局單調(diào)遞增，將其設(shè)為 pwin=2。當(dāng)對(duì)同一點(diǎn)相鄰兩次修正的效果de小于閾值Te時(shí)，會(huì)進(jìn)入到下一點(diǎn)，文中將其設(shè)為T(mén)e=1e-5 。抖動(dòng)因子 dh越小，單調(diào)遞增約束越緊，反之，單調(diào)遞增約束越松，允許結(jié)果存在 dh范圍內(nèi)的高程抖動(dòng)。因此，在數(shù)據(jù)質(zhì)量好且單調(diào)遞增約束緊時(shí)，其最小值為0。否則，可以對(duì)其進(jìn)行適當(dāng)放大。此外，對(duì)于MLS求解，本文選用二次基和高斯權(quán)函數(shù)。dmI初始化為原始點(diǎn)云地距分辨率ρy的3倍，且以2倍的乘因子進(jìn)行放大。參數(shù)β對(duì)重建結(jié)果也會(huì)有一定程度的影響，可以根據(jù)不同數(shù)據(jù)進(jìn)行合理設(shè)置，文中將其設(shè)為3。

綜上，本文中幾何約束MLS求解參數(shù)設(shè)置為：pwin=2,T e=1e-5,d mI=3ρy,β=3。其他參數(shù)就要根據(jù)數(shù)據(jù)進(jìn)行經(jīng)驗(yàn)設(shè)置。作為對(duì)比，LS多項(xiàng)式擬合階數(shù)設(shè)置為5，MLS的參數(shù)同本文方法。此外，3種方法的格網(wǎng)化節(jié)點(diǎn)一致，均為5 m。

4.2 仿真實(shí)驗(yàn)

4.2.1 數(shù)據(jù)生成

圖6給出了LS,MLS和本文方法對(duì)于4種不同地形的仿真數(shù)據(jù)處理結(jié)果。數(shù)據(jù)1為連續(xù)的S形地形，數(shù)據(jù)2為底部被遮擋住的地形，數(shù)據(jù)3和數(shù)據(jù)4示意包含一個(gè)和多個(gè)S形地形的復(fù)雜地形。數(shù)據(jù)1由3次樣條曲線(xiàn)生成，數(shù)據(jù)2到數(shù)據(jù)4地形更為復(fù)雜，是以AW3D30 DSM數(shù)據(jù)[28]為輸入的POV-RAY軟件[29]輸出數(shù)據(jù)。其中，POV-RAY軟件用于模擬地形遮擋情況。最后，通過(guò)在高程向上加入均值為0，標(biāo)準(zhǔn)差為2.5個(gè)高程向像素單元的高斯白噪聲模擬高程向定位誤差。仿真系統(tǒng)參數(shù)同下述實(shí)測(cè)實(shí)驗(yàn)。對(duì)于場(chǎng)景中心，地距分辨率約為0.5 m。

圖6 LS,MLS和本文方法的仿真結(jié)果Fig.6 Simulation results of LS,MLS and the proposed method

4.2.2 結(jié)果與分析

數(shù)據(jù)1和數(shù)據(jù)2的投影變換參數(shù)設(shè)置為：T h=10,Tnum=10,T rdis1=100,T hdis1=5,T rdis2=5,Thdis2=2,T nei=5,T r=20。數(shù)據(jù)3的不同參數(shù)為：Trdis1=10,T r=48。數(shù)據(jù)4的不同參數(shù)為：Trdis1=10,T r=12。此外，數(shù)據(jù)1的抖動(dòng)因子d h=0.3，其他數(shù)據(jù)為d h=1.2。

下面定性分析3種方法的重建性能?？梢钥闯?，全局多項(xiàng)式擬合的LS過(guò)平滑問(wèn)題非常突出，MLS可以對(duì)其進(jìn)行一定程度的改進(jìn)，但是在駐點(diǎn)鄰近區(qū)域內(nèi)會(huì)引入斜距誤差。本文方法只進(jìn)行高程向定位誤差修正，可以很好地避免引入斜距誤差，最接近于真值，尤其在駐點(diǎn)鄰近區(qū)域。

然后，以絕對(duì)測(cè)高精度為評(píng)估指標(biāo)[30]，給出定量化評(píng)估結(jié)果。3種方法處理結(jié)果的格網(wǎng)節(jié)點(diǎn)一致，且有與其一一對(duì)應(yīng)的真值。原始點(diǎn)云非均勻分布，利用線(xiàn)性插值進(jìn)行格網(wǎng)化。因此，絕對(duì)測(cè)高精度為重建結(jié)果與其對(duì)應(yīng)真值的高度向均方根誤差，下述稱(chēng)為精度。表2列出了4種數(shù)據(jù)的精度?？梢钥闯?，LS和MLS在追求平滑的同時(shí)，可能降低原始點(diǎn)云的精度，如數(shù)據(jù)2；MLS可以在一定程度上提高LS的重建精度，尤其是對(duì)于數(shù)據(jù)3和數(shù)據(jù)4這種復(fù)雜地形；本文方法均可以將MLS的定位誤差降低至原來(lái)的1/2左右。進(jìn)一步地，圖7呈現(xiàn)了不同方法的精度對(duì)比曲線(xiàn)。輸入為加入不同高程向誤差的數(shù)據(jù)1。以高程向像素?cái)?shù)為單位，10個(gè)切片的高斯白噪聲的標(biāo)準(zhǔn)差依次為：0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,均值為0?？梢钥闯鰧?duì)于同一切片，所提方法性能最優(yōu)；隨噪聲增加，3種方法的定位效果均下降，但是所提方法性能下降緩慢，且逐漸顯著優(yōu)于MLS。

圖7 基于數(shù)據(jù)1的不同方法的精度對(duì)比曲線(xiàn)Fig.7 The precision comparison curves of different methods based on data 1

表2 仿真結(jié)果精度(m)Tab.2 Precision of the simulated results (m)

最后，圖8呈現(xiàn)了本文方法對(duì)于切片5的主要中間結(jié)果。圖8(a)中原始點(diǎn)云被賦予顏色圖，值越接近1，由投影變換預(yù)測(cè)得到的地距值越大?？梢钥闯觯嬷档母叱毯偷鼐酀M(mǎn)足單調(diào)遞增關(guān)系，原始點(diǎn)云地距亂序分布。此外，MLS結(jié)果滿(mǎn)足單調(diào)遞增條件，但在駐點(diǎn)鄰近區(qū)域誤差較大。如圖8(b)所示，駐點(diǎn)定位實(shí)現(xiàn)正確分區(qū)，斜距修正后的MLS結(jié)果不再滿(mǎn)足單調(diào)遞增條件，需要進(jìn)行迭代修正。圖8(c)是斜距修正MLS結(jié)果的地距排序結(jié)果，這里有很多非單調(diào)遞增點(diǎn)。圖8(d)示意的是地距-高程平面的重建結(jié)果。通過(guò)迭代執(zhí)行單調(diào)遞增修正和斜距修正等處理，重建結(jié)果滿(mǎn)足單調(diào)遞增約束條件。最后，執(zhí)行大地坐標(biāo)系投影變換和格網(wǎng)化處理，即可得到最終的重建結(jié)果。

圖8 本文方法對(duì)于切片5的主要中間結(jié)果Fig.8 The main intermediate results of the proposed method for slice 5

4.3 實(shí)測(cè)實(shí)驗(yàn)

4.3.1 數(shù)據(jù)介紹

2019年，我們團(tuán)隊(duì)在中國(guó)四川省樂(lè)山市市域內(nèi)的黃泥巴山區(qū)開(kāi)展了機(jī)載陣列TomoSAR飛行實(shí)驗(yàn)。系統(tǒng)工作在X-波段。天線(xiàn)陣列由沿水平向均勻的12個(gè)天線(xiàn)構(gòu)成，基線(xiàn)間隔是0.2 m。詳細(xì)的系統(tǒng)參數(shù)見(jiàn)表3。在獲得回波數(shù)據(jù)后，依次執(zhí)行SAR圖像生成、圖像配準(zhǔn)、高程分辨和高程模糊求解[10]。一個(gè)高程周期的離散點(diǎn)數(shù)是128。

表3 主要實(shí)驗(yàn)參數(shù)Tab.3 Main experimental parameters

4.3.2 結(jié)果與分析

圖9是觀(guān)測(cè)場(chǎng)景的SAR圖像。場(chǎng)景的像素?cái)?shù)是7500×6800 (方位×地距)。方位向和斜距向像素大小約為0.11 m和0.14 m。圖10呈現(xiàn)了SAR圖像中3個(gè)典型切片的地距-高度平面和斜距-高程平面的處理結(jié)果，其中藍(lán)色點(diǎn)為T(mén)omoSAR點(diǎn)云數(shù)據(jù)。除了切片2的 Trdis1=80,T r=40，所有切片投影變換參數(shù)同仿真數(shù)據(jù)1。抖動(dòng)因子均設(shè)為1.2。AW3D30 DSM和1:10000 DEM兩種外源數(shù)據(jù)用作參考值。DSM格網(wǎng)間隔為30 m，標(biāo)稱(chēng)測(cè)高精度為5 m。DEM格網(wǎng)間隔為5 m，平地高程中誤差高達(dá)1 m，但高山地高程中誤差是6.7 m。為了對(duì)比更多的細(xì)節(jié)，本文呈現(xiàn)的兩種外源數(shù)據(jù)是經(jīng)global mapper軟件1 m格網(wǎng)化處理后的結(jié)果。相比于DEM，DSM包含了地表地物的高度。在無(wú)誤差的情況下，對(duì)于同一觀(guān)測(cè)目標(biāo)，DSM數(shù)據(jù)的高度應(yīng)該不小于DEM。然而，在局部區(qū)域，DSM數(shù)據(jù)的高度小于DEM，如矩形框2。這種情況可能由插值處理或者是兩種數(shù)據(jù)本身的測(cè)高誤差導(dǎo)致。因此，本文只利用兩種外源數(shù)據(jù)來(lái)定性分析3種方法的重建性能。此外，實(shí)測(cè)數(shù)據(jù)存在斜距誤差，執(zhí)行幾何約束MLS前需要對(duì)同一高程向上的點(diǎn)取斜距平均進(jìn)行抑制。

圖9 觀(guān)測(cè)場(chǎng)景的SAR圖像(顏色圖表示歸一化后的散射強(qiáng)度)Fig.9 SAR image of the observation scene (the color map indicates the normalized scattering intensity)

首先，分析TomoSAR點(diǎn)云精度。整體來(lái)看，原始點(diǎn)云走勢(shì)與兩種外源數(shù)據(jù)幾乎一致。而且對(duì)于DSM高于DEM的區(qū)域(如矩形框1)，同樣包含地物高度信息的原始點(diǎn)云更接近于DSM，符合實(shí)際情況。需要說(shuō)明的是，本文呈現(xiàn)的原始點(diǎn)云數(shù)據(jù)存在非均勻基線(xiàn)導(dǎo)致的局部偏差問(wèn)題，如圖10(b)中高程值為128的點(diǎn)。但這種局部偏差仍然可以等效為高程誤差，并不影響本文方法有效性的驗(yàn)證。

然后，分析3種方法的重建性能?？梢钥闯觯琇S結(jié)果顯著偏離原始點(diǎn)云和兩種外源數(shù)據(jù)，MLS和所提方法結(jié)果走勢(shì)和兩種外源數(shù)據(jù)幾乎一致。但是，如圖10(c)兩個(gè)局部放大矩形框所示，MLS和所提方法結(jié)果高度差高達(dá)10 m。進(jìn)一步地，根據(jù)圖10(b)、圖10(d)和圖10(f)中矩形框可以推斷，所提方法處理結(jié)果更貼近原始點(diǎn)云，符合疊掩地形的S形曲線(xiàn)特征，而MLS結(jié)果在追求平滑的同時(shí)帶來(lái)了斜距誤差。因此，對(duì)于山區(qū)疊掩地形，尤其是駐點(diǎn)鄰近區(qū)域，所提方法重建精度更高。

圖10 實(shí)測(cè)實(shí)驗(yàn)的處理結(jié)果Fig.10 Results of the real experiment

最后，圖11呈現(xiàn)了觀(guān)測(cè)場(chǎng)景的光學(xué)圖像。矩形框示意了重建區(qū)域，其大小是800 m×1400 m (方位×地距)。圖12呈現(xiàn)了本文方法重建的三維圖像。該場(chǎng)景的SAR圖像見(jiàn)圖9。可以看出，TomoSAR山區(qū)三維重建圖像與光學(xué)圖像基本一致：兩座主峰分別位于重建的北側(cè)和東側(cè)，兩峰之間有明顯的溝壑；主峰1處的中部有一處凹陷。此外，本文結(jié)果和DSM的高度變化范圍都是600～1200 m。

圖11 觀(guān)測(cè)場(chǎng)景的光學(xué)圖像Fig.11 Optical image of the observation scene

圖12 三維重建圖像(顏色圖表示歸一化后的散射強(qiáng)度)Fig.12 The three-dimensional reconstruction image (the color map indicates the normalized scattering intensity)

4.4 影響因素分析

首先，分析高程向噪聲對(duì)駐點(diǎn)定位和地距投影的影響。對(duì)于實(shí)測(cè)點(diǎn)云數(shù)據(jù)中可能存在的離群點(diǎn)和虛假目標(biāo)，需要在高程模糊求解步驟進(jìn)行去除[10]，本文不做分析。噪聲對(duì)駐點(diǎn)定位的影響可以分為定位誤差、虛警和漏檢3種。定位誤差和虛警對(duì)地距投影以及重建精度影響不大。然而，漏檢會(huì)導(dǎo)致本文方法的重建精度急劇下降，且可能會(huì)顯著低于MLS。這是因?yàn)?，本文方法的核心思想是利用TomoSAR點(diǎn)云地距和高程間的單調(diào)遞增關(guān)系進(jìn)行高程修正，而駐點(diǎn)定位和地距投影是正確預(yù)判點(diǎn)云地距相對(duì)順序的關(guān)鍵。

地形遮擋是山區(qū)三維重建中較常見(jiàn)的現(xiàn)象。如圖8(a)和圖8(b)所示，I區(qū)位于地形上部，II區(qū)位于地形中部，III區(qū)位于地形底部。那么，駐點(diǎn)A相鄰區(qū)域幾乎不會(huì)被遮擋住。雖然駐點(diǎn)B相鄰區(qū)域以及II區(qū)底部可能因?yàn)榍皞?cè)地形遮擋或者照射不完全等問(wèn)題導(dǎo)致對(duì)應(yīng)點(diǎn)云數(shù)據(jù)丟失，如圖6中數(shù)據(jù)2和圖10中切片3，但本文方法對(duì)這種情況很魯棒。同時(shí)，根據(jù)TomoSAR點(diǎn)云的高程單調(diào)遞增關(guān)系可知，雖然在I區(qū)相距駐點(diǎn)較遠(yuǎn)的位置可能出現(xiàn)點(diǎn)云空洞(如圖10(f))，但不會(huì)破壞單調(diào)遞增特性。因此，本文方法對(duì)地形遮擋現(xiàn)象很魯棒，且會(huì)通過(guò)格網(wǎng)化插值處理對(duì)這些空洞區(qū)進(jìn)行填充。

綜上，本文方法對(duì)地形遮擋現(xiàn)象較魯棒，但對(duì)駐點(diǎn)漏檢導(dǎo)致的地距預(yù)判順序嚴(yán)重錯(cuò)誤的情況非常敏感。因此，對(duì)于被噪聲嚴(yán)重破壞的原始點(diǎn)云，要通過(guò)合適的投影參數(shù)設(shè)置、去噪等預(yù)處理甚至是人工干預(yù)來(lái)避免駐點(diǎn)漏檢等情況的發(fā)生。

5 結(jié)論

TomoSAR山區(qū)三維重建是三維成像領(lǐng)域的前沿研究。本文針對(duì)TomoSAR山區(qū)點(diǎn)云高程向定位精度差的問(wèn)題，推導(dǎo)了TomoSAR點(diǎn)云幾何約束模型，進(jìn)而提出了一種基于幾何約束MLS的三維重建方法。通過(guò)駐點(diǎn)定位和地距投影，雷達(dá)坐標(biāo)系表征的點(diǎn)可以被投影到新的方位-地距-高程坐標(biāo)系，幾何約束得以利用。通過(guò)在MLS代價(jià)函數(shù)中引入單調(diào)遞增約束懲罰項(xiàng)，所提方法可以針對(duì)性地對(duì)斜距修正后不滿(mǎn)足單調(diào)遞增約束的MLS擬合結(jié)果進(jìn)行高程向定位誤差修正，進(jìn)而提高了三維重建精度。仿真和實(shí)測(cè)的機(jī)載陣列TomoSAR山區(qū)數(shù)據(jù)處理結(jié)果驗(yàn)證了本文方法的有效性，對(duì)于仿真數(shù)據(jù)結(jié)果，本文方法相比LS和MLS重建精度提高一倍以上。而且，相比于DSM數(shù)據(jù)和DEM數(shù)據(jù)，本文重建結(jié)果在水平分辨率和散射信息恢復(fù)等方面具有顯著的優(yōu)勢(shì)，有望用于一些需要精細(xì)化地形地貌支撐的軍事偵察和災(zāi)難評(píng)估等應(yīng)用。未來(lái)，我們也將關(guān)注原始點(diǎn)云高精度生成等處理，以得到更高精度的TomoSAR山區(qū)三維重建結(jié)果。