聯(lián)合對(duì)飛數(shù)據(jù)的單控制點(diǎn)機(jī)載干涉SAR定標(biāo)算法

汪丙南，向茂生,2，蔣 帥,2，付希凱,2

1. 中國科學(xué)院電子學(xué)研究所微波成像技術(shù)國家級(jí)重點(diǎn)實(shí)驗(yàn)室，北京 100190； 2. 中國科學(xué)院大學(xué)，北京 100049

干涉合成孔徑雷達(dá)(interferometric SAR,InSAR)利用交軌向雙天線或重復(fù)軌跡飛行對(duì)地面同一區(qū)域進(jìn)行成像，由于觀測(cè)目標(biāo)到雙天線存在波程差，從而可利用干涉形成的相位圖反演地面高程信息。InSAR技術(shù)將合成孔徑雷達(dá)技術(shù)和干涉測(cè)量相結(jié)合，具有測(cè)繪范圍廣、測(cè)量精度高，能夠全天時(shí)、全天候工作等優(yōu)點(diǎn)，在地形測(cè)繪方面具有良好的應(yīng)用前景。

在實(shí)際干涉SAR系統(tǒng)中，由于系統(tǒng)采樣延時(shí)誤差、通道間相位誤差、基線矢量測(cè)量誤差等因素的存在，使得三維重建幾何模型中的參數(shù)取值不準(zhǔn)確，需要通過干涉參數(shù)定標(biāo)來校準(zhǔn)各參數(shù)，以提高平面定位和高程反演的精度[1-2]。文獻(xiàn)[3—5]給出了單控制點(diǎn)或無控制點(diǎn)SAR圖像定位方法，針對(duì)的是單天線SAR系統(tǒng)定標(biāo)，通過單個(gè)控制點(diǎn)對(duì)雷達(dá)斜距或系統(tǒng)采樣延時(shí)進(jìn)行標(biāo)定，然后根據(jù)SAR成像幾何模型對(duì)圖像中每個(gè)像素點(diǎn)進(jìn)行定位。對(duì)于機(jī)載雙天線InSAR而言，不僅需要對(duì)雷達(dá)斜距進(jìn)行標(biāo)定，更為重要的是基線矢量、相位偏置等干涉參數(shù)測(cè)量誤差引起高程精度的下降。干涉參數(shù)定標(biāo)是InSAR高程測(cè)量精度的保證，是干涉SAR地形測(cè)圖流程的關(guān)鍵環(huán)節(jié)。

干涉參數(shù)定標(biāo)方法從單景定標(biāo)發(fā)展到多景大區(qū)域聯(lián)合定標(biāo)。文獻(xiàn)[6—10]利用平坦地形干涉相位頻率或外部DEM估計(jì)基線參數(shù)，這種方法要求成像場景中存在平坦地形或已知場景地形信息。文獻(xiàn)[11—14]針對(duì)初始相位偏置參數(shù)估計(jì)問題，提出基于外源粗精度DEM或交叉軌數(shù)據(jù)的相位偏置估計(jì)方法；這兩種針對(duì)基線和初始相位偏置兩種干涉參數(shù)單獨(dú)進(jìn)行估計(jì)的方法，均需要外部信息才能實(shí)現(xiàn)對(duì)干涉參數(shù)的估計(jì)，且依賴于外部數(shù)據(jù)的精度，主要研究對(duì)象還是單景干涉圖像對(duì)。文獻(xiàn)[15—16]通過事先布放的大量控制點(diǎn)信息來解算干涉參數(shù)，這是最為直接和精度最高的干涉參數(shù)定標(biāo)方法，但缺點(diǎn)是需要大量的控制點(diǎn)，在實(shí)際工程中難以實(shí)現(xiàn)。針對(duì)大區(qū)域干涉參數(shù)定標(biāo)的難題，國內(nèi)外學(xué)者們相繼提出利用相鄰圖像對(duì)重疊部分連接點(diǎn)信息來降低控制點(diǎn)數(shù)量，文獻(xiàn)[17—21]詳細(xì)闡述了將重疊區(qū)域同名點(diǎn)信息構(gòu)建到干涉定標(biāo)方程中，實(shí)現(xiàn)稀疏控制點(diǎn)InSAR區(qū)域網(wǎng)平差算法。為了解決同名點(diǎn)信息質(zhì)量問題，文獻(xiàn)[22—23]研究了賦予各同名點(diǎn)不同權(quán)值，從而進(jìn)一步提高區(qū)域網(wǎng)平差算法精度。文獻(xiàn)[24—26]在稀疏控制點(diǎn)加上同名點(diǎn)聯(lián)合平差的思想框架里，從InSAR構(gòu)像模型出發(fā)，進(jìn)一步完善了干涉參數(shù)聯(lián)合定標(biāo)方程。

