韋立登 李永杰 孫中昶 高建

摘要干涉合成孔徑雷達(dá)（Interferometric Synthetic Aperture Radar，InSAR）已成為獲取高精度數(shù)字正射影像圖（Digital Orthophoto Map，DOM）和數(shù)字表面模型（Digital Surface Model，DSM）的關(guān)鍵技術(shù)之一，其不受天氣狀況的影響，可以全天時(shí)、全天候進(jìn)行數(shù)據(jù)獲取.機(jī)載雙天線毫米波InSAR不受失相關(guān)的影響，具有小體積、高分辨率、機(jī)動(dòng)靈活等特點(diǎn)，可以實(shí)現(xiàn)大尺度、高精度成像.本研究利用機(jī)載雙天線毫米波InSAR通過干涉批處理、區(qū)域網(wǎng)平差、地理編碼、圖像拼接鑲嵌等步驟生成了貴州施秉試驗(yàn)區(qū)（丘陵、山地）和四川邛崍?jiān)囼?yàn)區(qū)（高山地）高精度的DOM與DSM，并利用地面控制點(diǎn)（GCPs，Ground Control Points）進(jìn)行了精度驗(yàn)證，驗(yàn)證結(jié)果表明獲取的DSM精度符合1∶5 000地形圖繪制要求，表明了機(jī)載雙天線毫米波InSAR具備生成不同地形的DOM/DSM的能力，為解決困難地區(qū)DOM/DSM數(shù)據(jù)缺失問題提供了新的技術(shù)手段.關(guān)鍵詞機(jī)載毫米波InSAR;區(qū)域網(wǎng)平差;高精度;數(shù)字表面模型;數(shù)字正射影像圖

中圖分類號(hào)P225.7

文獻(xiàn)標(biāo)志碼A

0引言

InSAR可以全天時(shí)、全天候快速獲取大面積高精度的地表三維信息，已經(jīng)成為獲取高精度數(shù)字表面模型（Digital Surface Model，DSM）的重要技術(shù)手段[1].機(jī)載雙天線毫米波InSAR不受失相關(guān)因素影響，具有體積小、分辨率高、實(shí)時(shí)獲取、機(jī)動(dòng)靈活等特點(diǎn)，并且其波長短、穿透性弱，更加有助于DSM的獲取.

獲取DSM的手段主要包括外業(yè)測(cè)繪、衛(wèi)星遙感、航空攝影測(cè)量、激光雷達(dá)、干涉合成孔徑雷達(dá)等技術(shù).傳統(tǒng)的外業(yè)測(cè)繪耗時(shí)耗力，并不合適獲取大范圍的DSM;激光雷達(dá)雖然可以獲取高精度的DSM，但是其技術(shù)門檻較高;通過衛(wèi)星遙感、攝影測(cè)量等手段獲取我國西南地區(qū)DSM時(shí)，常常受云、霧、雨等天氣狀況的影響，獲取數(shù)據(jù)的質(zhì)量和精度會(huì)受到很大程度的影響，造成部分地區(qū)DSM數(shù)據(jù)缺失;而InSAR可以克服這些天氣狀況的影響，可以全天時(shí)、全天候、大范圍進(jìn)行數(shù)據(jù)獲取，在對(duì)地觀測(cè)領(lǐng)域起著舉足輕重的作用[2].

星載InSAR雖然覆蓋區(qū)域較廣，但是其衛(wèi)星軌道相對(duì)固定，不能隨意改變研究區(qū)域，而機(jī)載雙天線毫米波InSAR可以機(jī)動(dòng)靈活地選取研究區(qū)域，不受失相關(guān)因素的影響，并且毫米波可用頻帶寬，應(yīng)用較小口徑天線即可獲取更窄的天線波束以及更高的天線增益，可以實(shí)現(xiàn)距離向和方位向的高分辨率成像.另外，機(jī)載毫米波InSAR系統(tǒng)使用毫米波，波長較短，相較于其他InSAR系統(tǒng)，機(jī)載毫米波系統(tǒng)的天線和微波器件也會(huì)相應(yīng)減小，更加有利于實(shí)現(xiàn)雷達(dá)系統(tǒng)小型化，可以實(shí)現(xiàn)在無人機(jī)飛行平臺(tái)上應(yīng)用.機(jī)載毫米波InSAR使用Ka波段的電磁波，相比C、S、X、Ku等波段其波長短，圖像細(xì)節(jié)更加清晰，同時(shí)其波長更加接近可見光，其雷達(dá)圖像的效果要優(yōu)于其他InSAR圖像[3].

在國外一些發(fā)達(dá)國家，機(jī)載InSAR地形測(cè)繪技術(shù)已經(jīng)很成熟，并且可以提供業(yè)務(wù)化運(yùn)行的服務(wù)和產(chǎn)品.1951年，美國Goodyear宇航中心Carl Wiley最先提出InSAR技術(shù)，并在1969年首次將該技術(shù)應(yīng)用到金星表面的測(cè)量[4];1986年，Zebker和Goldstein采用機(jī)載InSAR數(shù)據(jù)得到了高程精度為2～10 m的舊金山大橋地區(qū)的三維地形圖，證明了機(jī)載InSAR獲取高精度地形信息的能力[5];20世紀(jì)90年代初，采用STAR-3i系統(tǒng)通過Lear 36飛機(jī)對(duì)巴拿馬運(yùn)河及周圍區(qū)域?qū)嵤y(cè)圖任務(wù)，獲取的DSM精度為3 m，格網(wǎng)間距為10 m，雷達(dá)影像分辨率為2.5 m[6];1998年，德國機(jī)

