地基大視場SAR形變監測的非均勻大氣相位校正方法

白澤朝 王彥平 王振海 胡 俊 李 洋 林 赟

①(北方工業大學電氣與控制工程學院 北京 100144)

②(北方工業大學信息學院雷達監測技術實驗室 北京 100144)

③(中南大學地球科學與信息物理學院 長沙 410083)

1 引言

滑坡是我國主要地質災害之一，其中特大型滑坡，其滑坡體量巨大，往往會造成嚴重的財產損失和人員傷亡，滑坡地表變形是實現地質災害監測和預警的重要輔助信息。地基干涉合成孔徑雷達(Ground-Based Interferometric Synthetic Aperture Radar,GB-InSAR)是一種高精度的地表形變監測設備，雷達天線在線性軌道上移動，形成方位合成孔徑，有效提升了影像方位向分辨率，通過對同一場景的重復掃描和成像，采用干涉測量方法提取目標區域的形變信息[1,2]。與星載干涉合成孔徑雷達(Interferometric Synthetic Aperture Radar,InSAR)[3,4]相比，GB-InSAR具有觀測角度靈活、測量精度高(可達亞毫米)的優點，是局部區域高精度形變監測的重要監測手段。對于特大型滑坡監測，不僅要求GB-InSAR平臺實現遠距離覆蓋，而且在方位向上具有大視場角度成像能力，來滿足特大滑坡廣域范圍監測需求。

此外，觀測范圍在方位向上的擴大，使得場景內環境更為復雜，特別是觀測中大氣變化的影響，大氣會引起電磁波折射，從而造成電磁波傳播路徑和方向的變化。在兩次成像時刻大氣條件(溫度、氣壓和濕度)的變化會導致不同的傳播延遲，并且在相位中將存在大氣相位分量(通常稱為“大氣相位屏”)[5-7]。目前的GB-InSAR大氣相位校正方法主要分為3類，一類基于氣象數據的方法[5]，依據電磁波傳播的折射率與場景內溫度、氣壓和濕度數據之間的關系來校正大氣相位。一類基于函數模型方法[8-12]，通常大氣影響在空間上具有強相關性，當滿足大氣在空間均勻性假設時，穩定的永久散射體點相位差在距離方向上的投影呈現線性函數關系，通過求取函數關系即可校正大氣相位。此外，考慮了陡峭地形條件下大氣折射率會隨高程發生改變，有學者進一步考慮了高程對大氣延遲系數的影響，給出了距離-高程的多元回歸模型[13,14]或者兼顧水平位置和高程等信息的多參數函數模型[15,16]。另一類空間插值的方法[17-21]，當大氣在空間非均勻變化時，依據穩定的特征點，通過插值的方式進行大氣相位校正。

基于穩定永久散射體(Persistent Scatterers,PS)線性函數模型和距離-高程的多元回歸模型的方法易于實現，并且已廣泛應用在GB-InSAR形變監測領域。然而，在特大滑坡變形監測場景下，需要在方位向上具有較大的視場角度(>100°)，GB-In-SAR系統方位向范圍的擴大，使得觀測場景內環境更為復雜，可能存在明顯的湍流或者局部水汽變化[13]，導致大氣介質在方位方向呈現非均勻的特性[21]。因此，傳統的基于穩定PS點函數模型的方法無法在復雜大氣條件下準確校正大氣相位。

在方位向大視場角度山區復雜大氣條件下，大氣參數的隨機變化，可能在空間上存在線性和非線性兩個特定特征的異常相位波動。針對上述問題，本文提出一種兩階段半經驗模型校正線性和非線性大氣相位。選擇具有方位向較大視場角度和跨越長江觀測的三峽庫區新鋪與藕塘滑坡實驗區，詳細介紹了大氣校正的處理流程，并采用本文提出的兩階段半經驗模型和常規方法進行大氣相位校正實驗，驗證了算法的有效性。

2 兩階段半經驗模型大氣相位校正方法

2.1 PS點選擇