基于稀疏控制點(diǎn)和同名點(diǎn)連接的區(qū)域網(wǎng)聯(lián)合干涉參數(shù)定標(biāo)方法，理論上能夠有效實(shí)現(xiàn)控制點(diǎn)稀疏分布，但是由于干涉參數(shù)定標(biāo)方程的病態(tài)性[19]，要求增加較多GCP(ground control points)來減小干涉參數(shù)間的耦合性引起的高程誤差。然而，與攝影測(cè)量遙感不同的是，為了能夠在雷達(dá)圖像中定位識(shí)別，地面控制點(diǎn)為角反射器，在InSAR作業(yè)前需要外業(yè)人員現(xiàn)場布放，這給外定標(biāo)工作帶來巨大的工作量，特別是在山區(qū)、沙漠等人難以到達(dá)的區(qū)域。為了進(jìn)一步降低控制點(diǎn)數(shù)量，本文提出一種聯(lián)合對(duì)飛數(shù)據(jù)的單控制點(diǎn)干涉SAR定標(biāo)算法。

1 InSAR三維定位建模

干涉SAR三維重建幾何關(guān)系如圖1所示，A1、A2是雙天線相位中心。以主天線相位中心為坐標(biāo)原點(diǎn)建立SAR成像坐標(biāo)系A(chǔ)1-xyz，x軸指向平臺(tái)飛行方向，z軸朝下，y軸與xz構(gòu)成右手系。設(shè)平臺(tái)高度為H，成像區(qū)域中任意一目標(biāo)點(diǎn)P高程為h，主天線A1相對(duì)目標(biāo)P點(diǎn)的雷達(dá)斜距為r，副天線A2對(duì)應(yīng)的雷達(dá)斜距為r+Δr。副天線相位中心A2的坐標(biāo)為(Bx,By,Bz)，Bx是由于SAR斜視觀測(cè)(姿態(tài)擾動(dòng))引起的順軌基線，由于干涉高程測(cè)量對(duì)順軌基線不敏感[2]，這里直接忽略，令Bx=0，By和Bz是基線矢量在交軌平面上投影對(duì)應(yīng)的水平基線和垂直基線，即

By=Bcosα

(1)

Bz=Bsinα

(2)

圖1 InSAR三維定位基本原理Fig.1 Principle of InSAR three-dimensional location

式中，基線長度B=|A1A2|，α是A1A2矢量與水平面夾角，即基線傾角。干涉合成孔徑雷達(dá)通過視角差異引起的干涉相位來反演目標(biāo)高程信息，由圖1中三維重建幾何關(guān)系，目標(biāo)高程可表示為

h=H-rcosθ

(3)

可見雷達(dá)主天線相對(duì)目標(biāo)P下視角θ包含了目標(biāo)高程信息，而雙天線接收到信號(hào)的干涉相位差φ能夠真實(shí)反映入射角的變化

(4)

式中，λ是SAR發(fā)射信號(hào)波長，干涉相位差可通過雙天線得到的復(fù)圖像對(duì)共軛相乘來計(jì)算。式(3)、式(4)為InSAR高程測(cè)量的基本公式，表明高程測(cè)量與干涉相位、基線長度、基線傾角、平臺(tái)高度和SAR斜距5個(gè)特征參數(shù)測(cè)量有關(guān)，這5個(gè)參數(shù)的測(cè)量誤差均會(huì)引起高程測(cè)量性能的退化。

再來考慮目標(biāo)P點(diǎn)圖像平面幾何位置。如圖1所示，InSAR在對(duì)地成像時(shí)，通過距離方程、多普勒方程和干涉相位方程描述InSAR基本測(cè)量值與地面目標(biāo)位置之間的關(guān)系。干涉相位信息包含了第三維高程信息，平面位置信息可通過距離多普勒成像方程來計(jì)算

r=|A1-P|

(5)

(6)

式中，fd=2|v|sinθsq/λ是方位多普勒頻率；v是速度矢量；A1(x0,y0,H)是主天線相位中心空間三維位置；P(x′,y′,h)是待解算目標(biāo)點(diǎn)三維空間位置。首先來定義地理空間直角坐標(biāo)系x′y′z′，x′軸、y′軸分別指向地理東向和北向，z′軸和InSAR成像坐標(biāo)系z(mì)軸重合。將成像坐標(biāo)系xyz繞z軸逆時(shí)針旋轉(zhuǎn)角度?(InSAR成像參考航跡地理方位角)，并將坐標(biāo)原點(diǎn)移至地理坐標(biāo)系原點(diǎn)，就可得到空間直角坐標(biāo)系x′y′z′。根據(jù)圖1空間成像幾何及式(5)、式(6)，化簡可得圖像空間P點(diǎn)在地理空間直角坐標(biāo)系中坐標(biāo)

(7)

式中，γ是SAR成像斜視角在水平面上的投影，可由斜視角和雷達(dá)斜距計(jì)算得到

(8)

分析式(7)、式(8)得知，InSAR平面定位與SAR位置、雷達(dá)斜距、多普勒頻率、目標(biāo)高程和地理方位角有關(guān)。而SAR位置可由高精度定位定姿系統(tǒng)事先測(cè)量已知。多普勒頻率和地理方位角雖然也可由POS系統(tǒng)提供的姿態(tài)信息進(jìn)行計(jì)算，由式(7)可知，平面定位誤差對(duì)這兩個(gè)角度敏感度與斜距成正比，因此需要對(duì)二角度進(jìn)行精確的估計(jì)。

