基于相干系數的DInSAR地表形變信息提取方法

樂 穎，夏元平，劉媛媛，胡 昊

(東華理工大學測繪工程學院，江西 南昌 330013)

0 引言

20世紀70年代左右，合成孔徑雷達干涉(interferometric synthetic aperture radar，InSAR)技術出現，Rogers A E等于1969年首次利用 InSAR 技術研究了金星表面的地形特征[1]。隨后，Gabriel等于1989年首次提出合成孔徑雷達差分干涉測量(Differential InSAR，DInSAR)技術[2]。該技術是基于InSAR技術提取地表形變信息的一門極具潛力的新興技術，由于在一定程度上克服了光照和天氣的影響，可全天時全天候監測，具有高精度、大區域、快速等優勢，因此受到了學者的廣泛關注[3]。目前，DInSAR技術精度甚至可達毫米級，己經成為地震、火山、滑坡、地表形變信息等應用領域的主要監測手段[4]。

高分辨率SAR衛星的快速發展，使DInSAR技術已經成為目前獲取地表形變信息的研究熱點[5]。文獻[6]比較了兩種配準方法的特點，并將DInSAR技術應用于神木礦區的沉陷監測，獲取最新的礦區沉陷信息及對應的時間序列。文獻[7]利用DInSAR技術監測呂梁山區某礦區地表，通過兩個剖面數據分析地表沉陷規律，表明能夠得到對應的形變信息。文獻[8]將DInSAR技術運用于陜西彬長礦區，并與實際情況相結合，有效獲取了礦區的形變信息。文獻[9]基于SAR 影像的振幅信息提出了一種改進的PS探測方法，通過對上海市局部區域進行研究，準確地提取了研究區域的形變信息。文獻[10]基于DInSAR技術，結合GPS水平形變結果，可快速獲得研究區地表三維形變，為后期的研究提供有效的參考。文獻[11]針對不均勻地表形變問題，利用TS-DInSAR(time series DInSAR)技術對Sentinel-1A數據進行研究，結果表明該方法對大范圍地表形變監測具有可靠的精度。然而，以上文獻為了校正衛星軌道和相位偏移，去除處理過程中殘余相位，在運用DInSAR技術處理影像時憑借人為主觀經驗選取控制點進行軌道精煉，這將影響相位信息轉化為形變信息的準確性。若研究區域缺少地面控制點信息時，隨意選取控制點將影響變形信息結果提取的精度。本文針對此問題，提出基于相干系數的DInSAR地表形變信息提取方法。

1 DInSAR技術

1.1 基本原理

DInSAR形變監測技術是經由InSAR技術衍生發展而來，主要包括二軌、三軌與四軌三種。其中，“二軌法”差分干涉技術運用領域最廣，其表達式如下[12]：

φ=φflat+φtopo+φdefo+φatm+φnoise

(1)

式(1)中，φflat為平地相位；φtopo為地形起伏引起的地形相位，可以將地形信息恢復；φdefo為地表發生形變造成的相位；φatm為大氣發生延遲的大氣相位；φnoise為觀測時產生的噪聲相位[12]。簡而言之，就是將式(1)中φflat、φtopo、φatm與φnoise等相位去除，最后只保留φdefo形變相位的技術[13]。

1.2 數據處理流程圖

由于雷達數據與光學影像不同，屬于主動遙感，影像包含的是振幅、相位和極化等多種信息。目前，常用的技術方法主要有合成孔徑雷達干涉、極化分析、幅度追蹤、層析建模和立體量測等五種方法，而其中合成孔徑雷達干涉方法是最常用也是運用范圍最廣的技術。合成孔徑雷達差分干涉通過處理覆蓋同一地區不同時刻獲取的兩幅SAR影像的相位信息，能夠測量地表厘米級的形變。在DInSAR數據處理過程中，主要包括：影像配準、干涉圖生成、干涉圖自適應濾波和相干性計算、最小費用流相位解纏、控制點選取、軌道精煉、重去平、相位轉形變及地理編碼等，每個數據處理步驟過程中產生的誤差對干涉結果的精度和可靠性都具有顯著影響，具體的處理流程見圖1。DInSAR地表形變監測的數據準備主要包括SAR數據選取，外部DEM和精密軌道數據準備。其中，外部DEM是為了利用其模擬地形起伏相位分量，精密軌道數據則可以有效削弱SAR衛星軌道誤差導致的殘余參考相位。

圖1 DInSAR數據處理流程圖Fig.1 DInSAR data processing flow chart

2 基于相干性系數的 DInSAR 地表形變信息提取方法

2.1 干涉相干性測度

相干圖主要用于體現兩幅影像的相關程度，可以用于評價SAR影像精確配準后所得到干涉條紋圖的質量，主要是通過計算干涉條紋圖的相干系數實現。相干系數最初由普拉蒂于1993年提出，相干系數的取值范圍為[0,1]，值越接近1說明相干性越高，其理論模型表達式如下[14]：