PS點選擇是進行大氣相位校正的必要前提，需要根據實際觀測場景，避免植被波動和水域等造成誤選點情況，選擇合適的選點策略。一般根據點的幅度離散和相位標準差的統計特性，在高信噪比的情況下，目標點的穩定性直接由幅度離差指數DA表示為[22]

其中，δA是幅度標準偏差，mA是幅度平均值。幅度離差指數越小，幅度信息越穩定，然而場景中存在水域時，往往會出現誤選等情況，可以通過估算相干系數來估計每個像素點的干涉相位質量，濾除水域的誤選點，公式為[23]

其中，S1和S2是與干涉圖相同像素所對應的復值，|·|為絕對值運算符，E{}為期望值。將符合幅度離差閾值和相干性閾值的點選作PS點。

2.2 距離-高程模型初補償

GB-InSAR目標點的干涉相位可以表示為[1]

其中，k為整數表示模糊度，φdef表示形變相位分量，φatm表示大氣相位分量，φnoise表示誤差相位分量。為了獲取高精度的形變相位φdef，需要有效的大氣相位φatm校正方法。

通過沿路徑L積分得到雷達傳感器與目標點之間大氣相位為[10]

其中，n表示折射率的變化，它隨時間t和距離γ而變化，L表示信號的傳輸路徑，λ是波長。一般情況下，假設大氣在空間上是勻質的，n在距離γ上不發生變化，在時間t上隨機變化，因此可以將干涉圖中大氣相位分量表示為隨距離變化的公式[8]

其中，β0為常數項，β1表示與距離相關的線性系數，γ表示傳感器與目標點之間的距離。另外，考慮到地形對大氣相位的影響，兼顧距離和地形的大氣相位可以表示為[16]

其中，β0為常數項，β1表示與距離相關的線性系數，β2表示與距離、地形相關的系數，γ表示傳感器與目標點之間的距離，h表示目標點的高程值。可以建立大氣相位與距離和地形的方程組

其中，φ為n個PS點干涉相位構成的n×1維向量，X為常數1和n個PS點距離與高程值構成的n×3維向量，β為待估計的3×1維 向量，ε為相位誤差構成的n×1維 向量。采用最小二乘對β進行求解

則大氣相位的估計值為

將干涉相位φ減掉估計的大氣相位φatm即為大氣校正后相位。在線性大氣相位校正階段，基于相干性和幅度離差進行PS點選擇，設置較高的閾值選作高質量PS點，采用距離-高程函數模型估算大氣相位。此外，高質量PS點中可能存在一些不可靠點，比如形變區的點和相位跳變的點，為了確保估算結果的準確性，將大于二倍相位標準差σ的點剔除，迭代求解距離-高程函數模型參數，重新校正大氣相位。

2.3 空間維插值補償

在校正隨高程變化的大氣相位后，計算觀測場景內高質量PS點形變，設定形變閾值掩模選擇出穩定PS點。穩定PS點的相位主要包含大氣相位和噪聲相位，大氣相位在空間上表現出非勻質特征，但在較小的距離范圍內可以視為勻質特征，噪聲相位在空間上表現為隨機特征，不具備空間相關性。在一定的空間范圍內進行均值濾波，則可以有效濾除噪聲相位影響，同時基于穩定的PS點，采用反距離權重插值算法，插值出采用較低閾值選取的PS點大氣相位，從而有效校正非線性大氣相位，本文所研究的高質量和低閾值設置在4.2節詳細介紹。其中反距離權重插值算法是利用已知點與待插值點的距離來定義權重，根據周邊的已知點加權計算待插值點的相位值，距離越近權重越大，其公式為

其中，Z(x,y)為插值點結果，(x,y)為待插值點的空間坐標，Zi為第i個已知點相位值，|di|為第i個已知點與待插值點的空間距離，μ為權重的冪指數，一般取2[20]。此外，考慮到大氣相位的空間相關性和計算效率，只選取待插值點周邊最近的3個已知點參加計算。