至此，式(3)—式(8)描述了InSAR測(cè)量幾何中目標(biāo)點(diǎn)P(x′,y′,h)空間三維定位模型。總的說來，影響三維定位精度的因素包括平臺(tái)高度、雷達(dá)斜距r、基線長度B、基線傾角α、干涉相位φ、多普勒頻率和參考航跡地理方位角?。平臺(tái)高度參數(shù)的測(cè)量精度由高精度POS測(cè)量來保障；雷達(dá)斜距由系統(tǒng)初始采樣延遲(電纜延遲、大氣延遲等)誤差引起；基線參數(shù)、多普勒頻率和參考航跡地理方位角均受平臺(tái)姿態(tài)變換影響隨時(shí)間發(fā)生變化。

2 單控制點(diǎn)InSAR參數(shù)定標(biāo)算法

2.1 三維定標(biāo)方程降維

為了降低三維定標(biāo)方程的自由度，對(duì)干涉測(cè)量的特征參數(shù)的獨(dú)立性和相關(guān)性進(jìn)行系統(tǒng)的分析。SAR圖像各像素點(diǎn)三維定位是通過雷達(dá)天線相位中心絕對(duì)定位基準(zhǔn)通過成像幾何傳遞來實(shí)現(xiàn)，因此平臺(tái)高度H完全由POS高度向測(cè)量精度決定。雷達(dá)斜距誤差主要由系統(tǒng)初始采樣延遲(電纜延遲、大氣延遲等)誤差引起，初始采樣延遲誤差在環(huán)境穩(wěn)定雷達(dá)系統(tǒng)安裝方式固定情況下是穩(wěn)定的，因此可以利用1個(gè)控制點(diǎn)三維空間位置計(jì)算其與平臺(tái)之間的實(shí)際距離來校正

(9)

式中，(xgcp,ygcp,zgcp)是地面控制點(diǎn)三維位置；(xr,yr,zr)是雷達(dá)主天線相位中心三維位置；ri是地面控制點(diǎn)在SAR圖像中對(duì)應(yīng)的實(shí)際斜距值。

由式(7)、式(8)可知，多普勒頻率的估計(jì)精度可等效為斜視角在水平面上的投影角γ，而γ與參考航跡地理方位角?的差進(jìn)一步對(duì)平面定位精度產(chǎn)生影響，因此可將兩者之差γ-?作為一個(gè)特征參數(shù)同時(shí)解算，從而可降低控制點(diǎn)數(shù)量。因此，在平面定位方程(7)和式(8)中，首先利用一個(gè)控制點(diǎn)斜距誤差校正之后，平面位置只與高程和γ-?有關(guān)

(x′,y′)=f(h,γ-?)

(10)

高度誤差由POS系統(tǒng)保障，利用一個(gè)控制點(diǎn)的空間三維位置信息(x′,y′,h)可分別對(duì)雷達(dá)斜距和參考航跡地理方位角進(jìn)行獨(dú)立的校正。然而，剩下的3個(gè)特征參數(shù)(基線長度、基線傾角和初始相位偏置)存在相互耦合，難以獨(dú)立地分離出來。由于雷達(dá)天線相位中心是一種等效的邏輯中心，不在天線的幾何中心，基線長度不能在地面直接測(cè)量獲得；基線傾角測(cè)量雖然可以通過POS姿態(tài)測(cè)量獲得，但POS測(cè)量坐標(biāo)系與成像坐標(biāo)系存在固定的安裝誤差，因此也需要通過外定標(biāo)的方式校準(zhǔn)；干涉相位是通過復(fù)圖像對(duì)經(jīng)過配準(zhǔn)、濾波和相位展開后獲得，但由于三角函數(shù)的纏繞特性，經(jīng)過相位解纏后的干涉相位φunw不是絕對(duì)相位，仍然存在一個(gè)初始相位偏置φ0

φ=φ0+φunw

(11)

干涉相位測(cè)量誤差除了圖像信噪比引起的隨機(jī)相位誤差之外，還存在雷達(dá)系統(tǒng)初始相位漂移引起的固定誤差，因此需要通過外定標(biāo)的手段進(jìn)行校準(zhǔn)。由高程反演式(3)、式(4)可知，初始相位偏置估計(jì)與基線長度和基線傾角相耦合，難以獨(dú)立的分離出來。因此，7個(gè)自由度的InSAR干涉參數(shù)定標(biāo)方程可以降低到3個(gè)未知參數(shù)，需要地面至少有3個(gè)控制點(diǎn)才能解算出耦合的3個(gè)特征參數(shù)。

為了更進(jìn)一步降低定標(biāo)方程的維度，減少控制點(diǎn)數(shù)量，本文提出增加一條對(duì)飛航線對(duì)目標(biāo)區(qū)域重復(fù)觀測(cè)，如圖2所示。假設(shè)第1軌觀測(cè)稱為正飛，從相反方向?qū)ν粔K區(qū)域觀測(cè)稱為反飛。主要思想是兩次觀測(cè)擁有相同目標(biāo)點(diǎn)(同名點(diǎn))，就可以依賴圖像中更多的分布式目標(biāo)點(diǎn)建立大量的定標(biāo)方程，利用增加的反飛航線獲得額外的觀測(cè)信息，從而降低對(duì)控制點(diǎn)數(shù)量的要求。