載InSAR系統(tǒng)AeS-1對(duì)巴西北部區(qū)域和委內(nèi)瑞拉南部進(jìn)行區(qū)域測(cè)圖，獲取原始雷達(dá)影像分辨率為0.5 m，DSM高程精度達(dá)到0.5 m[7];2000年，美國“奮進(jìn)號(hào)”航天飛機(jī)在佛羅里達(dá)州發(fā)射升空，歷時(shí)13 d獲取了北緯60°到南緯58°全球性的DSM，并且可以免費(fèi)獲取[8];2003年，Intermap公司進(jìn)行了NextMap美國計(jì)劃，獲得的數(shù)據(jù)產(chǎn)品包括精度為1 m，格網(wǎng)間距為5 m的DSM，雷達(dá)正射影像分辨率為1.25 m[9].我國機(jī)載SAR研制和星載SAR同時(shí)起步，目前已經(jīng)可以生產(chǎn)多種型號(hào)的機(jī)載SAR系統(tǒng)，同時(shí)也進(jìn)行了大量的試驗(yàn).2004年，中國科學(xué)院電子學(xué)研究所成功研制了X波段機(jī)載雙天線InSAR系統(tǒng)，并進(jìn)行試驗(yàn)，在地勢(shì)平坦地區(qū)取得了較好的效果[10];2011年，在國家“863”計(jì)劃的支持下，中國科學(xué)院電子學(xué)研究所研制出我國第一個(gè)機(jī)載毫米波三基線InSAR原理樣機(jī)，并于2011年5月在運(yùn)12飛機(jī)上進(jìn)行了試驗(yàn)[11];2013年，中國航天科工集團(tuán)第二研究院二十三所也研制出了機(jī)載毫米波SAR原理樣機(jī)，并進(jìn)行了大量飛行試驗(yàn)[12].

本研究利用機(jī)載雙天線毫米波InSAR系統(tǒng)通過干涉處理、區(qū)域網(wǎng)平差、地理編碼、圖像拼接鑲嵌等步驟，生成了貴州省施秉縣和四川省邛崍市2個(gè)試驗(yàn)區(qū)的DOM（Digital Orthophoto Map，數(shù)字正射影像圖）和DSM，試驗(yàn)區(qū)地形涉及丘陵、山地、高山區(qū)等，對(duì)解決丘陵、山區(qū)等困難地區(qū)DSM繪制提供參考.通過控制點(diǎn)驗(yàn)證，生成的DSM滿足1∶5 000地形圖繪制要求，表明了機(jī)載毫米波InSAR具備生成不同地形條件下大比例尺DOM/DSM的能力.

本研究提出：1）基于GPU并行處理技術(shù)，進(jìn)行干涉批處理，同時(shí)利用最小平衡樹的方法快速、精確地進(jìn)行相位展開;2）基于SIFT算法全自動(dòng)生成實(shí)驗(yàn)區(qū)內(nèi)影像間的連接點(diǎn)，極大提高了效率和精度;3）通過敏感度分析，基于視向量正交分解算法構(gòu)建三維重建模型，利用區(qū)域網(wǎng)平差方法對(duì)基線長度、基線傾角、相位偏移量3個(gè)參數(shù)進(jìn)行了精確的定標(biāo).

1試驗(yàn)區(qū)概況與載機(jī)參數(shù)

1.1試驗(yàn)區(qū)概況

貴州施秉試驗(yàn)區(qū)（圖1a）坐落于貴州省東南部，位于黔東南州苗族侗族自治州西北部施秉縣，經(jīng)緯度范圍為：108.184°108.185°E，26.889°～26.919°N，試驗(yàn)區(qū)長14.99 km，寬3.35 km，海拔范圍為560～1198 m.試驗(yàn)區(qū)地貌主要以高山和丘陵為主，同時(shí)也包含一些城鎮(zhèn)、道路、農(nóng)田、樹林等地物，整體地勢(shì)比較崎嶇.四川邛崍?jiān)囼?yàn)區(qū)（圖1b）坐落于四川省中部，位于邛崍市西部山區(qū)地帶，經(jīng)緯度范圍為：103.140°～103.283°E，30.319°～30.493°N，試驗(yàn)區(qū)由西南至東北呈帶狀分布，長21.27 km，寬2.56 km，海拔范圍為：603～1 344 m.實(shí)驗(yàn)區(qū)內(nèi)主要以高山為主，同時(shí)也包含一些城鎮(zhèn)和農(nóng)田，整體地勢(shì)非常崎嶇，起伏非常明顯.

1.2試驗(yàn)區(qū)數(shù)據(jù)

使用機(jī)載毫米波InSAR系統(tǒng)由西向東方向獲取了貴州施秉試驗(yàn)區(qū)3個(gè)條帶共21景機(jī)載數(shù)據(jù)，并且在試驗(yàn)區(qū)首末端布置了26個(gè)GCPs（地面控制點(diǎn)）用于區(qū)域網(wǎng)平差與結(jié)果驗(yàn)證，其中14個(gè)用于區(qū)域網(wǎng)平差，12個(gè)用于DSM和DOM結(jié)果驗(yàn)證.

使用機(jī)載毫米波InSAR系統(tǒng)由西南至東北方向獲取了四川邛崍?jiān)囼?yàn)區(qū)1個(gè)條帶共13景機(jī)載數(shù)據(jù)，并且在試驗(yàn)區(qū)首末端布置了25個(gè)GCPs，其中12個(gè)用于區(qū)域網(wǎng)平差，13個(gè)用于DSM和DOM結(jié)果驗(yàn)證.

