地基干涉合成孔徑雷達(dá)圖像非線性大氣相位補償方法

胡 程 鄧云開 田衛(wèi)明 曾 濤

①(北京理工大學(xué)信息與電子學(xué)院雷達(dá)技術(shù)研究所 北京 100081)

②(北京理工大學(xué)衛(wèi)星導(dǎo)航電子信息技術(shù)教育部重點實驗室 北京 100081)

1 引言

作為一種高精度的形變測量儀器，地基干涉合成孔徑雷達(dá)(Ground-Based Interferometric Synthetic Aperture Radar,GB-InSAR)已經(jīng)在形變監(jiān)測領(lǐng)域得到了廣泛的應(yīng)用[1]。GB-InSAR通常是基于差分干涉測量技術(shù)，通過對同一位置、不同時刻獲取的兩幅SAR圖像進行差分干涉處理，基于相位信息來實現(xiàn)形變測量，其一般工作在X或者Ku波段，形變測量精度可以達(dá)到毫米或者亞毫米量級。GB-InSAR測量誤差的主要來源是大氣相位，由于不同時刻氣象條件的不同，電磁波在大氣中傳播的速度會發(fā)生改變，從而導(dǎo)致大氣延時誤差[2]。

為實現(xiàn)大氣相位的補償，國內(nèi)外學(xué)者主要提出了3種解決方法。第1種的補償方法是在觀測場景內(nèi)建立氣象站，基于大氣折射率模型，利用氣象數(shù)據(jù)(溫度、濕度、大氣壓)來對大氣相位進行定量估計[3]；第2種是在場景中人共布設(shè)或者選擇出一些高度穩(wěn)定的特征點，采用空間插值的方法實現(xiàn)對整幅圖像的大氣相位的補償[4,5]；第3種則是基于永久散射體(Permanent Scatterer,PS)技術(shù)，根據(jù)大氣相位的空間分布特征，建立描述大氣相位的方程，估計大氣相位參數(shù)，實現(xiàn)大氣相位的補償[6]。

在基于PS技術(shù)進行大氣相位補償時，首先需要選擇出未處在形變區(qū)域的PS點，然后對大氣相位進行合理的建模，并建立線性方程組，實現(xiàn)對大氣參數(shù)的粗估計，此后剔除與模型偏差較大的PS點，逐步迭代實現(xiàn)大氣參數(shù)的準(zhǔn)確估計。該方法不需要氣象參數(shù)及布設(shè)特征點，可以基于大量的PS點進行大氣參數(shù)的估計，估計精度較高。一般情況下，大氣在空間上均勻變化，可以將大氣相位建模為隨斜距線性變化的分量[7]。在地形陡峭的山區(qū)，大氣在空間上非均勻變化，可以采用兼顧斜距和高程的多參數(shù)模型[8]。

基于PS技術(shù)的大氣相位補償方法，已經(jīng)在GB-InSAR領(lǐng)域得到了廣泛的應(yīng)用，但還存在一些典型的問題。首先，該方法要求采用場景中的非形變PS點建立觀測方程組；其次，氣象條件一直在隨時間改變，大氣相位的時變性很強，在較差天氣條件(降雨、降雪、強風(fēng)等)下，大氣在空間上非均勻變化，導(dǎo)致大氣相位可能表現(xiàn)出復(fù)雜的空變性，無法建立合理的多參數(shù)模型來模擬大氣相位。因此，常規(guī)的基于PS技術(shù)的補償方法，在應(yīng)用于時序GB-InSAR圖像處理時，還存在較大的改進空間。

針對上述問題，本文借鑒在場景中布設(shè)若干個穩(wěn)定的地面控制點，通過空間維插值估計大氣相位的思想。首先采用常規(guī)補償方案對所有的干涉相位圖進行大氣相位補償，并分析PS點的相位序列的標(biāo)準(zhǔn)差，設(shè)定門限選擇出穩(wěn)定PS點；然后基于K均值聚類算法 (K-means clustering algorithm,K-means)將穩(wěn)定PS點劃分出一定數(shù)量的子區(qū)域，將每一個子區(qū)域的中心點設(shè)定為控制點，采用反距離加權(quán)插值算法估計所有PS點的大氣相位，從而實現(xiàn)非線性大氣相位的補償。