圖2 InSAR對(duì)飛觀測(cè)Fig.2 InSAR observing in two converse flights

2.2 對(duì)飛數(shù)據(jù)同名點(diǎn)定標(biāo)

由2.1節(jié)三維定標(biāo)方程中，平臺(tái)高度、雷達(dá)斜距、多普勒頻率和參考航跡地理方位角4個(gè)參數(shù)的測(cè)量誤差可以成功從InSAR三維定標(biāo)方程中分離出來，利用一個(gè)控制點(diǎn)即可完成這4個(gè)參數(shù)測(cè)量誤差的校正。根據(jù)圖1干涉SAR高程反演基本原理，目標(biāo)高程與基線參數(shù)和干涉相位的函數(shù)關(guān)系可寫成

(12)

式(12)中包含基線長度、基線傾角和初始相位偏置3個(gè)未知數(shù)，理論上至少需要3個(gè)控制點(diǎn)，這是前面提到的每個(gè)條帶至少包含3個(gè)控制點(diǎn)才能完成干涉參數(shù)定標(biāo)解算的理論來源。本文提出通過對(duì)同一成像區(qū)域在反方向增加一條反飛航線，利用增加的對(duì)飛與正飛干涉圖像對(duì)間的內(nèi)在聯(lián)系，建立新的參數(shù)定標(biāo)方程，從而減少控制點(diǎn)數(shù)量。

2.2.1 基于SIFT的同名點(diǎn)提取流程

首先需要通過SAR圖像匹配建立起正反飛干涉圖像對(duì)間聯(lián)系，必須對(duì)兩組干涉圖像對(duì)進(jìn)行匹配處理。由于兩組數(shù)據(jù)為不同方向和角度獲取的同一區(qū)域影像，且存在地物類型差異及SAR側(cè)視成像幾何，不能實(shí)現(xiàn)圖像中每個(gè)像素點(diǎn)完全匹配。本文考慮采用尺度不變特征變換(scale-invariant feature transform，SIFT)算法對(duì)正反兩幅數(shù)據(jù)中提取一定數(shù)量的同名點(diǎn)，SIFT算法提取的特征對(duì)圖像尺度變化、旋轉(zhuǎn)、仿射變換、噪聲污染、明亮度變化及視角變化具有良好的不變性。SIFT方法是一種提取圖像局部特征的有效算法，它能夠在尺度空間內(nèi)尋找到一些極值點(diǎn)，對(duì)圖像的亮度、平移、旋轉(zhuǎn)、尺度變化具有較強(qiáng)的適應(yīng)性，利用特征點(diǎn)周圍圖像提取該特征點(diǎn)的特征描述符，從而可以在特征描述符之間進(jìn)行匹配。因此雖然二次SAR圖像獲取存在較大的角度差異，但是SIFT算法目前能夠很好地適應(yīng)SAR圖像同名點(diǎn)提取應(yīng)用[27-28]。同名點(diǎn)提取算法處理流程如圖3所示，主要包含構(gòu)建尺度空間、計(jì)算主方向和描述符、計(jì)算最小歐氏距離、隨機(jī)抽樣一致性算法(RANSAC)剔除誤匹配點(diǎn)、基于質(zhì)量圖的低相干區(qū)同名點(diǎn)剔除。同名點(diǎn)提取完成后得到兩幅圖像中同一目標(biāo)點(diǎn)的圖像空間坐標(biāo)信息。

2.2.2 同名點(diǎn)定標(biāo)方程

在完成同名點(diǎn)信息提取后，然后針對(duì)兩組干涉復(fù)圖像對(duì)分別經(jīng)過配準(zhǔn)、濾波和相位解纏處理，得到兩組干涉數(shù)據(jù)解纏后的干涉相位矩陣φ1、φ2。根據(jù)提取的同名點(diǎn)在各自圖像空間中坐標(biāo)信息，在解纏后的干涉相位中提取相應(yīng)的干涉相位值，建立同名點(diǎn)約束方程

(13)

式中，k=1,2,…,N，且N是同名點(diǎn)數(shù)量；φ1k、φ2k是第k個(gè)同名點(diǎn)在兩組干涉數(shù)據(jù)中的解纏后相位；H1、H2是正飛和反飛航跡下平臺(tái)高度；r1k、r2k分別是第k個(gè)同名點(diǎn)對(duì)應(yīng)的斜距；B1、α1、φ1、B2、α2、φ2分別是正反飛條帶基線長度、基線傾角和初始相位偏置值，是待求解的6個(gè)未知參數(shù)。分析式(13)可知，理論上只要找到6個(gè)二次干涉圖像對(duì)間的同名點(diǎn)，不需要任何無控制點(diǎn)的情況下可解算出6個(gè)未知的干涉定標(biāo)參數(shù)。然而，干涉參數(shù)的反演精度受入射角、基線參數(shù)、干涉相位噪聲等因素的影響。定性的分析可知，兩組方程特征參數(shù)差異越大方程越穩(wěn)健，而實(shí)際中平臺(tái)高度、雷達(dá)斜距、基線長度和基線傾角4個(gè)參數(shù)受平臺(tái)飛行和設(shè)備安裝條件的限制往往難以做到較大的差異，因此只剩下可以通過干涉相位差異來提升定標(biāo)方程的穩(wěn)健性。

