張延冰，郭華東，韓春明

(1.中國科學(xué)院遙感與數(shù)字地球研究所，北京 100094;2.中國科學(xué)院大學(xué)，北京 100049)

0 引言

合成孔徑雷達干涉測量(interferometric synthetic aperture Radar，InSAR)技術(shù)具有全天時、全天候，測量精度和效率高等特點，已成為當(dāng)前快速獲取高精度地面三維信息的重要手段之一[1]。在熱帶雨林以及地形復(fù)雜的高山、極地地區(qū)，InSAR已成為獲取DEM的重要技術(shù)手段。雙天線InSAR系統(tǒng)具有時效性和自動化程度高等特點，且不存在時間失相關(guān)問題[2-3]，在地形制圖方面有著巨大的優(yōu)勢。自1974年Graham利用機載InSAR技術(shù)獲取了能滿足1∶25萬比例尺地形圖要求的高程數(shù)據(jù)以來，該技術(shù)越來越多地應(yīng)用到地形測量中，并被證實了具有獲取高精度DEM的能力[3-4]。

在科技發(fā)達國家，機載InSAR技術(shù)已被應(yīng)用于地形測繪方面，用以提供業(yè)務(wù)化運行的服務(wù)和產(chǎn)品。如InterMap公司的STAR-3i系統(tǒng)已實現(xiàn)商業(yè)運行多年，實現(xiàn)了巨大的經(jīng)濟效益。國內(nèi)機載InSAR系統(tǒng)研究相對較晚[5]，在測繪領(lǐng)域的應(yīng)用仍處于起步階段。2004年，中國科學(xué)院電子學(xué)研究所成功研制了“X波段機載雙天線InSAR系統(tǒng)”[6]，并進行了地形測繪飛行試驗，在地勢平坦地區(qū)取得了較好的效果。但在地形復(fù)雜區(qū)域，由于缺乏精確的干涉參數(shù)和影像外方位元素，在數(shù)據(jù)處理時遇到了干涉質(zhì)量、影像糾正等難題[7]，所獲得的DEM還沒有達到較高的精度。

本文在已有應(yīng)用機載InSAR數(shù)據(jù)生成DEM的技術(shù)流程[8]基礎(chǔ)上，引入精確干涉參數(shù)定標(biāo)和區(qū)域網(wǎng)平差處理技術(shù)，提出利用國產(chǎn)機載雙天線SAR數(shù)據(jù)生成丘陵地區(qū)大面積高精度DEM的方法。

1 試驗區(qū)概況及試驗設(shè)計

1.1 試驗區(qū)概況

江油試驗區(qū)位于四川盆地西北部及涪江上游，屬于四川省綿陽市;地貌以低山丘陵為主;海拔高度383～717 m。測區(qū)包含居民地、河流、耕地、林地、鐵路以及公路等地物類型，東西長59 km，南北寬4.6 km。經(jīng)緯度范圍介于 E104°30′00″～105°07′30″，N31°37′30″～ 31°40′之間。

1.2 機載InSAR系統(tǒng)及飛行參數(shù)

試驗采用X波段機載InSAR系統(tǒng)。該系統(tǒng)頻率為9.6 GHz，HH極化、乒乓模式;基線長度為2.33 m，脈沖重復(fù)頻率(PRF)為3 125 Hz，脈沖寬度為16 μs，平均功率約220 W;獲取的InSAR數(shù)據(jù)地面分辨率優(yōu)于0.5 m。飛機上安裝的高精度位置姿態(tài)測量系統(tǒng)(POS/AV510)能獲取傳感器的空間坐標(biāo)、旋轉(zhuǎn)角、速度及加速度等信息，參數(shù)如表1所示。

表1 飛行參數(shù)Tab.1 Flight parameters

1.3 野外數(shù)據(jù)獲取

為實現(xiàn)高精度區(qū)域網(wǎng)平差，便于成果精度檢驗，需要在試驗區(qū)精密布設(shè)角反射器控制點。根據(jù)任務(wù)需要及實際地形情況，在試驗區(qū)兩端及中部沿距離向分別布設(shè)3列(共28個)控制點。

另外，為實現(xiàn)高精度機載InSAR系統(tǒng)參數(shù)定標(biāo)，在測區(qū)中部北側(cè)沿距離向均勻布設(shè)9個定標(biāo)控制點。試驗區(qū)范圍及控制點布設(shè)如圖1所示。