本文詳細(xì)介紹了補償方法的實現(xiàn)流程，主要分為穩(wěn)定PS點選擇和空間維插值補償兩部分，然后分別采用改進方法和常規(guī)方法對460幅地基多輸入多輸出 (Multiple-Input Multiple-Output,MIMO)雷達(dá)圖像進行了分析，對比驗證了本文方法的有效性。

2 補償方法

2.1 穩(wěn)定PS點選擇

在利用像素點的相位信息進行形變測量時，差分干涉相位的質(zhì)量直接影響到形變測量的精度。因此，GB-InSAR差分干涉處理時，通常需要選擇出一些高質(zhì)量的像素點，即PS點，來進行形變分析。在GB-InSAR領(lǐng)域，廣泛采用幅度離差法來進行PS點的選擇，一個像素點的時序幅度序列的標(biāo)準(zhǔn)差與均值之比，即為幅度離差指數(shù)。通過設(shè)置合理的幅度離差門限，即可以實現(xiàn)PS點的選擇[9]。

一個PS點的差分干涉相位 Δφ可以建模為

其中，φdefo為形變相位分量；φatm為兩幅圖像獲取期間，由大氣條件改變所導(dǎo)致的大氣相位分量；φnoi為噪聲相位分量；2kπ為相位模糊度，k為整數(shù)。由于相位周期性的影響，Δφ處在-π～π的范圍內(nèi)[10]。在進行大氣相位補償前，要對干涉相位圖進行相位解纏，可以采用非均勻網(wǎng)格下的最小費用流算法，下文中采用 Δφ指代解纏相位。

大氣相位φatm可以建模為

其中，λ為信號的波長，ΔN表示折射率的變化，其隨著空間r和時間t發(fā)生變化，L表示信號的傳輸路徑。一般情況下，大氣的空間同質(zhì)性很好，可以假設(shè)ΔN在空間r上不發(fā)生變化，在時間t上隨機變化，因此可以將φatm建模為隨斜距線性變化的分量

其中，β0為常數(shù)分量，β1表示與斜距相關(guān)的線性系數(shù)，R表示雷達(dá)與目標(biāo)點之間的距離?；赑S技術(shù)進行大氣相位補償時，首先建立線性方程組

ΔΦ為n個PS點的解纏相位Δφ1,Δφ2,···,Δφn構(gòu)成的n×1維向量，X為常數(shù)1與n個PS點的斜距R1,R2,···,Rn構(gòu)成的n×2維矩陣，β為待估計的2×1維向量，ε為n×1維的隨機誤差向量，ε1,ε2,···,εn為各PS點的模型誤差相位分量。采用最小二乘法對β進行估計，可以得到

其中，T表示矩陣的轉(zhuǎn)置。大氣相位的估計分量為

ΔΦ和ΔΦA(chǔ)PS之間的差值即為補償后相位。首先基于所有的PS點對β進行估計，然后為提高估計的精度，剔除一些不可靠PS點。不滿足式(7)準(zhǔn)則的PS點即為不可靠PS點，ΔTatm的取值在0.1～0.2 rad范圍內(nèi)。剔除不可靠PS點后，基于剩余的PS點對大氣相位參數(shù)進行二次估計[11]

除了最基本的線性斜距模型外，常用的參數(shù)模型還包括高階斜距模型、斜距-方位角模型、斜距-高程模型等，如式(8)-式(10)。其中β0,β1和β2是各模型中待估計的未知參數(shù)，R,θ和h分別代表目標(biāo)點的斜距、方位角和高程。對線性方程組式(4)修正后，即可以進行大氣相位補償

基于PS技術(shù)的補償方法，其有效補償?shù)臈l件，首先是要求所采用的參數(shù)模型可以準(zhǔn)確模擬大氣相位，其次是迭代估計參數(shù)時，可以有效剔除不可靠PS點。采用上述方法對一幅干涉相位圖進行大氣相位補償后，一個PS點的補償后相位ΔφAPC可以表示為

2.2 空間維插值補償

2.2.1 子區(qū)域劃分

