利用參數(shù)獨立分解的星載SAR干涉測量檢校方法

范 軍，李 濤，左小清，陳乾福，張 祥，祿 競

1. 昆明理工大學國土資源工程學院，云南 昆明 650093； 2. 自然資源部國土衛(wèi)星遙感應(yīng)用中心，北京 100048

合成孔徑雷達干涉(interferometric synthetic aperture radar,InSAR)技術(shù)具有全天時、全天候、高精度等優(yōu)勢[1-2]，在地形測繪中具有良好的應(yīng)用前景。然而，高精度DEM的獲取，一方面取決于InSAR平臺指標的設(shè)計，包括載荷設(shè)計和姿軌控制等，另一方面則取決于幾何及干涉參數(shù)的檢校工作[3]。幾何檢校通過校準獲取的幾何參數(shù)來提升SAR影像的幾何定位精度，而干涉測量檢校技術(shù)是通過校準獲取的干涉參數(shù)來提升干涉結(jié)果的高程精度，二者共同決定了全球DEM的精度水平[4]。幾何檢校技術(shù)一般以距離-多普勒(Range-Doppler,R-D)模型為基礎(chǔ)解算方位向時間延遲和距離向時間延遲參數(shù)[5]。ALOS_PALSAR衛(wèi)星使用幾何檢校方法后，條帶模式的幾何定位精度達9.7 m[6]；Radarsat-2衛(wèi)星利用角反射器進行幾何檢校后，獲得了5.74 m的平面定位精度[7]；TerraSAR-X利用角反射器完成了系列衛(wèi)星的高精度幾何檢校，實現(xiàn)了0.3 m的平面定位結(jié)果[8-9]；2018年，文獻[10]針對TerraSAR-X/TanDEM-X影像和GF-3影像進行了幾何檢校，檢校后的TerraSAR-X/TanDEM-X影像幾何定位精度優(yōu)于3 m，GF-3影像幾何定位精度優(yōu)于7.5 m。因此，幾何檢校方法已經(jīng)在工程和科研上取得了成功應(yīng)用。干涉測量檢校技術(shù)較為復(fù)雜，國內(nèi)外已開展了大量研究工作，提出了很多干涉測量檢校方法。1999年，文獻[11]分析了影響干涉測高的系統(tǒng)參數(shù)，并對誤差源進行了敏感度分析，在平地幾何關(guān)系下推導出基于敏感度方程的機載干涉測量檢校方法；2001年，文獻[12]在此基礎(chǔ)上，推導出斜視下基于干涉敏感度方程的機載干涉測量檢校方法，通過仿真數(shù)據(jù)證明了不同敏感度矩陣條件數(shù)對干涉測量檢校產(chǎn)生不同的影響，但是沒有針對定標器布放問題進行研究；2003年，文獻[13]針對干涉敏感度方程的推導引入了三維重建模型，并提出一種使敏感度矩陣條件數(shù)最小化的定標器布放算法；2010年，文獻[14]提出一種考慮干涉相位偏置、基線長度和基線水平角參數(shù)的機載干涉定標模型，但是未考慮控制點的平面信息；文獻[15]在此基礎(chǔ)上提出利用控制點三維信息的機載雙天線SAR干涉定標方法，驗證了定標后DEM的平面精度和高程精度，但是并沒有獨立解算每一項干涉參數(shù)的修正量；2011年，文獻[16]基于三維重建模型的基礎(chǔ)上，提出了一種分布式衛(wèi)星SAR干涉測高模型，通過仿真數(shù)據(jù)針對基線矢量和相位誤差特性進行了分析；2014年，文獻[17]提出了星載InSAR立體基線的檢校方法，通過仿真數(shù)據(jù)表明該方法能完成高精度基線矢量檢校，但是缺少試驗數(shù)據(jù)的支撐；2016年，文獻[18]利用TanDEM-X數(shù)據(jù)驗證了基于InSAR三維重建模型、敏感度方程的干涉測量檢校方法，證明了星載InSAR干涉測量檢校方法的通用性?，F(xiàn)階段，商業(yè)軟件GAMMA在反演InSAR高程時只利用控制點的高程信息完成了對干涉參數(shù)的精估計，沒有考慮平面精度。另外，軟件中并未獨立分解干涉誤差源，沒有達到高精度的干涉測量檢校要求。到目前為止，從機載干涉測量檢校模型發(fā)展到星載中，盡管已經(jīng)開展大量研究工作，但是在干涉過程中，多數(shù)參數(shù)的誤差具有極強的耦合性，聯(lián)合解算卻無法進行準確的區(qū)分，從而導致干涉參數(shù)修正量不夠精確。本文針對這一問題，提出利用參數(shù)獨立分解的干涉測量檢校方法，在保證幾何定位精度的前提下，通過4對TanDEM-X/TerraSAR-X數(shù)據(jù)處理驗證了星載SAR干涉測量檢校方法的正確性。