圖3 基于SIFT算法的SAR圖像同名點(diǎn)提取流程Fig.3 Tie-points extraction flowchart based on SIFT

干涉相位值的差異與地形起伏程度和入射角差異關(guān)系密切。本文設(shè)計(jì)通過對(duì)向飛行實(shí)現(xiàn)入射角差異最大化，但是入射角受雷達(dá)波束寬度的限制不可能無限大，因此只剩下地形起伏程度來增大兩次觀測(cè)間的干涉相位差異。為了定量化的解釋此問題，引入定標(biāo)矩陣條件數(shù)的概念來定量化地說明地形起伏對(duì)于聯(lián)合對(duì)飛數(shù)據(jù)定標(biāo)方程穩(wěn)健性問題。矩陣條件數(shù)等于F的范數(shù)與F逆矩陣范數(shù)的乘積。矩陣條件數(shù)決定了對(duì)應(yīng)的線性方程組的數(shù)值解的好壞和可信程度。條件數(shù)越小，則考慮的線性方程組的數(shù)值解的精度越高。高程反演敏感度矩陣F的條件數(shù)可以表示為

(14)

仿真了X波段2 m基線機(jī)載雙天線干涉SAR情況下，在高程起伏標(biāo)準(zhǔn)差分別在0 m、20 m、40 m條件下聯(lián)合對(duì)飛數(shù)據(jù)的定標(biāo)方程矩陣條件數(shù)與完全利用控制點(diǎn)定標(biāo)情況下的矩陣條件數(shù)對(duì)比情況。由圖4可知，隨著地形起伏程度和入射角差異的增大，矩陣條件數(shù)變小，方程更為穩(wěn)健，這與前面定性分析的結(jié)論是一致的。同時(shí)在地形起伏標(biāo)準(zhǔn)差達(dá)到40 m時(shí)，矩陣條件數(shù)已逐漸逼近到完全依賴控制點(diǎn)的干涉參數(shù)定標(biāo)矩陣條件數(shù)變化曲線。

圖4 定標(biāo)方程矩陣條件數(shù)隨地形起伏變化Fig.4 Calibration matrix condition number under different terrain conditions

地形起伏越大，干涉定標(biāo)方程越穩(wěn)健。然而地形起伏越大，地表高程對(duì)基線參數(shù)估計(jì)影響越大，當(dāng)用于定標(biāo)的同名點(diǎn)存在匹配誤差時(shí)，將會(huì)引起干涉參數(shù)的反演精度下降。因此，在地形起伏較大區(qū)域?qū)νc(diǎn)匹配精度提出了更高的要求。實(shí)際中同名點(diǎn)選取應(yīng)綜合考慮平地和山區(qū)，使其兼顧穩(wěn)健性和精度的要求。在實(shí)際InSAR地形測(cè)圖應(yīng)用中，無法保證所有的區(qū)域均有較大的地形起伏，因而筆者提出單控制點(diǎn)干涉參數(shù)定標(biāo)方法，利用控制點(diǎn)絕對(duì)位置信息，降低算法對(duì)地形起伏的要求，從而提高算法地形適應(yīng)性能。

2.3 聯(lián)合對(duì)飛數(shù)據(jù)的單控制點(diǎn)參數(shù)定標(biāo)

完全依賴同名點(diǎn)定標(biāo)方程解算干涉參數(shù)在地形平坦區(qū)域性能無法保障，因此引入一個(gè)單控制點(diǎn)定標(biāo)方程，來提高算法的適應(yīng)性。假設(shè)控制點(diǎn)對(duì)應(yīng)的高程為hgcp，在正反兩組解纏后的干涉相位圖中相位值為φ1gcp和φ2gcp，該單控制點(diǎn)可在兩組干涉數(shù)據(jù)中分別建立一個(gè)定標(biāo)方程，如式(15)所示

(15)

式中，f(·)是干涉SAR高程反演模型，可參考式(12)。控制點(diǎn)具有外部的絕對(duì)的高程信息，因此每個(gè)方程解算是相對(duì)穩(wěn)健的。根據(jù)2.2.2節(jié)，還可以構(gòu)建N組同名點(diǎn)定標(biāo)方程

(16)