在選擇出穩(wěn)定PS點后，對于每一幅干涉相位圖，分別進行大氣相位補償。由于穩(wěn)定PS點的相位分量中包含噪聲相位和大氣相位，考慮到噪聲相位在干涉相位圖上隨機變化，不具備空間相關(guān)性，而大氣相位雖然會在整幅圖像范圍內(nèi)發(fā)生變化，但在較小的距離范圍內(nèi)可以視為一個常數(shù)。因此如果對空間上較小距離范圍內(nèi)的所有PS點進行相位平均，則噪聲相位可以得到很好的濾除。由于穩(wěn)定PS點是非均勻的分布在整幅圖像范圍內(nèi)，可以基于K-means算法對穩(wěn)定PS點進行簇劃分，將每一個簇定義為一個子區(qū)域。每一個簇的中心點作為控制點，并對該簇內(nèi)的所有PS點進行相位平均，作為當(dāng)前控制點的相位。

K-means算法的實現(xiàn)原理是對于給定的樣本集x，按照樣本之間的距離大小，將樣本集劃分為K個簇，讓各個簇內(nèi)的點的距離盡可能的小，而讓各個簇之間的點的距離盡可能的大[12]。假設(shè)將樣本集x劃分為K個簇(C1,C2,···,CK)，各簇間的平方誤差和E可以表示為

其中，∥·∥表示2階范數(shù)，即向量的模。μi是簇Ci的均值向量，可以表示為

其中，|·|表示1階范數(shù)，即簇中點的數(shù)量。

經(jīng)過K-means劃分后可以得到K個子區(qū)域，將每一個子區(qū)域中心點定義為控制點CPKM，其相位為該子區(qū)域內(nèi)所有PS點的相位均值。

2.2.2 大氣相位估計

基于這K個控制點CPKM的平均相位，進行空間維插值來估計所有PS點的大氣相位，可以采用反距離加權(quán)插值算法。反距離加權(quán)插值算法是利用已知點與待插點之間的距離來定義加權(quán)因子，然后加權(quán)計算待插點的相位，距離越近加權(quán)比重越大。其計算公式為

其中，Z(x,y)為插值結(jié)果，(x,y)為待插點的空間坐標(biāo)，Zi為第i(i=1,2,···,n)個參考點的數(shù)值，|di|表示待插點與第i個參考點之間的空間距離，μ表示加權(quán)因子的冪指數(shù)，一般取為2[13]。

將這些控制點視為一組離散點，可以基于Delaunay三角剖分法則來構(gòu)建一個三角形網(wǎng)絡(luò)。對于該網(wǎng)絡(luò)內(nèi)的任意一個三角形，其由3個控制點作為頂點，且外接圓中不包含其他控制點。在基于反距離加權(quán)插值算法來估計所有PS點的大氣相位時，如果一個PS點處在某一個三角形內(nèi)部，則以該三角形的3個控制點作為參考點，基于式(15)估計當(dāng)前PS點的大氣相位。如果一個PS點處在所有三角形的外部，則選擇最近的3個控制點，同樣基于式(15)來估計大氣相位。由于每一個控制點的相位，均是通過對一個子區(qū)域內(nèi)的所有PS點進行相位平均獲得，則上述方法可以視為采用了大量的PS點作為參考點，空間維插值精度可以得到保證。

3 實驗信息

本實驗選定的區(qū)域為一露天開采礦坑(E118°36′23′′,N40°06′44′′)，其位于河北省遷安市馬蘭莊鎮(zhèn)。圖1(a)所示為場景照片，其中紅色橢圓所示為場景的形變區(qū)域。圖1(b)所示為場景衛(wèi)星圖，礦坑正上方為橢圓形，其長軸約1050 m，短軸約680 m，其中黃色矩形代表雷達(dá)的布放位置，雷達(dá)的觀測視角范圍為 60°。礦坑邊坡為典型巖質(zhì)邊坡，基本無植被覆蓋，最大開采深度約400 m，邊坡傾角為38°～47°[14,15]。

實驗中采用一部MIMO雷達(dá)對該礦坑進行了監(jiān)測，連續(xù)獲取了460幅雷達(dá)圖像，時間從2018年10月26日17:30～2018年10月27日11:30。該MIMO雷達(dá)采用16個發(fā)射天線構(gòu)成兩個密集子陣列，16個接收天線構(gòu)成一個稀疏子陣列，如圖2所示。其工作在Ku波段，波長為λ=1.86 cm，等效合成孔徑長度為1.138 m，角分辨率為0.466°[16]。