1.3載機(jī)參數(shù)

采用機(jī)載雙天線毫米波InSAR系統(tǒng)獲取了2個(gè)試驗(yàn)區(qū)的機(jī)載數(shù)據(jù)，載機(jī)參數(shù)如表1所示.

載機(jī)搭載高精度的航空定位測(cè)姿系統(tǒng)（POS AV610），完美地繼承了全球?qū)Ш叫l(wèi)星系統(tǒng)和慣性導(dǎo)航系統(tǒng)，POS AV610在每秒鐘內(nèi)可以對(duì)航空傳感器進(jìn)行上百次的精確定位定向，能夠精確地獲取載機(jī)的三維坐標(biāo)、姿態(tài)角、速度、加速度等實(shí)時(shí)載機(jī)信息，同時(shí)使用戶在處理數(shù)據(jù)時(shí)無需地面參考信息，消除了繁瑣耗時(shí)的航空攝影空三解算帶來的麻煩，極大提高了數(shù)據(jù)獲取的精度與效率.

2地形三維重建原理

本研究基于視向量正交分解[13-14]進(jìn)行了地形三維重建，將視向量從載機(jī)移動(dòng)坐標(biāo)系轉(zhuǎn)換到地固坐標(biāo)系，三維重建模型如圖2所示.

圖2a為Madsen引入建立的移動(dòng)坐標(biāo)系（Madsen Moving Coordinates，MMC），V軸為平臺(tái)速度矢量方向，W軸為速度矢量與基線矢量的叉積，N軸由V軸和W軸由右手法則共同決定，，，分別為三軸的單位矢量.b表示基線矢量，bv與bn為b在順軌和交軌方向的分量，b與水平面的夾角為θb，bn與水平面的夾角為α.

移動(dòng)坐標(biāo)系中，單位視向量可由、、表示：

=μ+η+ξ. （1）

設(shè)MMC坐標(biāo)系下單位視向量表示為

=μηξ=sin βsin θ1-bvbnsin β-k2cos2 β-sin θ1-bvbnsin β2， （2）

其中，β為斜視角，表示為

sin β=λf/2v. （3）

θ1定義為回波到達(dá)方向角：

θ1=arcsinb22rbn+k1k2λφ2Qπbn-λ2φ28Q2π2rbn， （4）

Q為雷達(dá)工作模式，Q=1時(shí)為標(biāo)準(zhǔn)模式，Q=2時(shí)為乒乓模式;k1=1時(shí)表示右側(cè)視方式，k1=-1時(shí)表示左側(cè)視方式;k2=-1時(shí)表示主天線位于左側(cè)，k2=1時(shí)表示主天線位于右側(cè).

圖2b為平臺(tái)坐標(biāo)系，原點(diǎn)為載機(jī)POS相位中心，X軸指向載機(jī)機(jī)頭方向，Z軸為過原點(diǎn)垂直向下方向，Y軸垂直于X軸、Z軸構(gòu)成的平面，指向飛行方向右側(cè).載機(jī)在飛行過程中的姿態(tài)通過橫滾角θroll、俯仰角θpitch、偏航角θyaw來表示.旋轉(zhuǎn)矩陣表示為

Rroll=1000cos θroll-sin θroll0sin θrollcos θroll， （5）

Rpitch=cos θpitch0-sin θpitch010sin θpitch0cos θpitch， （6）

Ryaw=cos θyaw-sin θyaw0sin θyawcos θyaw0001. （7）

圖2c為機(jī)載InSAR三維重建模型，A1和A2分別表示載機(jī)主副天線相位中心的位置，H為主天線相位中心的高度.b表示基線矢量，由主天線指向副天線.P為地物點(diǎn)，高程為h，主天線相位中心到地物點(diǎn)P的矢量為r1，矢量長度為r，此方向上的單位矢量設(shè)為.O-XYZ為航跡坐標(biāo)系，X軸指向理想航機(jī)方向，Z軸指向正上方向，Y軸垂直于X、Z軸構(gòu)成的平面，并且指向左側(cè)方向.A和P分別為O到主天線相位中心和地物點(diǎn)的矢量.由移動(dòng)坐標(biāo)系到O-XYZ坐標(biāo)系的旋轉(zhuǎn)矩陣為

Γ=1000k2cos θbnk1k2sin θbn0-k1k2sin θbnk2cos θbn=

1000k21-b2b2n sin2θbk1k2bbnsin θb0-k1k2bbnsin θbk21-b2b2n sin2θb. （8）

由式（1）—（8）可得，矢量P表示為

P=A+r=A+rRyawRpitchRrollΓvnw， （9）

式（9）為理想航跡坐標(biāo)系下地物點(diǎn)位置的三維表達(dá).

在正側(cè)視的條件下，雷達(dá)波束中心面即為零多普勒面，此時(shí)β=0.假設(shè)主天線位于載機(jī)左側(cè)，并且右側(cè)視情況下，地形高程表示為

h=H-r（cos θpcos（α-θr）cos θ1+cos θpsin（α-θr）sin θ1）， （10）

其中，θ1為回波角，定義為

θ1=sin-1b2r-λ2πb-λ228π2rb. （11）

3數(shù)據(jù)處理流程