1 InSAR三維重建模型

InSAR三維重建模型基本原理是通過聯(lián)立距離方程、多普勒方程和干涉相位方程獲取地面目標點的三維坐標。圖1所示為地心坐標系下InSAR測量的幾何關(guān)系，其中S1、S2為主輔天線相位中心位置矢量，R1、R2表示主輔天線相位中心到地面點的位置矢量，R1、R2分別為主輔天線相位中心到地面點P的距離，P為地面目標點位置矢量，v為速度矢量，B為基線矢量，fdop為多普勒頻率，λ為雷達波長，φ為絕對干涉相位，φ為纏繞干涉相位，k為整周期數(shù)，Q為天線收發(fā)模式，Q=1表示“一發(fā)雙收”模式，Q=2表示“乒乓”模式。

本文研究主要集中于星載InSAR交軌干涉測量(cross track interferometry,XTI)模式[13]，則三維重建幾何關(guān)系表示為

P=S1+R1

(1)

(2)

(3)

圖1 星載InSAR在地形測繪中的幾何關(guān)系Fig.1 Geometric relationship of InSAR in topographic mapping

根據(jù)干涉SAR的幾何關(guān)系，目標點坐標在地心空間直角坐標系下表示為

P=S1+R1=S1+R1r

(4)

式中，R1為視向量；r為單位視向量。根據(jù)式(4)可知，對目標點坐標的求解一般采用視向量分解法[13]，即將視向量轉(zhuǎn)換為對單位視向量的求解。單位視向量需要在移動坐標系VPQ中分解為V、P和Q3個坐標軸方向[15,18]。VPQ坐標系是以主天線相位中心為坐標原點，速度矢量方向為V軸，速度和基線叉乘方向為Q軸，V軸、P軸和Q軸滿足右手坐標系關(guān)系。rv、rp和rq分別表示單位視向量r在V軸、P軸和Q軸上的分量，即

(5)

視向量分解法將視向量從移動坐標系向地心空間直角坐標系轉(zhuǎn)換，從而獲取目標點的三維坐標[19]。

P=S1+R1=S1+R1Tvpqrvpq=

(6)

式中，Tvpq是將視向量轉(zhuǎn)換為地心空間直角坐標系的變換矩陣；bv為基線在速度方向的分量；bpv為基線在垂直于速度方向的分量。

由三維重建模型表達式可知，影響干涉SAR測高精度的主要因素有：主星位置、主星速度、主星斜距、多普勒頻率、干涉相位誤差、基線矢量誤差等。主星位置、主星速度、主星斜距、多普勒頻率等參數(shù)在幾何檢校中校正，以提升SAR影像的平面精度[20-21]。而干涉相位誤差和基線矢量誤差采用本文提出的干涉測量檢校方法進行校正，以提升干涉結(jié)果的高程精度，最終獲取高精度的DEM產(chǎn)品。

2 干涉敏感度分析

本文在保證幾何檢校精度的前提下，本文針對干涉相位和基線矢量參數(shù)引起的誤差進行敏感度分析，即定量分析干涉參數(shù)的變化對高程誤差的影響，同時為后續(xù)試驗部分提供理論參考。

2.1 干涉相位誤差敏感度分析