圖3(a)所示為該露天礦坑在極坐標(biāo)系下的成像結(jié)果，邊坡區(qū)域內(nèi)像素點的幅值主要分布在-30～0 dB的范圍內(nèi)?；诜入x差法進行PS點的選擇，設(shè)置幅度離差門限0.15，幅度門限-25 dB，可篩選出61764個PS點，如圖3(b)所示?？梢钥闯觯麦w上大部分像素點的幅度穩(wěn)定性很高，僅中間的道路上由于雷達(dá)觀測視角的原因，未能有效選擇出PS點。

圖1 場景信息Fig.1 Scene information

圖2 MIMO雷達(dá)系統(tǒng)照片F(xiàn)ig.2 Photo of the MIMO radar system

以第1幅圖像為主圖像，最后一幅圖像為輔圖像，獲取差分干涉相位圖，如圖4(a)所示，將該相位圖反投到3維地形圖上，如圖4(b)所示。借助于3維地形，可以更加準(zhǔn)確地確定形變區(qū)域的發(fā)生位置。圖4(c)所示為PS點的干涉相位圖，顯然大部分PS點的相位在0 rad左右，而紅色橢圓標(biāo)識出的一塊區(qū)域，其干涉相位與其他區(qū)域的明顯不同，該區(qū)域為發(fā)生了較大形變的形變區(qū)域。圖4(d)所示為PS點的干涉相位隨其斜距變化的分布圖，如果將大氣相位建模為隨斜距線性變化的相位分量，其中紅色實線所示為線性大氣相位估計結(jié)果，顯然存在很大的誤差。因此，對較長時間基線的干涉相位圖進行大氣相位補償時，如果不能合理的剔除形變PS點，將會嚴(yán)重影響到對大氣參數(shù)的估計。在處理時序GB-InSAR圖像時，為提高大氣相位補償?shù)臏?zhǔn)確度，進行差分干涉的兩幅雷達(dá)圖像之間的時間基線不宜過大。

圖3 MIMO雷達(dá)圖像與PS圖Fig.3 MIMO radar image and PS map

圖4 長時間基線干涉相位圖Fig.4 Interferometric phase map with long temporal baseline

4 實驗結(jié)果

分析這460幅時序MIMO雷達(dá)圖像時，采用相鄰的兩幅圖像構(gòu)成一個干涉圖像對，則可以獲取459幅干涉相位圖。每一幅干涉相位圖的時間基線僅為2～3 min，通常在這么短的時間內(nèi)，大氣條件的變化很小，相應(yīng)的大氣相位誤差分量也很小。圖5(a)和圖5(b)所示為干涉相位圖A和其相位散點圖，所有PS點的干涉相位均在0 rad左右。在采用線性斜距模型補償大氣相位時，最大的補償分量不超過0.05 rad。但在部分時間段，大氣在空間上不再是均勻變化，導(dǎo)致大氣相位呈現(xiàn)出非線性變化。圖5(c)和圖5(d)所示為干涉相位圖B和其相位散點圖，顯然PS點的干涉相位變化比較復(fù)雜，采用線性斜距模型補償大氣相位時，最大會帶來約0.5 rad的補償誤差，相應(yīng)的形變測量誤差約為0.74 mm，會嚴(yán)重影響到MIMO雷達(dá)的形變測量精度。即使是采用其他模型對相位圖B進行分析，依然無法有效的補償大氣相位。

對459幅干涉相位圖進行相位解纏，并采用線性斜距模型進行大氣相位補償，然后計算每一個PS點的補償后相位序列的標(biāo)準(zhǔn)差。圖6(a)所示為所有PS點的相位標(biāo)準(zhǔn)差圖，顯然形變PS點的標(biāo)準(zhǔn)差很大，設(shè)置0.3 rad的標(biāo)準(zhǔn)差門限，篩選出4341個PS點。這些PS點的分布情況圖如圖6(b)所示，除紅色橢圓形變區(qū)域內(nèi)的形變PS點外，還有少量的形變區(qū)域外的噪聲PS點也被選擇，如紫色橢圓所標(biāo)識出的部分點。