根據(jù)機(jī)載InSAR獲取的SLC數(shù)據(jù)，采用自主開發(fā)的機(jī)載InSAR數(shù)據(jù)處理軟件AirborneInSARMap，使用高性能的GPU并行技術(shù)，通過粗配準(zhǔn)、精配準(zhǔn)、干涉處理、相位濾波、相位解纏等干涉批處理步驟，批量、快速地獲取測(cè)區(qū)內(nèi)每景影像對(duì)應(yīng)的解纏相位數(shù)據(jù).基于尺度不變特征變換（Scale-Invariant Feature Transform，SIFT）算法，自動(dòng)選取測(cè)區(qū)內(nèi)影像間的連接點(diǎn)，建立影像之間的約束關(guān)系，使用最小二乘方法，通過區(qū)域網(wǎng)平差干涉定標(biāo)方法來校正所有影像的干涉參數(shù)（基線長度、基線傾角、相位偏移量）、求解連接點(diǎn)的三維坐標(biāo)、求解影像的三維坐標(biāo).由設(shè)置的DOM/DSM像元大小，通過cubic convolution三次卷積插值方法對(duì)求得的平面坐標(biāo)以及高程值進(jìn)行插值處理，生成每景影像對(duì)應(yīng)的DOM/DSM.通過拼接鑲嵌對(duì)所有影像進(jìn)行拼接處理后，根據(jù)國家標(biāo)準(zhǔn)比例尺進(jìn)行影像輸出.

利用機(jī)載InSAR生成DSM數(shù)據(jù)處理流程如圖3所示.

3.1干涉處理

干涉處理是生成DSM過程中的第一步，主要包括復(fù)圖像配準(zhǔn)、預(yù)濾波、干涉圖生成、去平地效應(yīng)、干涉圖濾波、相位解纏等步驟，本研究在處理過程中省略了預(yù)濾波與去平地效應(yīng)這兩部分.

3.1.1復(fù)圖像配準(zhǔn)與干涉圖生成

在生成DSM的整個(gè)過程中，主、輔圖像配準(zhǔn)是最基礎(chǔ)、最關(guān)鍵的一步，其核心思想在于確定兩幅圖像之間匹配位置的相對(duì)偏移量.SAR影像配準(zhǔn)不精確，將導(dǎo)致較大的干涉相位誤差.配準(zhǔn)精度達(dá)到1/10像素才對(duì)干涉條紋圖的質(zhì)量沒有明顯影響[15].常用匹配測(cè)度主要有3種：基于復(fù)數(shù)值的相似性測(cè)度、基于相位值的相似性測(cè)度和基于強(qiáng)度值的相似性測(cè)度[16].本研究利用基于FFT的復(fù)相關(guān)精確配準(zhǔn)獲取同名點(diǎn)的精確偏移量，對(duì)匹配數(shù)據(jù)進(jìn)行多項(xiàng)式擬合后計(jì)算主輔圖像對(duì)的坐標(biāo)轉(zhuǎn)換關(guān)系，最后對(duì)輔圖像進(jìn)行重采樣完成配準(zhǔn)[17].配準(zhǔn)后將主圖像的復(fù)數(shù)值與輔圖像的復(fù)數(shù)值進(jìn)行共軛相乘后得到干涉圖.

3.1.2干涉圖濾波

由于SAR固有的斑點(diǎn)噪聲、雷達(dá)系統(tǒng)熱噪聲、雷達(dá)陰影、配準(zhǔn)誤差等存在使得干涉圖存在很多相位噪聲，表現(xiàn)為相位不連續(xù)、周期性不明顯、相位條紋圖不明顯，給后續(xù)解纏造成很大影響[18-19].常用的濾波方法包括空間域?yàn)V波、頻率域?yàn)V波、時(shí)頻分析濾波，本研究利用Goldstein濾波[20]對(duì)干涉相位進(jìn)行濾波，選取相互重疊的相位塊在頻率域采用平滑濾波器處理后，采用濾波參數(shù)對(duì)其功率譜進(jìn)行處理.

3.1.3相位解纏

經(jīng)干涉處理得到的干涉圖相位為相位主值，其值域在[0，2π]或[-π，π].為了獲取精確的地形高程信息，需要得到干涉圖對(duì)應(yīng)的絕對(duì)干涉相位值.相位解纏采用一定的數(shù)學(xué)方法或者計(jì)算方法對(duì)干涉圖進(jìn)行處理，得到各干涉相位之間相差的整周期數(shù)，從而獲取連續(xù)變化的干涉相位的過程[21-22].本研究提出了一種基于最小平衡樹（Minimum Balanced Trees，MBTs）的易于并行實(shí)現(xiàn)的快速、精確的解纏方法，利用高性能的GPU并行技術(shù)，構(gòu)建最小平衡樹/林和推算解纏優(yōu)先順序圖進(jìn)行相位解纏[22-23].

3.2區(qū)域網(wǎng)平差

獲取DSM過程中，從POS中得到的系統(tǒng)參數(shù)會(huì)存在誤差，這些參數(shù)誤差會(huì)通過相高轉(zhuǎn)換、地理編碼等步驟傳遞到DSM中，影響DSM的精度，所以需要對(duì)這些參數(shù)進(jìn)行定標(biāo)處理.傳統(tǒng)的干涉參數(shù)定標(biāo)方法基于GCPs對(duì)系統(tǒng)參數(shù)進(jìn)行定標(biāo)處理.假設(shè)對(duì)基線長度、基線傾角、相位偏移量進(jìn)行定標(biāo)，每景影像至少需要3個(gè)GCPs，適合小范圍并且易于布設(shè)控制點(diǎn)的區(qū)域.對(duì)于大范圍、多條帶或者不易布設(shè)控制點(diǎn)的區(qū)域（高山、丘陵、沼澤等），傳統(tǒng)的干涉參數(shù)定標(biāo)的方法就不適用了.