本文算法的大氣相位校正流程如圖1所示，首先，利用幅度離差閾值和相干性閾值選取PS點，設置較高的閾值選作高質量PS點，參與后續的線性大氣相位估計。然后，根據經驗判斷干涉圖中大氣誤差影響范圍和特點，基于高質量PS點采用距離-高程的函數模型對干涉相位校正高程相關的大氣影響，在計算中剔除大于二倍相位標準差PS點，迭代求解距離-高程函數模型參數。最后計算序列圖像形變，設定閾值選擇出穩定PS點，采用反距離權重插值算法估計低閾值選取的PS點大氣相位，從而有效校正非線性大氣相位。

圖1 本文算法大氣相位校正流程Fig.1 The workflow for atmospheric phase screen processing steps

理論上，在地形復雜區域距離-高程的函數模型可以提高大氣校正的準確性。然而，在實踐中干涉圖可能同時存在不同特征的大氣相位，簡單的函數模型可能校正效果無法滿足監測需求。為了在復雜大氣條件下保持形變測量的可靠性，拓展GB-In-SAR在大視場觀測的應用場景，本文方法同時兼顧了線性和非線性大氣相位，實現不同特征大氣相位的有效校正。

3 研究區描述和系統參數

3.1 研究區描述

研究區為新鋪和藕塘滑坡實驗區，位于長江三峽庫區重慶市奉節縣安坪鎮。圖2為研究區范圍，其中雷達布設在長江北岸，跨江觀測長江南岸的橘色半橢圓區，觀測區長度約為3500 m，寬度約為2500 m，雷達在方位向上的觀測視角為120°。新鋪滑坡位于長江南岸(右岸)斜坡地帶，北抵長江河床，南至泰山廟山脊，總體斜長約2 km，相對高差約300 m，斜坡坡度15°～20°。藕塘滑坡位于新鋪滑坡的西側，屬淺中切割單斜低山河谷地貌，巖層傾向與坡向近于一致，長江從藕塘北邊由西流向北東，與巖層走向夾角為10°～15°。

圖2 研究區范圍Fig.2 The study area

3.2 地基大視場SAR系統

實驗采用自主研制的地基大視場SAR系統對新鋪和藕塘滑坡進行了觀測，連續獲取了29幅雷達圖像，時間從2021年7月27日17:44～22:56。該GBInSAR系統采用角度關聯合成大視場成像技術，實現方位向120°以上大視場角度范圍的形變監測。如圖3所示，其工作在Ku波段，波長為1.74 cm，其方位分辨率在1 km處為5.4 m，距離分辨率為0.37 m。

圖3 GB-InSAR系統Fig.3 The GB-InSAR system

4 實驗分析

4.1 大氣相位分析

本文提出的大氣相位校正模型應用于三峽庫區新鋪和藕塘滑坡形變觀測中，并對大氣校正效果進行了評估。首先，對相鄰圖像進行干涉處理，每個干涉圖的時間基線約為10 min。根據觀測區情況，選擇同時滿足幅度離差指數和相干系數的點作為PS點，其中設置幅度離差指數0.15，相干系數0.9，篩選出25837個高質量PS點。圖4(a)和圖4(b)分別為平均幅度圖和平均相干系數圖，圖4(c)為高質量PS點的幅度離差圖，可以看出坡體上大部分像素點的幅度穩定性和相干系數都很高，不存在水域和植被區的誤選點。

圖5顯示了受大氣影響小的干涉圖，圖6和圖7分別顯示了受到隨高程變化的大氣和非勻質大氣影響的干涉圖，使用不同的色標表示-π～π范圍內觀察到的相位變化。干涉對A的影像采集時間分別為19:10與19:20，干涉對B的影像采集時間分別為20:04與20:14，干涉對C的影像采集時間分別為19:53與20:04。在數據采集的間隔時間內，圖5(a)中的干涉相位沒有超過色標范圍，也沒有顯示出重復的色標變化，圖5(b)的高質量PS相位散點也都在0 rad附近波動，表明在2500 m的最大范圍內沒有觀察到相位纏繞現象。另外，由于GB-InSAR的成像幾何，可能在雷達近處存在有強目標的干擾，會在部分干涉圖零方位角上出現一條相位變化十分均勻的條紋，后續通過PS點選擇可以避免對大氣相位校正的影響。

圖6 干涉對B相位Fig.6 Interferogram B