經(jīng)過相位標(biāo)準(zhǔn)差篩選后，得到57423個穩(wěn)定PS點，設(shè)定每個子區(qū)域中的PS點平均數(shù)量為200個，則可以劃分出287個子區(qū)域。以圖5(c)所示的干涉相位圖B為例來說明本文所提方法，圖7(a)中紅色方形點標(biāo)識出了每一個子區(qū)域中心，即287個控制點，以干涉相位圖B作為背景，圖7(b)為局部放大圖，其中每一個黑色多邊形，代表各個子區(qū)域中PS點的最小外接多邊形。之后對于各個子區(qū)域中的PS點進行相位平均，獲取控制點的平均相位，采用反距離加權(quán)法估計所有PS點的大氣相位，如圖7(c)所示，可以很直觀地看出，估計圖中PS點的相位隨其斜距呈現(xiàn)出非線性變化，和圖5(c)高度相似。圖7(d)所示則為補償結(jié)果，所有PS點的補償后相位均在0 rad左右，空變性的大氣相位得到了有效的補償。

圖5 短時間基線干涉相位圖A和BFig.5 Interferometric phase maps of A and B with short temporal baselines

圖6 形變PS點選擇Fig.6 Selection of deformation PS

圖7 非線性大氣相位補償Fig.7 Non-linear atmospheric phase compensation

采用本文所提非線性補償方法和常規(guī)的線性斜距模型補償方法，對這459幅相位圖分別分析，然后從不同斜距處選出5個幅度離差指數(shù)最小的PS點作為參考點，來對比說明本文方法的有效性。圖8(a)所示為采用本文方法獲取的累積相位圖，由459幅補償后的干涉相位圖累加獲取，該幅圖像可以用來反映整個監(jiān)測周期內(nèi)，場景中PS點的相位變化情況。形變區(qū)域PS點的相位變化達(dá)到了約-3 rad，非形變區(qū)域PS點的相位變化在±1 rad范圍內(nèi)，直觀上很難看出非形變區(qū)域PS點的相位是否存在空變性。圖8(a)中紅色方形點標(biāo)識出的點1～點5為非形變區(qū)域的參考點。圖8(b)所示為這5個參考點的時序相位變化曲線，點1～點5的相位變化趨勢高度相似，這是由于大氣相位具有較高的空間相關(guān)性。圖8(c)和圖8(b)所示分別為采用本文方法和常規(guī)方法對點1～點5的補償結(jié)果。常規(guī)方法的補償結(jié)果中，點1的相位隨著時間逐漸增大，這說明了大氣相位殘余誤差隨著時間逐漸積累了起來。本文方法的補償結(jié)果中，各點的相位序列變化更為隨機。

圖9(a)和圖9(b)所示為采用常規(guī)方法和本文方法獲取的累積相位圖，為了更直觀地對比場景中PS點的相位變化情況，兩幅圖像中均忽略了形變區(qū)域?？梢院苤庇^地看出，常規(guī)方法的補償相位圖中，PS點的相位依然呈現(xiàn)出明顯的空變性，尤其是圖像中500～600 m范圍內(nèi)的區(qū)域，該區(qū)域PS點的相位達(dá)到了1 rad，相應(yīng)的形變量為1.48 mm。由于大氣相位補償誤差隨著時間積累起來的原因，這部分區(qū)域很容易被誤認(rèn)為在雷達(dá)監(jiān)測周期內(nèi)，也出現(xiàn)了形變。

圖9 累積相位圖Fig.9 Cumulative phase maps

5 結(jié)束語

本文提出了一種適用于時序GB-InSAR圖像非線性大氣相位補償?shù)母倪M方法，解決了常規(guī)方法中多參數(shù)模型無法準(zhǔn)確模擬大氣相位，從而無法對部分干涉相位圖進行有效補償?shù)膯栴}。采用本文方法與常規(guī)方法，對一露天礦坑的460幅時序地基MIMO雷達(dá)圖像分別進行了處理，對比了參考點的時序相位曲線和累積相位圖，驗證了本文方法的有效性。

本文方法還存在著一些不足。首先本文方法要求先采用常規(guī)方法對時序干涉相位圖進行補償，基于PS點的補償后相位序列的標(biāo)準(zhǔn)差大小來選擇穩(wěn)定PS點，但相位標(biāo)準(zhǔn)差門限是人為設(shè)定的；其次如果部分PS點僅在少量圖像中發(fā)生了形變，這些PS點的相位標(biāo)準(zhǔn)差可能較小，無法有效地篩選出；最后本文方法目前是應(yīng)用于對時序GB-InSAR圖像的事后大氣相位補償，實際形變監(jiān)測中更關(guān)注的是實時處理，需要進一步研究實時形變處理中的非線性大氣相位補償方法。