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

胡 程 鄧云開 田衛明 曾 濤

①(北京理工大學信息與電子學院雷達技術研究所 北京 100081)

②(北京理工大學衛星導航電子信息技術教育部重點實驗室 北京 100081)

1 引言

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

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

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

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

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

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

2 補償方法

2.1 穩定PS點選擇

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

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

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

大氣相位φatm可以建模為

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

其中，β0為常數分量，β1表示與斜距相關的線性系數，R表示雷達與目標點之間的距離。基于PS技術進行大氣相位補償時，首先建立線性方程組

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

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

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

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

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

2.2 空間維插值補償

2.2.1 子區域劃分

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

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

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

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

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

2.2.2 大氣相位估計

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

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

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

3 實驗信息

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

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

圖3(a)所示為該露天礦坑在極坐標系下的成像結果，邊坡區域內像素點的幅值主要分布在-30～0 dB的范圍內。基于幅度離差法進行PS點的選擇，設置幅度離差門限0.15，幅度門限-25 dB，可篩選出61764個PS點，如圖3(b)所示。可以看出，坡體上大部分像素點的幅度穩定性很高，僅中間的道路上由于雷達觀測視角的原因，未能有效選擇出PS點。

圖1 場景信息Fig.1 Scene information

圖2 MIMO雷達系統照片Fig.2 Photo of the MIMO radar system

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

圖3 MIMO雷達圖像與PS圖Fig.3 MIMO radar image and PS map

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

4 實驗結果

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

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

經過相位標準差篩選后，得到57423個穩定PS點，設定每個子區域中的PS點平均數量為200個，則可以劃分出287個子區域。以圖5(c)所示的干涉相位圖B為例來說明本文所提方法，圖7(a)中紅色方形點標識出了每一個子區域中心，即287個控制點，以干涉相位圖B作為背景，圖7(b)為局部放大圖，其中每一個黑色多邊形，代表各個子區域中PS點的最小外接多邊形。之后對于各個子區域中的PS點進行相位平均，獲取控制點的平均相位，采用反距離加權法估計所有PS點的大氣相位，如圖7(c)所示，可以很直觀地看出，估計圖中PS點的相位隨其斜距呈現出非線性變化，和圖5(c)高度相似。圖7(d)所示則為補償結果，所有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

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

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

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

5 結束語

本文提出了一種適用于時序GB-InSAR圖像非線性大氣相位補償的改進方法，解決了常規方法中多參數模型無法準確模擬大氣相位，從而無法對部分干涉相位圖進行有效補償的問題。采用本文方法與常規方法，對一露天礦坑的460幅時序地基MIMO雷達圖像分別進行了處理，對比了參考點的時序相位曲線和累積相位圖，驗證了本文方法的有效性。

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