同名點(diǎn)約束方程中，由于沒有絕對(duì)的高程信息參考，但是兩組干涉參數(shù)間仍然存在相互耦合，特別是地形平坦的區(qū)域。雖然正反對(duì)飛入射角的差異能夠一定程度上緩解這個(gè)問題，但是仍然存在一定的誤差耦合?；诖耍瑢⑹?15)和式(16)聯(lián)立起來，方程穩(wěn)健性將得到很大的提升。因此，干涉參數(shù)定標(biāo)的問題就歸結(jié)為N+2個(gè)方程組(聯(lián)立式(15)、式(16))解6個(gè)未知數(shù)的數(shù)學(xué)問題。由于該方程組是超定的非線性方程組無解析解，擬采用最優(yōu)化函數(shù)的方法進(jìn)行求解。最優(yōu)化函數(shù)能夠獲取全局最優(yōu)解，避免方程組取得局部最優(yōu)解，從而降低干涉參數(shù)的定標(biāo)精度。聯(lián)立控制點(diǎn)定標(biāo)方程式(15)和同名點(diǎn)定標(biāo)式(16)，構(gòu)建如下的最優(yōu)化函數(shù)

(17)

式中,θ1gcp、θ2gcp是該單個(gè)控制點(diǎn)相對(duì)雷達(dá)入射角，可通過控制點(diǎn)高程值計(jì)算。A0、Ak分別是控制點(diǎn)方程和同名點(diǎn)方程不同的權(quán)值，受InSAR相位噪聲和同名點(diǎn)匹配精度影響，一般A0>Ak。使得最優(yōu)化函數(shù)M達(dá)到最小值的干涉參數(shù)(B1,α1,φ1,B2,α2,φ2)即為方程組的解，為最終的干涉參數(shù)定標(biāo)結(jié)果。由于同名點(diǎn)為非人工布設(shè)的定標(biāo)點(diǎn)，其信噪比直接決定了干涉參數(shù)定標(biāo)精度。評(píng)價(jià)同名點(diǎn)相位噪聲質(zhì)量通常采用干涉相關(guān)系數(shù)來表征，因此針對(duì)每個(gè)同名點(diǎn)對(duì)最優(yōu)化函數(shù)的貢獻(xiàn)進(jìn)行了加權(quán)處理，加權(quán)系數(shù)與該點(diǎn)的相干系數(shù)值相關(guān)

(18)

式中，S1和S2是干涉圖像對(duì)的單視復(fù)圖像。

式(17)是解定標(biāo)方程構(gòu)建的最優(yōu)化函數(shù)，即找到一組干涉定標(biāo)參數(shù)，使得在全局范圍內(nèi)最優(yōu)化函數(shù)達(dá)到最小值。在求解此最優(yōu)化問題時(shí)，最優(yōu)化理論中包括最速下降法、牛頓法、共軛梯度法、遺傳算法、模擬退火算法等。采用經(jīng)典的最速下降法進(jìn)行求解定標(biāo)方程，求解步驟如下：

為了驗(yàn)證添加一個(gè)控制點(diǎn)能夠提高定標(biāo)方程組的穩(wěn)健性，這里給出了聯(lián)立單控制點(diǎn)和對(duì)飛定標(biāo)方程的矩陣條件數(shù)仿真結(jié)果，如圖5所示，在使用了地形起伏程度在0 m情況下，單控制點(diǎn)定標(biāo)在入射角差異為30°時(shí)已經(jīng)逼近于全控制點(diǎn)定標(biāo)的結(jié)果，可見引入單個(gè)控制點(diǎn)可以大大降低本文方法對(duì)地形起伏程度的依賴。

圖5 單控制點(diǎn)定標(biāo)敏感度矩陣條件數(shù)Fig.5 Calibration matrix condition number using single corner reflector

3 機(jī)載實(shí)際數(shù)據(jù)試驗(yàn)驗(yàn)證

采用中國科學(xué)院電子學(xué)研究所研制的機(jī)載干涉合成孔徑雷達(dá)系統(tǒng)在某機(jī)場開展的飛行試驗(yàn)來進(jìn)行實(shí)際數(shù)據(jù)驗(yàn)證。該InSAR系統(tǒng)工作在X波段，具備全極化干涉數(shù)據(jù)獲取能力，最高分辨率0.5 m，搭載獎(jiǎng)狀Ⅱ型機(jī)載平臺(tái)，如圖6所示。雙天線水平安裝在機(jī)腹下小型吊艙中，物理基線長度達(dá)2.2 m，系統(tǒng)設(shè)計(jì)高程測(cè)量精度滿足1∶10 000比例尺地形測(cè)圖精度要求。

圖6 機(jī)載InSAR系統(tǒng)及飛行平臺(tái)Fig.6 Airborne InSAR system and flying platform