大氣引起的相位變化如圖6(a)和圖7(a)所示，可以將大氣根據空間分布分為兩種不同情況。圖6(a)為受隨高程變化大氣影響的干涉相位，從其相位散點圖6(b)中可以看出大氣在1000～1500 m間具有恒定的相位遞減，由于大氣條件變化存在明顯大氣相位分量，該距離范圍內干涉圖中高質量PS點相位值位于正相位區域的0～1 rad。相位變化沿距離和地形具有恒定的遞減，可以通過應用距離-高程模型對大氣進行校正。

如圖7(a)所示，干涉對C存在更為復雜的非勻質大氣影響。由于地基大視場SAR觀測場景的擴大，在受到風或者湍流大氣影響時，場景內勻質大氣的假設不再成立，從圖7(b)相位散點可以看出，大氣在1000～1500 m處呈現0～1 rad的相位遞減，而在1500～2500 m處存在明顯的非線性相位波動，因此采用傳統的距離-高程的大氣相位校正方法建模將難以有效去除該區域大氣影響。

傳統的基于函數模型的大氣校正方法僅適用于小范圍的勻質大氣，即當氣壓、溫度和濕度沿距離或者地形均勻變化時。因此，傳統的大氣校正方法不能作為所有干涉圖大氣校正的通用方法，必須設計一套能夠處理復雜大氣條件下的算法，使得可以適應地基大視場SAR觀測場景，尤其是跨江等水汽復雜的觀測中。

4.2 大氣相位校正效果分析

由于觀測范圍廣且地形復雜，大氣條件變化明顯，相應的大氣相位誤差也較大(如圖6和圖7)。在第1階段采用距離-高程模型校正大氣相位，圖8和圖9為在1000～1500 m范圍校正隨高程變化大氣后的干涉圖和對應高質量PS點相位散點圖，PS點相位在0 rad波動。與圖6和圖7對比，干涉對B和干涉對C都有效校正了距離-高程相關的線性大氣相位，最大的校正相位分量約為1 rad。但是，由于環境復雜，大氣在空間上呈現不均勻變化情況，即使校正線性大氣相位后，兩個干涉圖仍然存在非線性大氣相位變化。如圖9(b)所示，隨著距離的變化，高質量PS點的相位呈現非線性變化趨勢。

圖7 干涉對C相位Fig.7 Interferogram C

圖8 距離-高程校正后干涉對B相位Fig.8 Interferogram B after distance-elevation correction

圖9 距離-高程校正后干涉對C相位Fig.9 Interferogram C after distance-elevation correction

在第2階段采用非線性大氣相位校正方法，對選取的高質量PS點計算累計形變量，設定形變閾值為5 mm，選擇出穩定PS點，然后選取同時滿足幅度離差指數0.25和相干系數0.8的75104個較低閾值的PS點，基于穩定的PS點通過插值的方式估算所有PS點的大氣相位。如圖10和圖11為校正后的干涉對B和干涉對C，在校正前干涉對B在沿江邊距離向約1000 m位置存在約1 rad相位波動，在距離向約1500 m位置存在約-0.7 rad相位波動(圖10(a))。經過非線性大氣校正后，補償了明顯的大氣相位，干涉對B整體相位在0 rad附近波動。另外干涉對C存在更為復雜的大氣相位，在校正前，不僅沿江邊約1000 m,1500 m位置存在明顯大氣相位，在距離向約2500 m最遠距離位置也存在約-0.5 rad相位波動(圖11(a))，非線性大氣校正后，干涉對C得到了有效相位補償。

圖10 非勻質大氣校正后干涉對B相位Fig.10 Interferogram B after correcting nonlinear atmospheric

圖11 非勻質大氣校正后干涉對C相位Fig.11 Interferogram C after correcting nonlinear atmospheric