圖1 試驗區(qū)范圍及地面控制點分布Fig.1 Location of study area and ground control points distribution

外業(yè)實地踏勘并選取控制點，記錄控制點點名和點位略圖，用GPS設(shè)備組網(wǎng)靜態(tài)觀測控制點。機載雷達飛行獲取數(shù)據(jù)時，在控制點上精確布設(shè)鋁制3面等腰直角角反射器(邊長35 cm)。差分GPS基準(zhǔn)站采用測區(qū)中央雙基站布放，利用高精度雙頻GPS接收機進行連續(xù)同步觀測。基站在飛機起飛前0.5 h開機，飛機降落后0.5 h關(guān)機。

2 機載雙天線InSAR獲取DEM原理

圖2為機載雙天線InSAR干涉測量成像示意圖。

圖2 機載雙天線InSAR干涉成像示意圖Fig.2 Imaging geometry of airborne dual－antenna InSAR

式(2)—(4)中:c為光速;τ為絕對時間延遲;ρ為距離向采樣率;j為P點方位向影像坐標(biāo);λ為電磁波波長;φ為該點的相位。

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

采用自主開發(fā)的機載InSAR地形制圖處理系統(tǒng)InSARmap軟件處理數(shù)據(jù)，其流程如圖3所示。

圖3 機載雙天線InSAR數(shù)據(jù)生成DEM處理流程Fig.3 Data processing flowchart of airborne dual－antenna InSAR DEM generation

3.1 干涉參數(shù)定標(biāo)

干涉參數(shù)定標(biāo)是高精度測圖的關(guān)鍵技術(shù)，也是測圖誤差控制的關(guān)鍵步驟。干涉參數(shù)定標(biāo)就是要對系統(tǒng)設(shè)備的相位偏移和干涉參數(shù)偏差進行嚴(yán)格校正，以提高生成DEM的精度[9]。

本文基于敏感度方程模型進行干涉參數(shù)定標(biāo)。首先通過目標(biāo)點高程方程對待定標(biāo)的干涉參數(shù)求微分，得到敏感度矩陣，建立高程誤差對干涉參數(shù)偏差的關(guān)系;再利用GPS測量定標(biāo)控制點高程，用最小二乘法求解各個干涉參數(shù)偏差估計值，進而校正相應(yīng)的干涉參數(shù)。試驗選用基線長度B、基線傾角α、絕對延遲時間τ和相位偏移量φoffset4個參數(shù)建立敏感度方程，即

高程差和干涉參數(shù)偏差的關(guān)系為[10]

式中:Δ為L×1高程誤差向量;F為L×N敏感度矩陣;△X為N×1待定標(biāo)的干涉參數(shù)偏差向量;L為定標(biāo)控制點數(shù)目;N為待定標(biāo)參數(shù)個數(shù)，此處N=4。采用迭代方法進行干涉參數(shù)定標(biāo)，用初始的干涉參數(shù)和相位生成DEM，通過解算敏感度方程得到各干涉參數(shù)偏差值，校正相應(yīng)干涉參數(shù)，并重新計算DEM，直到2次計算的DEM均方差小于0.01為止。

3.2 干涉處理

干涉處理步驟包括復(fù)圖像配準(zhǔn)、干涉紋圖生成、干涉紋圖濾波以及相位解纏等。本研究采用基于快速傅里葉變換(fast fourier transform，F(xiàn)FT)的復(fù)相關(guān)法[8]進行復(fù)圖像精配準(zhǔn)，配準(zhǔn)圖像對干涉處理后生成干涉紋圖和相干圖，并采用中值濾波和Goldstein濾波相結(jié)合的濾波方法對復(fù)干涉紋圖進行濾波處理，采用最小費用流法對濾波后的干涉紋圖進行相位解纏，獲得影像上每點的相位數(shù)據(jù)。通過InSAR-map軟件中的干涉處理模塊可實現(xiàn)整條帶數(shù)據(jù)的自動批量干涉處理，能節(jié)省大量處理時間。

圖4 區(qū)域網(wǎng)連接點分布示意圖Fig.4 Control points distribution in bundle adjustment

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

由于POS系統(tǒng)測量的航偏角及雷達系統(tǒng)提供的多普勒中心頻率精度不夠高，在測區(qū)稀少控制點條件下僅依靠POS數(shù)據(jù)及干涉定標(biāo)參數(shù)難以實現(xiàn)大面積、多航帶SAR影像的高精度校正。區(qū)域網(wǎng)平差技術(shù)可以建立影像之間的約束關(guān)系，通過調(diào)整和精化部分外方位元素，實現(xiàn)稀少控制點條件下的多航帶SAR影像的高精度定位[11]。