干涉相位誤差一般由隨機誤差、系統(tǒng)誤差和絕對相位偏差組成。其中隨機誤差是由InSAR系統(tǒng)各項相干源造成，對應(yīng)的相位呈隨機分布特性，在本文中不參與檢校；系統(tǒng)誤差成分主要由硬件設(shè)備隨工作溫度變化的剩余漂移相位造成[22-23]，呈緩慢特性，在本文中不參與檢校；絕對相位偏差是指絕對干涉相位與解纏相位和平地相位之間相差的固定常數(shù)，因此，干涉相位誤差的主要成分為絕對相位偏差。干涉相位誤差Δφ引起的DEM高程誤差Δh為

(7)

式中，HoA為高程模糊度。TanDEM-X采用了兩期的HoA設(shè)置，第1期為40～55 m，第2期為35 m，即2π相位變化引的HoA在35～55 m之間[3,24]。在HoA=35 m的情況下，10°的相位誤差引入高程誤差約為0.97 m，在HoA=55 m的情況下，10°的相位誤差引入高程誤差約為1.54 m。因此，約10°的相位誤差引起的高程誤差變化范圍在0.97～1.54 m之間。設(shè)定干涉相位誤差范圍為[0°，10°]，結(jié)合三維重建模型表達式評定其對高程測量結(jié)果的影響，如圖2所示，10°相位誤差引起高程誤差約為1.42 m，隨著干涉相位誤差的增加，高程測量誤差呈線性遞增的趨勢。如果要進行1∶50 000比例尺地形圖測繪，那么丘陵地的高程精度需要達到5 m，相位誤差的影響已經(jīng)約占28.4%，因此，干涉相位誤差是InSAR地形測繪中的重要誤差因子。

圖2 干涉相位誤差敏感度分析Fig.2 Sensitivity analysis of interferometric phase

2.2 基線誤差敏感度分析

在星載InSAR地形測繪中，基線矢量直接影響著高程模糊度、高程誤差的傳遞系數(shù)等諸多因素。基線是一個三維矢量，為了更方便的表達基線參數(shù)，一般采用TCN坐標系，以主天線相位中心為坐標原點，N軸由主天線相位中心指向地心，N軸與速度矢量叉乘方向為C軸正方向，T、C、N三軸構(gòu)成右手坐標關(guān)系[18]。本文將基線矢量分解為TCN坐標系下的基線初值和與時間相關(guān)的基線速率，考慮到交軌干涉過程中，所有的計算都是投影到主影像的零多普勒平面中進行的，因此沿軌基線恒為零，即Bt=0[25]。

(8)

式中，Bt、Bc和Bn分別表示TCN坐標系的基線分量；bc0和bcv表示基線Bc在C方向分解的基線初值和基線速率；bn0和bnv表示基線Bn在N方向分解的基線初值和基線速率。由于三維重建模型是在地心空間直角坐標系表示InSAR的幾何關(guān)系，因此，本文需要將基線由TCN坐標系轉(zhuǎn)換為地心空間直角坐標系，即

(9)

式中，Ttcn為TCN與空間直角坐標系之間的變換矩陣。由于基線誤差對DEM高程精度影響較為復(fù)雜，本文在TCN坐標系下分析基線分量bc0、bcv、bn0和bnv對高程精度的影響，在地心空間直角坐標系下分析基線長度B對高程精度的影響。

若對C、N方向基線初值添加1 mm誤差，如圖3(a)、(b)所示，C方向引起1.2 m的高程誤差，N方向引起-1.6 m的高程誤差；若對C、N方向基線速率添加0.1 mm/s誤差，如圖3(c)、(d)所示，C方向引起0.5 m高程誤差，N方向-0.6 m的高程誤差?？梢钥闯?，數(shù)毫米級的基線初值誤差可引起米級的高程誤差，數(shù)亞毫米級的基線速率誤差可引起米級的高程誤差。此外，與基線初值誤差相比，基線速率誤差對高程測量精度影響更為敏感，與C方向基線矢量誤差相比，N方向?qū)Ω叱虦y量精度影響更為敏感。