采用的一組對(duì)飛干涉數(shù)據(jù)獲取于2013年，對(duì)飛兩組數(shù)據(jù)為一個(gè)架次飛行試驗(yàn)中獲取，兩條設(shè)計(jì)航線觀測(cè)同一塊成像區(qū)域，平臺(tái)飛行相對(duì)地面高度3000 m，雷達(dá)中心下視角45°，測(cè)繪幅寬3 km。圖7是對(duì)飛兩次觀測(cè)獲取的主天線SAR斜距影像。平臺(tái)東往西飛行時(shí)，雙天線SAR工作在乒乓收發(fā)模式，有效基線約2.2 m；西往東飛時(shí)，雙天線SAR工作在一發(fā)雙收模式，有效基線減半。成像區(qū)域中布放了5個(gè)地面控制點(diǎn)，GCP1—GCP5分布如圖7(a)所示?？刂泣c(diǎn)布放范圍主要集中在波束中心附近區(qū)域，沒有完全覆蓋整個(gè)測(cè)繪帶，因而這里僅將控制點(diǎn)覆蓋區(qū)域(圖7中虛線方框)作為精度驗(yàn)證的有效區(qū)域。該成像區(qū)域地形較為平坦，高程落差絕大部分小于5 m。針對(duì)平面定位定標(biāo)和干涉參數(shù)定標(biāo)實(shí)驗(yàn)結(jié)果分別進(jìn)行分析。其中干涉參數(shù)定標(biāo)分別針對(duì)多控制點(diǎn)、對(duì)飛無控制點(diǎn)和對(duì)飛單控制點(diǎn)3種定標(biāo)方案的定標(biāo)結(jié)果進(jìn)行比較。

圖7 InSAR對(duì)飛獲取的兩景影像Fig.7 Two images of InSAR in two converse flights

首先分析平面位置定標(biāo)精度(如表1所示)，由于高度測(cè)量誤差直接由POS系統(tǒng)決定，一般其測(cè)量精度優(yōu)于0.1 m，誤差量較小。雷達(dá)初始斜距誤差受通道信號(hào)傳輸延遲、大氣傳輸?shù)纫蛩氐挠绊?，在機(jī)載大氣環(huán)境相對(duì)穩(wěn)定情況下，主要由通道信號(hào)傳輸過程中引起的時(shí)間延遲決定，因此在安裝方式固定情況下是相對(duì)穩(wěn)定值。結(jié)合POS數(shù)據(jù)及地面控制點(diǎn)信息，計(jì)算得到斜距誤差約為42.1 m。下面來看γ-?參數(shù)定標(biāo)情況，正飛角度誤差估計(jì)值0.015°，引起的平面定位誤差如圖8所示，東西向引起最大1 m的平面定位誤差，而南北向由于位于距離向，因此對(duì)方位角的估計(jì)誤差不是很敏感?？梢?，微小的角度差γ-?估計(jì)誤差會(huì)引起較大的平面定位誤差，而傳統(tǒng)基于InSAR成像模型干涉定標(biāo)方法[17-21]只考慮了干涉高程反演參數(shù)的校正，忽略了多普勒中心頻率和參考航跡地理方位角參數(shù)的校正，因此必然導(dǎo)致平面定位精度的下降。

表1 平面位置參數(shù)定標(biāo)結(jié)果

圖8 傳統(tǒng)方法平面位置未校正誤差Fig.8 Planar location error using Conditional method

再對(duì)高程測(cè)量特征參數(shù)反演精度進(jìn)行分析。這里針對(duì)多控制點(diǎn)定標(biāo)、對(duì)飛無控制點(diǎn)定標(biāo)和對(duì)飛單控制點(diǎn)定標(biāo)3種方案進(jìn)行了比較。多控制點(diǎn)定標(biāo)表示利用場景中的5個(gè)控制點(diǎn)對(duì)正飛InSAR數(shù)據(jù)進(jìn)行干涉參數(shù)定標(biāo)；對(duì)飛無控制點(diǎn)定標(biāo)指的是完全利用對(duì)飛數(shù)據(jù)中同名點(diǎn)定標(biāo)方程對(duì)干涉參數(shù)進(jìn)行反演，在地形平坦的試驗(yàn)區(qū)理論性能將不高；對(duì)飛單控制點(diǎn)定標(biāo)是本文提出的干涉參數(shù)定標(biāo)算法，相比于對(duì)飛無控制點(diǎn)定標(biāo)方案能夠適應(yīng)平坦區(qū)域的干涉參數(shù)定標(biāo)解算。首先對(duì)兩組對(duì)飛的干涉圖像進(jìn)行同名點(diǎn)提取，同名點(diǎn)提取結(jié)果分布如圖9所示，同名點(diǎn)主要集中在波束中心區(qū)域，完全覆蓋了控制點(diǎn)區(qū)域，且在距離向形成一定的分布。經(jīng)過RANSAC剔除誤匹配點(diǎn)后共產(chǎn)生335組同名點(diǎn)，再根據(jù)干涉圖像對(duì)的質(zhì)量圖信息剔除相干性差的點(diǎn)(相位誤差大)，仍然能夠提取到大量的同名點(diǎn)信息(82組)?？梢娂词故莾纱潍@取的圖像不同方向獲取，由于基于SIFT的同名點(diǎn)提取具備尺度、旋轉(zhuǎn)不變性，提取正確度和精度都較高。

圖9 同名點(diǎn)提取結(jié)果Fig.9 Tie-points extraction results