本文基于 Range-Doppler-Phase(RDP)模型[11]建立載機位置、像點及地面點之間的關(guān)系，確定平差參數(shù)系數(shù)矩陣，建立高精度三維模型誤差方程組。利用地面控制點和連接點數(shù)據(jù)、定標(biāo)數(shù)據(jù)、POS數(shù)據(jù)和相位信息，按最小二乘原理求解各景影像準(zhǔn)確的外方位元素(平差參數(shù))和連接點的地面坐標(biāo)。確定的平差參數(shù)包括解纏相位φ偏移量、多普勒中心頻率fDC及偏航角θy。由于高程和平面坐標(biāo)對平差參數(shù)敏感程度不同，將高程和平面平差參數(shù)分開迭代求解。首先求解高程平差參數(shù)(解纏相位φ偏移量改正)，再將改正的高程平差參數(shù)代入平面誤差方程，求解平面平差參數(shù)(多普勒中心頻率fDC偏移量，偏航角θy改正數(shù))。

在試驗區(qū)4個條帶內(nèi)及條帶間影像重疊區(qū)域選取6 254個連接點，組成區(qū)域網(wǎng)，如圖4所示。

在試驗區(qū)布設(shè)的37個角反射器控制點(其中28個測區(qū)控制點，9個定標(biāo)點)中，21個(12個測區(qū)控制點，9個定標(biāo)點)用于參與平差計算，16個作為檢查點(圖1)。

3.4 相高轉(zhuǎn)換

用平差計算方法獲得各景影像準(zhǔn)確的外方位元素后，就可進行相高轉(zhuǎn)換，以獲得影像上各點的高程。由式(1)—(4)可知，獲取像點的高程需已知該點的航高H、基線傾角α、基線長度B、絕對時間延遲τ和相位φ。其中，H，α，B和τ延遲由參數(shù)定標(biāo)獲得，相位由干涉處理數(shù)據(jù)及區(qū)域網(wǎng)平差參數(shù)共同確定。通過相高轉(zhuǎn)換獲得影像上各點的高程。

3.5 地理編碼

為便于應(yīng)用，需將原始斜距幾何下的InSAR高度數(shù)據(jù)投影到某一參考系統(tǒng)下。本研究基于正側(cè)視模型，利用載機飛行數(shù)據(jù)信息，對高度圖像進行地理編碼。由POS數(shù)據(jù)獲取方位向上每掃描行數(shù)據(jù)的載機位置(XS，YS)，坐標(biāo)轉(zhuǎn)換得到其高斯坐標(biāo)(xs，ys)。為了能更準(zhǔn)確地確定目標(biāo)點位置，將通過平差計算得到的多普勒中心頻率fDC引入斜視角β[1]，即

由式(1)求出側(cè)視角θ。設(shè)θy為航偏角(由區(qū)域網(wǎng)平差結(jié)果得到)，利用幾何關(guān)系得到該像點在高斯平面上的坐標(biāo)(xp，yp)[10]，即

通過以上過程便可獲得每一像點的地理坐標(biāo)，再進行重采樣，將影像的高程h(i，j)賦到重采樣后的像點上，得到地理編碼的DEM。

3.6 影像鑲嵌

本文將試驗區(qū)分為4個條帶進行數(shù)據(jù)獲取。條帶間重疊率≥40%;條帶內(nèi)相鄰影像重疊率為20%左右。在生成各景DEM之后，需在統(tǒng)一的坐標(biāo)系統(tǒng)(WGS-84坐標(biāo)系高斯-克呂格投影3°分帶，1985國家高程基準(zhǔn))鑲嵌DEM。具體步驟如下:

1)對于非重疊區(qū)域的像素，其高度值由覆蓋該像素的DEM直接賦與。

2)對于二度重疊區(qū)域的像素，其高度值等于覆蓋該像素的2個高度值的平均值。

3)對于重疊度≥3區(qū)域的像素，采用加權(quán)平均的方法[3]確定，將該點上的所有高度值進行排序，把偏離中值超出閾值(設(shè)為σ=±1.2 m)的高度值從序列中剔除，剩余高度值按其干涉相干度值進行定權(quán)并進行加權(quán)平均，得到該點的高度值。