與小范圍觀測不同，GB-InSAR在遠距離、大視場的觀測中顯示出更為復雜的大氣相位變化，但可以通過選取穩定的有效測量點，提高形變信息的準確性。圖12為獲取的觀測區累計形變量，如圖12(a)所示，未校正大氣相位時在江沿岸約1000 m和遠距離2000～2500 m位置處存在明顯大氣相位引起的約5 mm形變誤差，如圖12(b)所示，僅使用距離-高程函數方法在沿江和遠距離處大氣相位有所削弱，但沿江邊區域仍然存在大氣相位，造成了約-3 mm的形變誤差，如圖12(c)所示，使用本文方法有效補償了沿江和遠距離區域大氣相位，表明本文算法的大氣相位校正結果相位誤差顯著降低。理論上，在復雜的場景根據大氣相位的特征本文方法可以提高大氣估計的準確性和形變提取的魯棒性。從累計形變量圖分析，監測區較穩定，在距離向約2000 m位置處，如圖12(c)中紅色矩形框所示，存在約4.2 mm的形變區，結合已有研究[24]，推測在監測時間段內，滑坡的后緣區存在明顯的形變，即在藕塘滑坡存在向下滑動的情況。此外，在形變區外存在零星分布的少量形變點，推測可能是植被擾動造成的。

圖12 累計形變量Fig.12 Cumulative deformation

為驗證大氣相位校正的準確性，挑選了在不同距離向和不同滑坡位置具有代表性的4個像元點變形曲線進行分析，圖12(c)中標記了像元點位置，圖13為像元點變形曲線圖。P1位于新鋪滑坡的高程較大位置，經過距離-高程函數校正后反而誤差更大，最大誤差約為2.4 mm。P2位于藕塘滑坡的高程較大位置，在未進行大氣相位校正前，形變的波動幅度較大，最大誤差約為-1.6 mm，經過距離-高程函數校正后，形變波動有所變小，但仍然存在約為-1.1 mm誤差，本文方法校正后，變形序列變化更為隨機且在0.5 mm范圍內波動。P3位于監測到的藕塘形變區，總體呈現向下滑動的趨勢，由于距離-高程函數方法主要校正沿江和遠距離區域隨高程變化的大氣相位，該區域未進行大氣校正，而本文方法大氣校正后沒有將真實的形變補償掉，仍然保留了原有的形變趨勢。P4位于新鋪滑坡靠近長江的沿岸，經過距離-高程函數校正后最大誤差由1.5 mm減小為1.1 mm，本文方法校正后，變形序列變化更為隨機且在0.5 mm范圍內波動。因此，本文方法不僅保留了真實的形變區，而且在非形變區實現了誤差的有效降低。

圖13 典型點變形曲線對比Fig.13 Cumulative deformation map

通過三峽庫區新鋪和藕塘特大滑坡的地基大視場雷達圖像實驗，驗證和分析了本文方法的有效性與準確性。在校正江沿岸和遠距離位置隨高程變化的大氣相位，本文方法可以達到與距離-高程函數方法相同甚至更優的精度。在校正非均勻大氣相位，函數模型方法往往會失效，采用大量的穩定PS點作為參考點，空間維插值的方式，可以有效保證非均勻大氣相位的校正精度。

5 結語

本文采用地基大視場SAR系統識別新鋪和藕塘滑坡的形變區，觀測到了由大氣條件變化引起的干涉相位變化。由于該地區地形和大氣條件復雜，傳統的大氣相位校正方法通常會提供不準確的形變結果。復雜大氣環境下的干涉相位不僅存在隨距離-高程變化的線性大氣相位成分，還存在非線性大氣相位成分。因此，提出了一種兩階段半經驗模型估計各層大氣相位，以消除空間域中的大氣相位變化。該方法基于空間模型和統計補償方法，在不使用溫度、氣壓或濕度等大氣參數的情況下，對空間非均勻性大氣相位進行補償，對于場景環境信息較少的GB-InSAR觀測，可以有效地消除大氣效應和準確地獲取形變信息。本文方法不需要在站點和實驗階段收集與處理氣象數據，拓展了GB-InSAR的適用性且能更好地對大范圍下的大氣誤差進行校正。

本文方法還存在著一些不足，由于地基SAR系統工作在Ku波段，容易受到植被擾動影響，在形變區外仍然存在零星分布的由于植被擾動造成的誤差點，需要在后續工作中進一步研究如何避免植被擾動帶來的誤差。