因此本研究基于視向量正交分解算法構(gòu)建三維重建模型，利用少量控制點(diǎn)和影像間的連接點(diǎn)，建立影像之間的約束關(guān)系，進(jìn)行區(qū)域網(wǎng)平差處理[24-26]，校正干涉參數(shù)與連接點(diǎn)坐標(biāo).

3.2.1連接點(diǎn)選取

基于SIFT算法[27-28]自動(dòng)、快速選取測(cè)區(qū)內(nèi)影像間的連接點(diǎn).同一條帶內(nèi)，滿足相干系數(shù)大于0.9，誤差小于0.5個(gè)像素的要求;條帶間，滿足誤差小于3個(gè)像素，同時(shí)可以通過人工復(fù)檢對(duì)誤差大于1個(gè)像素的連接點(diǎn)進(jìn)行剔除.

3.2.2敏感度分析

由式（10）可得，影響DSM高程精度的系統(tǒng)參數(shù)有主天線相位中心高度H、斜距r、基線長度b、基線傾角α、絕對(duì)干涉相位φ、橫滾角θroll、俯仰角θpitch等.設(shè)這些參數(shù)的誤差分別為dh，dH，dr，db，dα，dφ，dθroll，dθpitch，假設(shè)系統(tǒng)參數(shù)之間相互不影響，由協(xié)方差傳播定律可以得：

σ2h=dhdH2σ2H+dhdr2σ2r+dhdb2σ2b+dhdα2σ2α+dhdφ2σ2φ+dhdθroll2σ2θroll+dhdθpitch2σ2pitch， （12）

其中：

dhdH=1， （13）

dhdr=-cosθ+lλ2φ28π2rb2-b2r， （14）

dhdb=rl12r+λφ2πb2+λ2φ28π2rb2， （15）

dhdα=rlcos2β-sin2θ1， （16）

dhdφ=rl-λ2πb-λ2φ4π2rb， （17）

dhdθroll=-rlcos2β-sin2θ1， （18）

dhdθpitch=r[sinθpitchcos（α-θroll）cos2β-sin2θ1-sinθpitchsin（α-θroll）sinθ1]， （19）

其中系數(shù)l為

l=cosθpitchcos（α-θroll）sinθ1+sin（α-θroll）cosθ1cosθ1. （20）

為了分析各個(gè)參數(shù)與地面點(diǎn)高程之間的敏感度，設(shè)橫滾角為0.023°，俯仰角為0.649°，偏航角為2.254°，視角范圍為20°～70°，平臺(tái)高度為4 043.116 m，地面點(diǎn)參考高程為900 m，敏感度分析結(jié)果如圖4所示.

如圖4所示，隨著視角的增加，基線長度、基線傾角、橫滾角敏感度值為103級(jí)，遠(yuǎn)大于其他參數(shù).由于載機(jī)搭載高精度的POS AV610，可以精確地計(jì)算載機(jī)的相位中心高度H、斜距r.同時(shí)俯仰角誤差對(duì)DSM高程影響很小，可以忽略不計(jì).橫滾角可以通過校正基線傾角以消除它對(duì)高程的影響.同時(shí)由于系統(tǒng)噪聲影響，解纏相位與真實(shí)相位之間存在常數(shù)差，即相位偏移量，必須對(duì)其進(jìn)行精確定標(biāo)處理.因此需要對(duì)每個(gè)干涉像對(duì)的基線長度、基線傾角、相位偏移量進(jìn)行區(qū)域網(wǎng)平差定標(biāo)處理.

3.2.3區(qū)域網(wǎng)平差

區(qū)域網(wǎng)平差流程如圖5所示.

平差處理過程中認(rèn)為地面控制點(diǎn)（角反射器）三維坐標(biāo)精確，無需進(jìn)行校正，對(duì)于地面控制點(diǎn)，設(shè)

FGCP（i，j）（b，α，φ）=[0;0;0]. （12）

對(duì)于連接點(diǎn)，其三維坐標(biāo)為平差前給定的初值，需要對(duì)其三維坐標(biāo)進(jìn)行校正，設(shè)

FTP（i，k）（b，α，φ，X，Y，Z）=[0;0;0]， （13）

其中，GCP（i，j）表示第i個(gè)干涉像對(duì)上第j個(gè)控制點(diǎn)，TP（i，k）表示第i個(gè)干涉像對(duì)上第k個(gè)連接點(diǎn).

因式（12）、（13）為非線性方程，根據(jù)泰勒公式對(duì)其線性化得：

V=AΔx1+BΔx2-L， （14）

式（14）中，A為干涉參數(shù)系數(shù)矩陣，B為控制點(diǎn)和連接點(diǎn)的系數(shù)矩陣，Δx1為系統(tǒng)干涉參數(shù)的改正值，Δx2為連接點(diǎn)三維坐標(biāo)的改正值，L為常數(shù)項(xiàng).

對(duì)于i個(gè)干涉像對(duì)上第j個(gè)控制點(diǎn)：

VGCP_x（i，j）VGCP_y（i，j）VGCP_z（i，j）=

FGCP_x（i，j）bFGCP_x（i，j）αFGCP_x（i，j）

FGCP_y（i，j）bFGCP_y（i，j）αFGCP_y（i，j）

FGCP_z（i，j）bFGCP_z（i，j）αFGCP_z（i，j）·