根據(jù)提取的同名點(diǎn)信息，聯(lián)合控制點(diǎn)和同名點(diǎn)聯(lián)合定標(biāo)方程，利用構(gòu)建的最優(yōu)化函數(shù)(式(17))求解干涉參數(shù)。3種定標(biāo)方案輸出的高程誤差分析情況如圖10所示。圖10(a)中傳統(tǒng)的多控制點(diǎn)定標(biāo)結(jié)果與本文提出的對(duì)飛單控制點(diǎn)定標(biāo)算法的高程反演值吻合較好，而對(duì)飛無控制點(diǎn)高程反演曲線在高度向存在一個(gè)較大的常數(shù)誤差，這也證明了在地形平坦區(qū)域完全依賴同名點(diǎn)間聯(lián)系，是難以獲得高精度的干涉參數(shù)定標(biāo)結(jié)果。以傳統(tǒng)的多控制點(diǎn)定標(biāo)結(jié)果作為真值，分析另外2種方案的定標(biāo)誤差，如圖10(b)所示。對(duì)飛單控制點(diǎn)方法趨近于0，已經(jīng)很好地逼近多控制點(diǎn)的定標(biāo)結(jié)果。以正飛航線干涉測(cè)量為例，干涉參數(shù)定標(biāo)結(jié)果及高程誤差統(tǒng)計(jì)情況如表2。本文提出的定標(biāo)方案與傳統(tǒng)多控制點(diǎn)干涉參數(shù)反演結(jié)果有較大的差異，基線長度估計(jì)值差異2 mm，基線傾角估計(jì)值差異0.05°，而初始相位偏置更在1 rad附近。這主要是由于3個(gè)干涉參數(shù)誤差是相互耦合的，干涉參數(shù)反演的結(jié)果不是其真實(shí)值，而是在高程誤差最小情況下的最優(yōu)值。從表2中高程誤差的統(tǒng)計(jì)結(jié)果可看出，平坦場景情況下對(duì)飛無控制點(diǎn)定標(biāo)方案存在高程誤差均值為3.39 m，高程誤差標(biāo)準(zhǔn)差為0.14 m，無法滿足該InSAR系統(tǒng)指標(biāo)設(shè)計(jì)要求。由于無絕對(duì)參考控制點(diǎn)條件下對(duì)飛定標(biāo)方法，無同名連接點(diǎn)高程真值，即使定標(biāo)方程兩側(cè)出現(xiàn)同樣常數(shù)誤差，定標(biāo)方程仍然成立，這是對(duì)飛無控制點(diǎn)定標(biāo)方法產(chǎn)生較大的常數(shù)誤差主要原因。因此，本文提出單控制點(diǎn)定標(biāo)，該定標(biāo)方案高程誤差均值和標(biāo)準(zhǔn)差均在0.05 m附近，比該機(jī)載InSAR系統(tǒng)精度設(shè)計(jì)指標(biāo)低1個(gè)量級(jí)，因此是能夠滿足高精度的InSAR地形測(cè)繪應(yīng)用需求的。

圖10 3種定標(biāo)方案高程解算值比較Fig.10 DEM comparison of three calibration methods

誤差正飛多控制點(diǎn)對(duì)飛無控制點(diǎn)對(duì)飛單控制點(diǎn)基線長度/m2.19242.19852.1944基線傾角/(°)2.9612.9952.911初始相位偏置/rad695.134697.135696.229高程誤差均值/m03.390.053高程誤差標(biāo)準(zhǔn)差/m00.140.050

經(jīng)過平面定位、高程反演和地理編碼處理后，圖11給出了單控制點(diǎn)定標(biāo)后經(jīng)過處理獲取的場景區(qū)域數(shù)字正射影像圖(DOM)和數(shù)字高程模型(DEM)圖。更進(jìn)一步的，利用沒有參與定標(biāo)的剩余4個(gè)控制點(diǎn)對(duì)高程精度進(jìn)行檢查，高程誤差值分別為-0.048 m、-0.01 m、0.047 m和0.052 m。綜上所述，本文提出的方法與傳統(tǒng)多控制點(diǎn)的定標(biāo)結(jié)果差異小于0.05 m，比標(biāo)稱的InSAR系統(tǒng)干涉高程測(cè)量精度(0.5 m)低一個(gè)數(shù)量級(jí)，能夠滿足該機(jī)載InSAR系統(tǒng)地形測(cè)圖精度要求。

圖11 單控制點(diǎn)DOM和DEM反演結(jié)果Fig.11 DOM and DEM results using single corner reflector

4 結(jié) 論

本文提出了一種聯(lián)合對(duì)飛數(shù)據(jù)的單控制點(diǎn)干涉SAR參數(shù)定標(biāo)算法，在構(gòu)建InSAR平面和高程三維定標(biāo)模型的基礎(chǔ)上，提出了單個(gè)控制點(diǎn)加對(duì)飛數(shù)據(jù)同名點(diǎn)約束的最優(yōu)化聯(lián)合定標(biāo)方程。以某機(jī)場區(qū)域X波段2.2 m基線機(jī)載InSAR獲取的實(shí)際數(shù)據(jù)進(jìn)行驗(yàn)證，結(jié)果表明，本文提出的定標(biāo)算法與傳統(tǒng)多控制點(diǎn)干涉定標(biāo)高程差異小于0.05 m，比該機(jī)載InSAR測(cè)圖系統(tǒng)設(shè)計(jì)精度低一個(gè)量級(jí)，較好地逼近傳統(tǒng)的多控制點(diǎn)干涉定標(biāo)精度。