(2)

式(2)中，E[·]表示數學期望，u*表示共軛復數，u1與u2為主副圖像的信號。在實際的InSAR處理過程中，在同一像素中無法得到計算所需數量的采樣值。因此，根據式(2)的理論模型，基于SAR影像復數數據計算獲得相干系數的標準表達式為[15]：

(3)

式(3)中，u1(n,m)，u2(n,m)分別表示主、副影像數據塊內某個坐標(n,m)處的復數值；|?|2表示對應的二階范數；M與N表示計算相干性的數據塊尺寸大小；m與n表示數據塊內對應的行列號。

2.2 軌道精煉控制點的選取

利用DInSAR技術獲取地表形變信息時，引入穩定的地面控制點可以有效估計殘留相位并去除平地效應，從而提高結果的可靠性。一般而言，選取控制點要遠離形變區，認為形變為0，陡峭的地形區域和有殘余地形相位區域，地形起伏大的山區，最好是選擇山谷底部的平地區域[16]。為了更加快速準確地選擇軌道精煉的控制點，本文基于相干系數閾值選取所需的控制點，相干性系數圖的像元分布情況見圖2。

從圖2可看出，2015—2020年干涉對像元個數在相干性系數為0.5時達到了峰值，2015、2017和2020年的整體相干系數數值比其他年份高。由于相干系數圖上的像元個數數量眾多，選取軌道精煉時需要的像元數量很少。相干性閾值設置過低會導致相關性高像元數目劇增，大大增加數據計算的時間，從而增加后續篩選工作的難度。因此，在保證高相干性又考慮數據量大小的基礎上，通過反復對比像元個數進行設置合適的閾值，最終將相干系數閾值設定為0.98，每年干涉相干性系數的像元分布情況見表1。

圖2 2015—2020年相干性系數圖像元分布情況Fig. 2 Distribution of image elements of coherence coefficient from 2015 to 2020

表1 每年干涉相干性系數的像元分布情況

2.3 控制點選取驗證

為了更系統地驗證本文方法選取控制點的可靠性，首先統計2015—2020年期間每年干涉相干性系數圖符合要求的像元數量，然后根據相干系數閾值分別選取25～35個軌道精煉控制點，最后計算各點絕對殘余誤差，其計算公式為：

(4)

式(4)中，γ是干涉相干性測度，H是高程模糊度。

以2019年為例，對比分析DInSAR方法與本文方法用于選取軌道精煉控制點的絕對殘余誤差，結果見圖3。本文方法用于選取軌道精煉控制點整體的絕對殘余誤差較小，只有個別點的絕對殘余誤差偏大，誤差符合要求；而DInSAR方法選取控制點的絕對殘余誤差整體偏大，后續需要篩選絕對殘余誤差較大的點。此外，文中為了對比分析DInSAR方法與本文方法用于選取軌道精煉控制點的準確性，引入了標準差指標，并計算得到兩種方法絕對殘余誤差的標準差分別為1.33和0.36，而指標值越小說明準確性越高。因此，在缺少地面控制點信息或研究區形變先驗認識的情況下，若需要選擇可靠的地面控制點，可借助相干性系數圖進行剔除低相干性點。

圖3 2019年選取軌道精煉控制點絕對殘余誤差Fig.3 Absolute residual error at the control point of orbit refining in 2019

3 實例驗證

3.1 研究區域概況

贛州市的稀土儲量占全國最多，而定南縣是贛州市離子型稀土主產縣之一，足以表明定南縣稀土礦礦區扮演著重要的作用。基于此，本文以定南縣作為研究區域，其地理坐標為東經114°47′49″E～115°22′48″E，北緯24°33′37″N～ 25°03′21″N，總面積為1 321 km2。縣區內分布眾多河流，主要隸屬于定南水和桃江兩水系[17]。同時，定南縣地處南嶺與武夷成礦帶的中間部位，故而擁有良好的成礦條件，礦產資源非常豐富。在礦山開采過程中，露天礦坑開挖會直接破壞地表植被，嚴重影響著當地的可持續發展和生態功能恢復。因此，研究并提取定南縣的地表形變信息尤為重要[18]。

3.2 實驗數據

由于雷達數據與光學影像不同，屬于主動遙感，影像包含的是振幅、相位和極化等多種信息。目前，合成孔徑雷達干涉方法是最常用也是運用范圍最廣的技術。本文使用的雷達影像數據為C波段的Sentinel-1A衛星的SLC數據，選取了時間跨度為2015—2020年期間每年10月份共12景的影像數據，具體雷達影像數據的基本參數見表2。

表2 雷達影像數據的基本參數

3.3 結果分析