ΔbΔαΔφ-LGCP_xLGCP_yLGCP_z. （15）

對(duì)于i個(gè)干涉像對(duì)上第k個(gè)連接點(diǎn)：

VTP_x（i，k）VTP_y（i，k）VTP_z（i，k）=

FTP_x（i，k）bFTP_x（i，k）αFTP_x（i，k）

FTP_y（i，k）bFTP_y（i，k）αFTP_y（i，k）

FTP_z（i，k）bFTP_z（i，k）αFTP_z（i，k）

ΔbΔαΔ-LTP_x（i，k）LTP_y（i，k）LTP_z（i，k）-ΔXTP（i，k）ΔYTP（i，k）ΔZTP（i，k）. （16）

基于最小二乘原理，求解方程時(shí)先消去未知數(shù)較多的連接點(diǎn)三維坐標(biāo)的改正數(shù)，求解系統(tǒng)干涉參數(shù)的改正量，再消去系統(tǒng)干涉參數(shù)改正量，求解連接點(diǎn)的坐標(biāo)，采用迭代以及逐步趨近的方法求解模型參數(shù).

本研究利用少量控制點(diǎn)和大量連接點(diǎn)對(duì)貴州施秉試驗(yàn)區(qū)、四川邛崍?jiān)囼?yàn)區(qū)進(jìn)行區(qū)域網(wǎng)平差處理，求解得到了每個(gè)區(qū)域中每景影像對(duì)應(yīng)的系統(tǒng)干涉參數(shù)與連接點(diǎn)三維坐標(biāo).

3.3地理編碼與拼接鑲嵌

根據(jù)區(qū)域網(wǎng)平差與相高轉(zhuǎn)換得到的每景影像對(duì)應(yīng)的三維坐標(biāo)，根據(jù)DOM/DSM所需分辨率（本研究設(shè)置為0.5 m），去除影像周邊50像素的邊緣值，根據(jù)相干閾值（本研究設(shè)置為0.6），采取re_ts6p插值方式對(duì)DOM/DSM進(jìn)行插值計(jì)算，得到分辨率為0.5 m的DOM和格網(wǎng)間距為0.5 m的DSM.

對(duì)測(cè)區(qū)內(nèi)所有的DOM/DSM進(jìn)行拼接鑲嵌處理，對(duì)于拼接后的DOM/DSM存在拼接縫的情況要進(jìn)行去縫處理;對(duì)于拼接后DOM在拼接處兩側(cè)顏色不均勻的情況要根據(jù)直方圖匹配的方法進(jìn)行勻色.

4試驗(yàn)結(jié)果和分析

4.1試驗(yàn)結(jié)果

本研究利用自主研發(fā)的機(jī)載InSAR數(shù)據(jù)處理軟件AirborneInSARMap，基于干涉處理、區(qū)域網(wǎng)平差、地理編碼、拼接鑲嵌等流程快速、精確地生成了貴州施秉試驗(yàn)區(qū)和四川邛崍?jiān)囼?yàn)區(qū)0.5 m分辨率的DOM以及格網(wǎng)間距為0.5 m的DSM，分別如圖6、圖7所示.

4.2精度驗(yàn)證

為了驗(yàn)證生成DSM的精度，本研究分別利用12、13個(gè)地面控制點(diǎn)（角反射器）對(duì)生成的DSM的X，Y以及H 3個(gè)方向進(jìn)行驗(yàn)證，2個(gè)試驗(yàn)區(qū)驗(yàn)證結(jié)果如圖8所示.

由圖8可得，貴州施秉試驗(yàn)區(qū)3個(gè)方向的平均值分別為0.588、-0.534和0.340 m，RMSE分別為1.490、0.938和1.433 m;四川邛崍?jiān)囼?yàn)區(qū)3個(gè)方向的平均值分別為-0.049、0.075和-0.304 m，RMSE分別為1.410、2.156和1.846 m.2個(gè)試驗(yàn)區(qū)的平面定位誤差和高程中誤差均滿足1∶5 000地形圖制圖要求[29].

5結(jié)論

本研究利用自主研發(fā)的機(jī)載InSAR數(shù)據(jù)處理軟件AirborneInSARMap生成了貴州施秉試驗(yàn)區(qū)（丘陵、山地）和四川邛崍?jiān)囼?yàn)區(qū)（高山地）高精度的DOM與DSM，并利用控制點(diǎn)進(jìn)行了精度驗(yàn)證與分析，得出以下結(jié)論：

1）通過對(duì)2個(gè)實(shí)驗(yàn)區(qū)精度驗(yàn)證，高程中誤差可以滿足1∶5 000地形圖制圖精度要求，表明機(jī)載雙天線毫米波InSAR具備生成不同地形條件的DOM/DSM的能力，為困難區(qū)域DOM/DSM獲取提供了新的技術(shù)手段.

2）基于GPU并行處理技術(shù)，利用最小平衡樹解纏方法進(jìn)行相位展開，極大提高了解纏的效率.

3）基于視向量正交分解算法，利用SIFT算法自動(dòng)選取測(cè)區(qū)內(nèi)的連接點(diǎn)以及測(cè)區(qū)內(nèi)少量控制點(diǎn)，通過區(qū)域網(wǎng)平差對(duì)各干涉像對(duì)的基線長度、基線傾角和相位偏移量進(jìn)行了定標(biāo)，提高了DOM/DSM的精度.