若對基線長度添加1 cm誤差，如圖3(e)所示，則引入高程誤差約為1.61 m。如果要進行1∶50 000比例尺地形圖測繪，那么丘陵地的高程精度需要達到5 m，基線誤差的影響已經(jīng)約占32%，綜上所述，基線矢量是InSAR地形測繪中一個重要的誤差因子。

3 InSAR干涉測量檢校模型

本文方法基于干涉相位誤差與基線矢量誤差到高程誤差的不同傳遞特性進行建模。通過敏感度分析可知，干涉相位引入的高程誤差為1.42 m，基線在方位向和距離向引入高程誤差為1.61 m，由此可見，干涉相位誤差和基線矢量誤差對地形測繪的影響不可忽視。因此，本文提出一種利用參數(shù)獨立分解的干涉測量檢校方法。首先解算干涉相位誤差，隨后解算基線矢量誤差。為了保證兩種干涉參數(shù)的解算方法完全獨立，本文采用計算絕對干涉相位與解纏相位和平地相位之間固定偏差的方法，來獲取干涉相位誤差。首先根據(jù)地面控制點來求取每一個點的絕對干涉相位，然后依據(jù)該局部絕對干涉相位對解纏相位和平地相位進行補償，最后求取整景相位數(shù)據(jù)應(yīng)該補償?shù)墓潭ㄆ?，即可獲取每一景影像的干涉相位誤差。

圖3 基線誤差敏感度分析Fig.3 Sensitivity analysis of baseline errors

(10)

式中，R1、R2表示主輔影像的斜距；φunw為解纏相位；φfla為濾波前的平地相位；Δφ為干涉相位誤差。

對基線矢量參數(shù)進行干涉測量檢校時，首先從InSAR三維重建模表達式出發(fā)，然后建立敏感度方程，最后采用平差模型基于最小二乘法準則解算干涉參數(shù)的修正量。為了提高解算精度，本文采用迭代方式，將迭代前后坐標差值的均方根定義為檢校精度，當檢校精度小于指定閾值(本文取ε=0.05)時，則迭代收斂，最終獲取高精度的干涉參數(shù)，從而保證生成高精度DEM產(chǎn)品，干涉測量檢校流程如圖4所示。

利用三維重建模型表達式在地心空間直角坐標系推導出3個坐標軸方向的三維重建方程，形式如下[19]

(11)

圖4 干涉測量檢校流程Fig.4 Flow chart of interferometric calibration

根據(jù)三維重建方程對基線矢量計算偏導數(shù)，建立干涉敏感度方程，形式如下

(12)

建立誤差方程式

(13)

在干涉測量檢校過程中，由于不同觀測條件下，觀測值的精度往往不同，因而誤差也不同，不同控制點位置的數(shù)據(jù)質(zhì)量也不同。由于相干性系數(shù)是評價SAR影像質(zhì)量的重要標準，因此對控制點的相干性系數(shù)定權(quán)[22]，會在后續(xù)研究中考慮其他定權(quán)方法。此時目標函數(shù)為

Jmin=VTWV

(14)

式中，W是對角權(quán)陣

(15)

式中，wi(i=1,2,…,n)為每個控制點的相干性系數(shù)。

給定干涉參數(shù)的初值X0，解算干涉參數(shù)的修正量，得

ΔX=(ATWA)-1ATWl

(16)

從而得到檢校后的干涉參數(shù)

X=X0+ΔX

(17)

4 試驗與結(jié)果分析

4.1 試驗數(shù)據(jù)

試驗數(shù)據(jù)中心區(qū)域位于陜西渭南市，中心經(jīng)緯度為109°34′25″E、35°01′48″N，東西向范圍約達63 km，南北向范圍約達108 km。區(qū)域內(nèi)兼有平地和山地，海拔156～1338 m，其地理位置如圖5(a)所示。選取了83個地面控制點，通過采用SF-3050星站差分GNSS接收機實測，坐標系為大地坐標系，平面測量誤差為15 cm，可滿足1∶25 000比例尺地形圖測繪精度要求；高程測量誤差為20 cm，引入的干涉相位誤差小于1.12°，可用于進行高精度測量檢校。本文選取的地面控制點主要分布在SAR影像中具有明顯成像特征的區(qū)域，如道路交叉口、房屋拐角點等[26]，詳細點位分布如圖5中綠色三角形所示。本文在使用ICESat激光測高數(shù)據(jù)前進行了數(shù)據(jù)篩選[18,27]，將一些受云層、樹木、建筑物等遮擋引起高程誤差達到幾十米甚至到幾百米的數(shù)據(jù)剔除，最終共選擇激光點7480個，詳細點位分布如圖5中紅色圓點所示，圖5(b)為檢校后DEM結(jié)果。