4 試驗結(jié)果及精度分析

4.1 試驗結(jié)果

基于本文提出的高精度DEM的生成技術(shù)流程，對試驗區(qū)的1個架次、4條帶及76景SLC影像進行地形制圖處理，利用自主開發(fā)的機載InSAR地形制圖處理系統(tǒng)軟件快速實現(xiàn)上述處理過程，生成了覆蓋超過500 km2的一整幅高精度DEM圖像(圖5)。

圖5 處理生成的DEMFig.5 Generated DEM from test data

4.2 精度分析

為了評價上述流程所生成的DEM精度，本研究將精確布設(shè)在試驗區(qū)但未參與平差的16個角反射器(圖1中◆所示)作為檢查點，將其平面坐標(biāo)及高程的實測值與生成值進行比較，評定所生成DEM的內(nèi)部精度(圖6)。

利用在試驗區(qū)實測的185個GPS檢查點進行精度檢核。由于這些檢查點在SAR影像上的位置特征不明顯，刺點將會產(chǎn)生較大的位置偏差，所以未統(tǒng)計其平面誤差，僅檢核高程。發(fā)現(xiàn)其中有4個GPS點落入到陰影區(qū)域中，在后續(xù)比較中排除這4個點。GPS檢查點多位于開闊區(qū)域且其測量精度較高，根據(jù)GPS檢查點的平面坐標(biāo)將其與DEM數(shù)據(jù)疊加，用InSAR DEM高程值HDEM減去檢查點高程值HGPS得到外部檢核的高程誤差值，如圖7所示。

圖6 X，Y方向的平面誤差分布Fig.6 Error distribution of X and Y direction

圖7 高程誤差平面分布及其統(tǒng)計直方圖Fig.7 Distribution and statistical histogram of elevation error between GPSs and InSAR DEM

表2列出了16個角反射器檢查點(JC01—JC16)、181個地面 GPS高程檢查點(GZ01—GZ181)對DEM的檢測精度，并對檢查點的平面點位中誤差和高程中誤差分別進行了統(tǒng)計分析。

表2 DEM精度統(tǒng)計Tab.2 Precision statistics of DEM (m)

可以看出，圖6中X方向誤差和Y方向誤差在0值附近分布，其均值分別為0.15 m和0.01 m;在圖7左側(cè)，高程誤差在正負(fù)區(qū)間均有分布，且大部分高程誤差集中在±1 m之間，其均值為0.057 m，標(biāo)準(zhǔn)差為±0.547 m;在圖7右側(cè)，高程誤差頻率分布曲線略呈左偏態(tài)，而非正態(tài)分布。高程誤差在正區(qū)間0.2處兩側(cè)分布的點明顯偏多，即InSAR生成DEM略高于實測GPS檢查點的高程。這主要是由于X波段微波的表面穿透性較弱，InSAR獲取的并不是裸露地形表面的高度[2]。

通過精密布設(shè)的角反射器檢查點的檢核，DEM的平面點位中誤差為±1.188 m，高程中誤差為±0.508 m;動態(tài)GPS測量檢查點的高程中誤差也達到了±0.550 m?？梢姡瑧?yīng)用機載InSAR數(shù)據(jù)生成的DEM的精度達到了1∶1萬丘陵地區(qū)一級高程中誤差(表3)的要求。

表3 1∶1 萬 DEM 技術(shù)指標(biāo)[12]Tab.3 Technical indicator of 1∶10 000 DEM[12] (m)

5 結(jié)論

本文提出了利用國產(chǎn)機載雙天線InSAR數(shù)據(jù)生成大面積高精度DEM的技術(shù)流程?；谧灾鏖_發(fā)的機載InSAR地形測圖軟件處理了高分辨率X波段機載雙天線SAR數(shù)據(jù)，獲得了大面積高精度DEM，并對生成的DEM進行精度分析，得出如下結(jié)論:

1)對影響高程精度的基線長度、基線傾角、相位偏置等參數(shù)，利用布放合理、精密測量的角反射器作為地面控制點進行外定標(biāo)處理，這對于保證制圖精度非常重要。

2)區(qū)域網(wǎng)平差不僅可以減少地面控制點數(shù)量，還能提高連接點坐標(biāo)的精度和區(qū)域網(wǎng)整體性，是解決面積大、控制點稀少條件下地形測圖的關(guān)鍵技術(shù)。