本研究在數(shù)據(jù)處理過程中也存在一些不足：由于山區(qū)地帶存在大量陰影及疊掩，會(huì)導(dǎo)致解纏效果不理想，最終影響生成DSM的精度;對(duì)于由陰影與疊掩導(dǎo)致的黑洞問題，可以采用插值和對(duì)飛處理將DSM補(bǔ)全.

參考文獻(xiàn)

References

[1]Sun Z C，Guo H D，Li X W，et al.DEM generation and error analysis using the first Chinese airborne dual-antenna interferometric SAR data[J].International Journal of Remote Sensing，2011，32（23）：8485-8504

[2]孫中昶，郭華東，李新武.機(jī)載雙天線InSAR數(shù)據(jù)生成高精度DEM的誤差分析[J].高技術(shù)通訊，2012，22（2）：171-179

SUN Zhongchang，GUO Huadong，LI Xinwu.Error analysis of high-precision DEM generated from airborne dual-antenna interferometric SAR data[J].Chinese High Technology Letters，2012，22（2）：171-179

[3]張琦，張雷，郭俊棟，等.基于無人飛行平臺(tái)的小型Ka波段合成孔徑雷達(dá)系統(tǒng)研制[J].地球信息科學(xué)學(xué)報(bào)，2019，21（4）：524-531

ZHANG Qi，ZHANG Lei，GUO Jundong，et al.Development of Ka-band miniature synthetic aperture radar based on UAV[J].Journal of Geo-Information Science，2019，21（4）：524-531

[4]Rogers A E E，Ingalls R P.Venus：mapping the surface reflectivity by radar interferometry[J].Science，1969，165（3895）：797-799

[5]Zebker H A，Goldstein R M.Topographic mapping from interferometric synthetic aperture radar observations[J].Journal of Geophysical Research Atmospheres，1986，91（B5）：4993-4999

[6]牛瑞，袁軍.國外機(jī)載InSAR系統(tǒng)發(fā)展現(xiàn)狀及應(yīng)用分析[J].測(cè)繪科學(xué)與工程，2011（3）：70-74

NIU Rui，YUAN Jun.Development and application analysis of airborne InSAR system abroad[J].Geomatic Science and Engineering，2011（3）：70-74

[7]Nuri A L N.Mapping of large areas in tropical countries by using high resolution airborne interferometric radar[C]∥ISPRS，2000

[8]王聰.編隊(duì)干涉SAR對(duì)地測(cè)繪任務(wù)規(guī)劃方法研究[D].哈爾濱：哈爾濱工業(yè)大學(xué)，2014

WANG Cong.Research on the scheduling method of the formation flying InSAR mapping mission[D].Harbin：Harbin Institute of Technology，2014

[9]Mercer B.DEMs created from airborne IFSAR：an update[J].International Archives of Photogrammetry and Remote Sensing，2004，35（part B）

[10]Xiang M S，Wu Y R，Li S E，et al.Introduction on an experimental airborne InSAR system[C]∥Proceedings of IEEE International Geoscience and Remote Sensing Symposium，2005，7：4809-4812

[11]Li D J，Liu B，Pan Z H，et al.Airborne MMW InSAR interferometry with cross-track three-baseline antennas[C]∥9th European Conference on Synthetic Aperture Radar，2012：301-303

[12]師君，馬龍，韋順軍，等.基于導(dǎo)航數(shù)據(jù)的Ka波段InSAR成像處理與分析[J].雷達(dá)學(xué)報(bào)，2014，3（1）：19-27

SHI Jun，MA Long，WEI Shunjun，et al.Ka-band InSAR imaging and analysis based on IMU data[J].Journal of Radars，2014，3（1）：19-27

[13]Madsen S N，Zebker H A，Martin J.Topographic mapping using radar interferometry：processing techniques[J].IEEE Transactions on Geoscience and Remote Sensing，1993，31（1）：246-256

[14]Li Y W，Xiang M S，Lü X L，et al.Joint interferometric calibration based on block adjustment for an airborne dual-antenna InSAR system[J].International Journal of Remote Sensing，2014，35（17）：6444-6468

[15]張登榮，俞樂.一種高精度的干涉雷達(dá)復(fù)數(shù)影像配準(zhǔn)方法[J].遙感學(xué)報(bào)，2007，11（4）：563-567

ZHANG Dengrong，YU Le.A high-precision co-registration method for InSAR image processing[J].Journal of Remote Sensing，2007，11（4）：563-567

[16]馬凱，王仁禮，楊慶慶，等.機(jī)載InSAR圖像配準(zhǔn)算法比較研究[J].地理信息世界，2016，23（4）：50-53

MA Kai，WANG Renli，YANG Qingqing，et al.The comparison of co-registration algorithms of airborne InSAR image[J].Geomatics World，2016，23（4）：50-53

[17]孫中昶，郭華東，焦孟梅，等.機(jī)載雙天線InSAR復(fù)圖像自動(dòng)配準(zhǔn)研究[J].國土資源遙感，2010，22（1）：24-29

SUN Zhongchang，GUO Huadong，JIAO Mengmei，et al.The automatic registration of airborne dual-antenna interferometric SAR complex images[J].Remote Sensing for Land & Resources，2010，22（1）：24-29

[18]張斌，胡慶榮，韋立登，等.改進(jìn)的形態(tài)學(xué)干涉圖濾波方法[J].系統(tǒng)工程與電子技術(shù)，2018，40（10）：2230-2236

ZHANG Bin，HU Qingrong，WEI Lideng，et al.Improved morphological filtering algorithm of interferograms[J].Systems Engineering and Electronics，2018，40（10）：2230-2236