圖5 試驗區(qū)域地理位置及DEM的分布圖Fig.5 Geographical location and DEM distribution map of study area

試驗數(shù)據(jù)為4對TanDEM-X/TerraSAR-X干涉數(shù)據(jù)，如圖5中白色線框所示，獲取時間為2013年9月3日與2013年9月25日。選擇90 m格網(wǎng)分辨率SRTM作為外部DEM數(shù)據(jù)。試驗影像基本參數(shù)如表1所示。

表1 試驗影像參數(shù)

4.2 結(jié)果與分析

在保證幾何精度的條件下，應(yīng)用三維重建模型解算控制點高程，與已知控制點高程進行精度評價，如表2所示。

表2 檢校前干涉結(jié)果的高程精度

Tab.2 Elevation accuracy of interferometric results before interferometric calibration

影像編號影像類型高程精度/m696_8檢校景20.016696_9檢校景21.591696_13檢校景9.444696_14檢校景8.925

由表2可知，不同軌道干涉圖結(jié)果的高程精度存在較大差異。例如同軌696_13、696_14的高程精度低于9.44 m，而另一同軌696_8、696_9的高程精度低于21.59 m。高程存在較大誤差的原因在于，首先，由于地形的限制，導致控制點本身精度降低；其次，存在幾何檢校的殘余誤差，導致在干涉測量檢校中引入幾何誤差。

考慮到對基線誤差解算會引入干涉相位誤差，再對干涉相位誤差解算反而增加高程誤差。本文首先對干涉相位進行檢校，將每一景影像的干涉相位誤差分解后，根據(jù)三維重建模型解算控制點高程，與已知控制點高程進行精度評價，結(jié)果如表3所示。

表3 干涉相位檢校后干涉結(jié)果的高程精度

Tab.3 Elevation accuracy of interferometric results after interferometric phase calibration

影像編號影像類型干涉相位誤差/(°)高程精度/m696_8檢校景-120.7682.057696_9檢校景-127.0592.848696_13檢校景-50.0992.726696_14檢校景-52.1052.312

由圖2干涉相位誤差敏感度線性關(guān)系可知，以影像696_8為例，當干涉相位誤差增加120.768°時，則引起約17.20 m的高程誤差敏感度。本試驗中，120.768°的相位誤差引起高程誤差變化量是17.96 m，與敏感度分析結(jié)果基本一致。由表3可知，通過對干涉相位誤差檢校后，所有干涉結(jié)果的高程精度均優(yōu)于2.85 m。在此基礎(chǔ)上，對基線矢量展開檢校，精度評價結(jié)果如表4所示。

表4 基線矢量檢校后干涉結(jié)果的高程精度

Tab.4 Elevation accuracy of interferometric results after baseline vector calibration

影像編號影像類型基線誤差修正量Δbc0/mΔbn0/mΔbcv/(m/s)Δbnv/(m/s)高程精度/m696_8檢校景 0.0597-0.0473 -0.0092 0.0072 1.469696_13檢校景-0.05520.04200.0057-0.0043 2.535

本文將基線矢量分解為TCN方向的基線初值和基線速率，其共同影響干涉結(jié)果的高程精度。以影像696_8為例，C、N方向基線初值誤差為0.059 7 m、-0.047 3 m，基線速率誤差為-0.009 2 m/s、0.007 2 m/s。由式(18)計算得到沿方位向時間變化的基線誤差在[0.029，0.122]之間，基于檢校前后的高程精度，發(fā)現(xiàn)基線誤差引起高程誤差為19.16 m。根據(jù)圖3(e)基線誤差敏感度線性關(guān)系，將基線誤差增加到0.122 m時，高程誤差已經(jīng)達到19.63 m，與檢校結(jié)果基本一致。通過表4可以看出，對基線矢量檢校后，所有干涉結(jié)果的高程精度均優(yōu)于2.54 m。綜上所述，干涉測量檢校技術(shù)可提高干涉結(jié)果的高程精度。