3)對于丘陵地區(qū)，試驗證明了基于國產(chǎn)機載雙天線SAR數(shù)據(jù)的DEM生成技術(shù)可滿足1∶1萬比例尺的地形圖制圖精度要求;機載InSAR技術(shù)可作為復(fù)雜地區(qū)地形測圖制取的一種技術(shù)手段。

志謝:感謝中國科學(xué)院電子學(xué)研究所提供了機載雙天線SAR數(shù)據(jù)。

[1]孫中昶，郭華東，李新武.機載雙天線InSAR數(shù)據(jù)生成高精度DEM 的誤差分析[J].高技術(shù)通訊，2012，22(2):171-179.Sun Z C，Guo H D，Li X W.Error analysis of high precision DEM generated from airborne dual-antenna interferometric SAR data[J].Chinese High Technology Letters，2012，22(2):171-179.

[2]Li X P，Baker A B，Hutt T.Accuracy of airborne IFSAR mapping[C]//Proceedings of the American Society of Photogrammetry and Remote Sensing，XXII International Congress，Washington，USA.2002.

[3]Wimmer C，Siegmund R，Schwabisch M，et al.Generation of high precision DEMs of the Wadden Sea with airborne interferometric SAR[J].IEEE Transactions on Geoscience and Remote Sensing，2000，38(5):2234-2245.

[4]Mercer B.National and regional scale DEMs created from airborne InSAR[C]//ProcPIA，2007，36(3):W49A.

[5]黃國滿，張繼賢，趙 爭，等.機載干涉SAR測繪制圖應(yīng)用系統(tǒng)研究[J].測繪學(xué)報，2008，37(3):277-279.Huang G M，Zhang J X，Zhao Z，et al.Research on airborne SAR interferometry mapping system[J].Acta Geodaetica et Cartographica Sinica，2008，37(3):277-279.

[6]Xiang M S，Wu Y R，Li S E，et al.Introduction on an experimental airborne InSAR system[C]//Proceedings of the 2005 IEEE International Geoscience and Remote Sensing Symposium，Seoul Korea:IEEE，2005(7):4809-4812.

[7]劉艷華，趙 爭，黃國滿.機載InSAR數(shù)據(jù)自動生成DEM技術(shù)及其在內(nèi)蒙古豐鎮(zhèn)地區(qū)的應(yīng)用[J].地學(xué)前緣，2006，13(3):104-107.Liu Y H，Zhao Z，Huang G M.Automatic DEM generation technology by airborne InSAR data and its application to Fengzhen，Inner Mongolia[J].Earth Science Frontiers，2006，13(3):104-107.

[8]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.

[9]楊懷寧，郭華東，韓春明.機載InSAR敏感度方程定標(biāo)限制條件的仿真實驗[J].高技術(shù)通訊，2010，20(10):1049-1054.Yang H N，Guo H D，Han C M.Imitation experiment on the parameter requirements of the sensitivity equation-based calibration for airborne InSAR[J].Chinese High Technology Letters，2010，20(10):1049-1054.

[10]張 薇，向茂生，吳一戎.基于正側(cè)視模型的機載雙天線干涉SAR 外定標(biāo)方法[J].遙感技術(shù)與應(yīng)用，2008，23(3):346-350.Zhang W，Xiang M S，Wu Y R.Studies on outside calibration method based on the boresight model for dual-antenna airborne interferometric SAR[J].Remote Sensing Technology and Application，2008，23(3):346-350.

[11]馬 婧，尤紅建，胡東輝.FLeberl模型與干涉測量模型相結(jié)合的InSAR影像區(qū)域網(wǎng)平差[J].紅外與毫米波學(xué)報，2012，31(3):271-276.Ma J，You H J，Hu D H.Block adjustment of InSAR images based on the combination of FLeberl and interferometric models[J].Journal of Infrared and Millimeter Waves，2012，31(3):271-276.

[12]國家測繪局測繪標(biāo)準(zhǔn)化研究所.CH/T1008—2001基礎(chǔ)地理信息數(shù)字產(chǎn)品1∶10 000，1∶50 000 數(shù)字高程模型[S].北京:國家測繪局，2001.National Institute of Surveying and Mapping Standardization.CH/T1008—2001digital products of fundamental geographic information 1∶10 000，1∶50 000 digital elevation models[S].Beijing:State Bureau of Surveying and Mapping，2001.