[19]Ambrosino R，Baselice F，F(xiàn)erraioli G，et al.Extended Kalman filter for multichannel InSAR height reconstruction[J].IEEE Transactions on Geoscience and Remote Sensing，2017，55（10）：5854-5863

[20]Goldstein R M，Werner C L.Radar interferogram filtering for geophysical applications[J].Geophysical Research Letters，1998，25（21）：4035-4038

[21]Goldstein R M，Zebker H A，Werner C L.Satellite radar interferometry：two-dimensional phase unwrapping[J].Radio Science，1988，23（4）：713-720

[22]Gao J，Sun Z C.Phase unwrapping method based on parallel local minimum reliability dual expanding for large-scale data[J].Journal of Applied Remote Sensing，2019，13（3）：038506

[23]Jian G.Reliability-map-guided phase unwrapping method[J].IEEE Geoscience and Remote Sensing Letters，2016，13（5）：716-720

[24]張薇，向茂生，吳一戎.基于三維重建模型的機(jī)載雙天線干涉SAR外定標(biāo)方法及實(shí)現(xiàn)[J].遙感技術(shù)與應(yīng)用，2009，24（1）：82-87

ZHANG Wei，XIANG Maosheng，WU Yirong.Realization of outside calibration method based on the sensitivity equation for dual-antenna airborne interferometric SAR[J].Remote Sensing Technology and Application，2009，24（1）：82-87

[25]Jin G W，Xiong X，Xu Q，et al.Baseline estimation algorithm with block adjustment for multi-pass dual-antenna insar[J].ISPRS-International Archives of the Photogrammetry，Remote Sensing and Spatial Information Sciences，2016，XLI-B7：39-45

[26]Yue X J，Han C M，Dou C Y，et al.Research on block adjustment of airborne InSAR images[J].IOP Conference Series：Earth and Environmental Science，2014，17：012199

[27]Lowe D G.Distinctive image features from scale-invariant keypoints[J].International Journal of Computer Vision，2004，60（2）：91-110

[28]Schwind P，Suri S，Reinartz P，et al.Applicability of the SIFT operator to geometric SAR image registration[J].International Journal of Remote Sensing，2010，31（8）：1959-1980

[29]孫曉云，王曉東，苗小利.INSAR干涉處理獲得DSM的產(chǎn)品質(zhì)量分析[J].礦山測(cè)量，2015（4）：72-75

SUN Xiaoyun，WANG Xiaodong，MIAO Xiaoli.Product quality analysis of DSM obtained by INSAR interference processing[J].Mine Surveying，2015（4）：72-75

On topographic surveying in tough areas using airborne millimeter-wave InSAR

WEI Lideng1LI Yongjie2，3SUN Zhongchang2，4GAO Jian5

1Beijing Institute of Radio Measurement，Beijing100854

2Key Laboratory of Digital Earth Science，Aerospace Information Research Institute，Chinese Academy of Sciences，Beijing100094

3School of Land Science and Technology，China University of Geosciences，Beijing100083

4Hainan Provincial Key Laboratory of Earth Observation，Sanya572029

5School of Geographic and Biologic Information，Nanjing University of Posts and Telecommunications，Nanjing210003

AbstractInterferometric synthetic aperture radar （InSAR） has become one of the key technologies to obtain high-precision Digital Orthophoto Map （DOM） and Digital Surface Model （DSM）.It is not affected by weather conditions and can acquire data in all day and weather conditions.Airborne dual-antenna millimeter-wave InSAR is independent of loss-of-correlation which has the characteristics of small size，high resolution，high flexibility，etc.And it can obtain large-scale and high-precision images.This paper uses airborne dual-antenna millimeter-wave InSAR to acquire high-precision DOM and DSM in Shibing experimental area of Guizhou （mountainous） and Qionglai experimental area of Sichuan （alpine） through interference processing，block adjustment，geocoding，and image mosaic.Besides，the ground control points （GCPs） are used to verify the accuracy.The results show that the accuracy of DSM obtained meets the requirement of 1：5000 terrain mapping.It has been proven that the airborne dual-antenna millimeter-wave InSAR has the ability to generate DOM/DSM with different terrains，which provides a new technical means for solving the lack of DOM/DSM data in tough areas.

Key wordsairborne millimeter-wave InSAR;block adjustment;high precision;digital surface model （DSM）;digital orthophoto map （DOM）

收稿日期2019-10-21

資助項(xiàng)目海南省重點(diǎn)研發(fā)計(jì)劃項(xiàng)目（ZDYF2019008）

作者簡介韋立登，男，博士，研究員，主要研究方向?yàn)楦缮婧铣煽讖嚼走_(dá)成像處理技術(shù)、運(yùn)動(dòng)補(bǔ)償技術(shù)、自動(dòng)配準(zhǔn)技術(shù)、干涉處理技術(shù)等.koridge@163.com

孫中昶（通信作者），男，博士，副研究員，主要研究方向?yàn)槌鞘羞b感、微波遙感等.sunzc@aircas.ac.cn

1北京無線電測(cè)量研究所，北京，100854

2中國科學(xué)院空天信息創(chuàng)新研究院/數(shù)字地球重點(diǎn)實(shí)驗(yàn)室，北京，100094

3中國地質(zhì)大學(xué)土地科學(xué)技術(shù)學(xué)院，北京，100083

4海南省地球觀測(cè)重點(diǎn)實(shí)驗(yàn)室，三亞，572029

5南京郵電大學(xué)地理與生物信息學(xué)院，南京，210003