(18)

為了驗證同軌影像之間檢校參數(shù)的適用性，本文進行檢校補償驗證。分別將影像696_8、696_13檢校結(jié)果補償于同軌影像696_9、696_14，最終對驗證景數(shù)據(jù)進行精度評價，如表5所示。

表5 驗證景數(shù)據(jù)精度評價

由表5可知，利用檢校參數(shù)對驗證景補償后，無控制點的干涉結(jié)果高程精度均優(yōu)于2.23 m，滿足山地地區(qū)1∶25 000地形圖測繪要求。

本文通過提出了利用參數(shù)獨立分解的干涉測量檢校方法，分解了星載InSAR系統(tǒng)中的主要誤差源，同時這些因素也是引起DEM產(chǎn)品誤差的主要來源。本文選用冰、云和陸地高程衛(wèi)星(Ice,Cloud,and land Elevation Satellite,ICESat)[27-28]數(shù)據(jù)分別對未檢校和檢校后DEM的絕對高程精度進行了評估。如表6所示，影像696_8和696_13位于平地地區(qū)，高程精度提升到1.21 m；影像696_9和696_14位于山地地區(qū)，高程精度提升到3.11 m，其滿足了我國1∶25 000比例尺平地地區(qū)1.5 m、山地地區(qū)4 m的測繪精度。

綜上所述，本文提出的干涉檢校方法均可以提高DEM產(chǎn)品精度。但是檢校后平地地區(qū)獲取DEM產(chǎn)品的精度始終優(yōu)于山地地區(qū)，這主要是因為山地地形起伏較大，且植被茂盛，相干性較差，對干涉結(jié)果產(chǎn)生一定的影響，從而導致該區(qū)域的高程誤差較大，誤差分布比較離散。

表6 DEM精度評價

5 結(jié) 論

目前，隨著國產(chǎn)SAR衛(wèi)星的陸續(xù)發(fā)射，迫切需要對星載SAR干涉測量檢校技術(shù)深入研究，提高國產(chǎn)星載SAR系統(tǒng)獲取高精度DEM的能力。由于國內(nèi)外InSAR干涉測量檢校工作一般從機載工作開始研究，以此作為平臺逐步將干涉檢校技術(shù)擴展到星載平臺。本文針對傳統(tǒng)干涉測量檢校模型中存在參數(shù)耦合性過強的問題，提出了一種利用參數(shù)獨立分解的星載SAR干涉測量檢校方法，采用兩種獨立檢校算法分別對相位誤差和基線矢量誤差進行校正，通過獲取4對TanDEM-X/TerraSAR-X干涉數(shù)據(jù)開展干涉測量檢校試驗。結(jié)果表明，對于該試驗數(shù)據(jù)，干涉測量檢校后干涉結(jié)果的高程精度均優(yōu)于2.54 m。通過ICESat數(shù)據(jù)分別對未檢校、檢校后的DEM產(chǎn)品精度進行評估，檢校后平地地區(qū)DEM產(chǎn)品的絕對高程精度優(yōu)于1.21 m，山地地區(qū)DEM產(chǎn)品的絕對高程精度優(yōu)于3.11 m。由此可見，本文提出的干涉測量檢校方法能提高DEM產(chǎn)品精度，證明本文的方法是有效可行的。本文工作仍存在一些問題有待深入研究，例如未對控制點選取的數(shù)量以及分布情況進行探討，特別是如何減少檢校點的數(shù)量，實現(xiàn)稀少檢校點的干涉測量檢校等。另外，由于缺少高精度角反射器，本文未能對衛(wèi)星在軌檢校問題進行驗證，后續(xù)研究工作將主要針對這兩方面開展。

致謝：感謝由德國宇航中心提供TerraSAR-X、TanDEM-X數(shù)據(jù)，美國航空航天局提供ICESat數(shù)據(jù)。