通過對12景覆蓋研究區的SAR影像數據進行DInSAR處理，借助相干系數閾值獲取相干性高的像元，并將其作為軌道精煉的控制點引入后續的處理中，最終可得到贛州市定南縣鄉鎮在2015—2020年的地表形變信息，其形變結果圖如圖4所示。監測結果表明：研究區在監測時段內地表形變變化相對穩定，整體形變量主要集中在-0.07～0.1 m，部分區域出現了一定范圍的地表沉降，沉降量集中在-0.11～0.11 m，定南縣在2017、2019和2020年內的地表形變較2015、2016和2018年更穩定。

定南縣共有7個鎮，分別為鵝公鎮、黃香鎮、巋美山鎮、老城鎮、歷市鎮、嶺北鎮和天九鎮。根據雷達SAR影像數據處理得到的結果，得到定南縣各鎮每年的地表形變變化情況，所得結果如圖5所示。從圖5的監測結果可以得出：定南縣各個鄉鎮在2017、2018、2019和2020年內變化趨勢基本一致，2015年形變量較大的鄉鎮為嶺北鎮，形變量達到了-25.26 mm；2016年形變量較大的鄉鎮為嶺北鎮、黃香鎮和巋美山鎮，形變量分別為5.26、5.26和26.29 mm；整個研究區域在監測時段內最大沉降量和最大上升量分別為-84.82和84.816 mm，分別位于嶺北鎮和鵝公鎮。經過查找資料得知，贛州市定南縣內具有豐富的離子型吸附性稀土礦礦產資源，且這些礦區位置大多位于嶺北鎮，其中包含了木子山、甲子背、長坑尾等稀土礦區[17]。由于嶺北鎮的平均形變量與其他鄉鎮的平均形變量相差-20.20～-10.10 mm，稀土礦的開采可能會引起少量的地表形變。

圖4 2015年-2020年研究區域地表形變結果圖Fig.4 Results of surface deformation in the study area from 2015 to 2020

圖5 2015—2020年定南縣鄉鎮地表形變信息Fig.5 Surface deformation information of towns in dingnan county from 2015 to 2020

為了探尋稀土礦的開采與地表形變之間的關系，并獲取稀土礦開采區域的地表形變細節信息和地表形變發展趨勢信息，本文結合稀土礦權在已有礦權范圍內尋找了一些相關性高且能代表形變規律的特征點，特征點選取分布見圖6，最后對特征點進行時間序列分析。

分別計算選取的30個特征點在2015—2020年間的地表形變量，圖7繪制了地表形變細節信息和發展趨勢。

圖6 30個特征點選取分布圖Fig.6 Distribution of 30 feature points

圖7 30個特征點的地表形變情況Fig.7 Surface Deformation of 30 feature points

從圖7中不難發現，選取的30個特征點在同年份的變化趨勢基本一致，在2015年的地表形變發展整體為下降趨勢，沉降量范圍在-25.5～-15.5 mm，結果表明這些稀土礦區均存在開采現象；2016年和2017年的地表形變在0附近波動，地表變化幅度分別為-5.26～15.77 mm和-12.82～4.27 mm，由于這段時間國家開始加大對非法開采行為的打擊力度，使得一些稀土私礦被查封；2018年整體呈下降趨勢，部分點為向上抬升，地表形變在-32.984～14.136 mm之間，贛州地區對稀土礦山環境進行了綜合治理以解決歷史遺留問題，將稀土礦區平整為建設用地投入到當地經濟發展中去；2019年和2020年的地表變化幅度分別為-26.03～18.07 mm和-3.81～3.81 mm，說明贛州定南地區開展的廢棄稀土礦山治理取得效果顯著，對當地的生態環境有了明顯的修復。監測結果表明：由于離子型吸附性稀土礦的礦區范圍不大，而且分布稀疏，從而導致稀土礦的開采引起的地表形變較緩慢。由于稀土礦的開采嚴重破壞了生態環境，江西贛州自2013年來不斷推進稀土產業整頓，從加大稀土礦非法開采行為打擊力度到對生態環境進行綜合治理，很大程度上促進了礦山生態環境保護與稀土資源開發利用的協調發展。

4 結論

本文提出了基于相干性系數的 DInSAR地表形變信息提取方法。該方法通過設置一定的相干系數閾值，篩選出適合的高相干性像元作為控制點進行軌道精煉，最后提取出地表形變信息。以贛州市定南縣作為研究區域，借助 2015—2020 年間的Sentinel-1A衛星數據，利用相干性系數形變信息提取方法，監測出研究區范圍內-25.5～18.07 mm的地表形變信息。實驗結果表明，該方法相較于常規的DInSAR方法，用于提取控制點的標準差提高了72.9%，定量表明了基于相干性系數的DInSAR方法用于選取控制點的準確性更高，能夠實現對研究區域地表形變信息的動態監測。由于本文采用的數據量不足，研究過程中沒有考慮大氣因素的影響，導致得到的結果仍表現出一定的偏差。將雷達影像數據與光學遙感影像相結合，并與外業水準測量結果進行對比，將是下一步要進行的工作。