結合InSAR的礦區大梯度變形監測技術進展

陳 磊，趙學勝

(中國礦業大學(北京)地球科學與測繪工程學院，北京 100083)

傳統礦區形變監測技術主要有水準測量與GPS測量[1]，可精確獲取地表離散點的三維形變信息；但無法直接獲取面區域連續的形變信息，需布設大量的水準點或地表特征點，且在后期處理中還需對離散點數據進行內插擬合操作，容易引入更多的計算誤差；傳統形變監測技術均需進行實地測量，需要較多的人力和物力，有時可能因為試驗區環境惡劣無法完成實地監測。因此，使用傳統監測技術難以獲取礦區充足的地表形變信息。

合成孔徑雷達干涉測量(InSAR)技術是近幾十年發展起來的一種有效的地表形變監測技術，該技術可以全天時、全天候、大范圍、高精度地獲取地表影像，廣泛應用于提取高精度的DEM數據[2]。差分干涉測量(D-InSAR)技術是基于InSAR技術發展起來的，能有效獲取礦區微小的地表形變信息[3-5]，其理論精度可達到厘米級甚至毫米級。時序InSAR技術的提出是為消除傳統D-InSAR技術中存在的時間去相干、空間去相干及大氣延遲相位對測量精度的影響[6-9]，可高精度提取礦區沉降邊緣微小的地表形變信息[10-12]。但當地表形變梯度較大，超過InSAR技術所能測量的梯度范圍時，會造成干涉圖失相干而無法得到清晰的干涉條紋，影響后續的相位解纏操作，從而無法獲取精確的地表形變信息。因此，亟需引入新的技術來提取礦區大梯度的地表形變信息。

針對InSAR技術提取礦區大梯度形變信息的試驗，國內外許多學者進行了很多研究，經過整理可分為下述兩種：①結合像元追蹤法與InSAR技術；②結合礦山開采沉陷模型與InSAR技術。下面分別對兩種方法進行詳細描述。

1 結合像元追蹤法與InSAR技術

該技術是結合SAR影像的相位信息與幅度信息，采用InSAR技術，分析相位信息，提取試驗區的微小形變值，采用像元追蹤法，分析幅度信息，提取試驗區的大梯度形變值。該技術可有效彌補InSAR技術在大梯度形變區的相位失相干現象。

像元追蹤法是利用SAR影像的幅度信息，通過對兩景SAR影像進行配準，提取研究區方位向與距離向的像元偏移量，再根據影像方位向與距離向的像元分辨率，獲取研究區二維地表形變值。該技術可彌補InSAR技術中由于失相干因素引起的相位缺失現象，可有效提高在大梯度形變區的形變提取精度；且可彌補InSAR技術只能獲取視線向形變，以及視線向形變對水平向位移不敏感的缺陷[13-14]。

Gray等[15]在1998年首次提出采用影像配準方法獲取形變信息的可能，并利用其提取冰川移動信息。結果表明，該方法可提取冰川的大移動形變信息，但由于缺少高程數據和基線信息導致測量精度較低。Michel等[16]在1999年利用互相關技術實現SAR影像的高精度配準，并應用該方法提取美國Landers地震在距離向與方位向引起的形變值，但由于影像分辨率較低，其測量精度明顯比GPS精度要低。后來，國內外許多學者將該方法應用于提取地震、冰川、滑坡等活動引起的形變監測中[17-22]。

Zhao Chaoying等[23]在2013年首次采用像元追蹤法與InSAR相結合的技術監測礦區大梯度形變，利用SAR影像中的幅度信息對影像進行互相關配準得到兩影像間的斜距偏移，建立垂直向形變Δh與斜距偏移值Δr之間的函數關系為

(1)

式中，θ為入射角。該方法忽略水平形變，認為礦區形變主要是垂直向沉降，由此提取礦區大梯度沉降值。Liu D等[24]和Niu Yufen等[25]分別在2014年和2016年通過建立角反射器來提高礦區點目標密度，提高影像配準精度并將試驗結果與水準實測數據進行對比后表明建立角反射器有助于提高測量結果的準確性。

Fan Hongdong等[26]與Ou Depin等[27]結合像元追蹤法與InSAR技術，利用干涉條紋圖將變形區分為大變形區與小變形區兩類，分別提取礦區大梯度形變與微小形變信息。Chen Bingqian等[28]在2015年對該方法進行改進，采用相干系數值對變形區進行分類，首次建立最大可檢測變形值Dmax與相干系數γ間的函數關系為

Dmax=dmax+0.002(γ-1)

(2)

式中，dmax為理論最大可檢測變形值，一般采用相干系數0.5作為大小變形區域的劃分依據。該方法首先利用影像對的干涉條紋圖或相干系數值，將變形區分為大變形區與小變形區兩類。采用像元追蹤法獲取大梯度形變，利用雙線性多項式函數獲取試驗區的軌道相位，并將其去除得到大變形區域的沉降值；采用InSAR技術獲取微小形變值，將兩者結合得到試驗區的地表形變信息。

Huang Jilei等[29]與Yan Shiyong等[30]在2016年對像元追蹤法進行改進，對影像配準后的總偏移Ototal細分為以下幾個部分

Ototal=Odis+Otop+Oorb+Oatt+Oion+Onoi+Ores

(3)

式中，Odis為采礦活動引起的變形偏移；Otop為地形變化引起的偏移；Oorb為衛星軌道誤差引起的偏移；Oatt為傳感器姿態引起的偏移；Oion為電離層引起的偏移；Onoi為噪聲偏移；Ores為總體殘余誤差偏移。為獲取采礦活動偏移Odis，需要將后面6種偏移量去除。基于試驗區的穩定點，建立一個二階多項式來去除軌道誤差、傳感器姿態、電離層影響以及部分地形起伏引起的誤差偏移Oerr，公式為

Oerr=a0+a1x+a2y+a3xy+a4x2+a5y2

(4)

式中，a0,…,a5為未知參數；x與y分別為影像中點位置的行、列號。提取穩定點的偏移值采用最小二乘擬合估算上述5個參數值，確定誤差模型，提取每個像元位置處的誤差偏移。去除得到的誤差偏移之后，利用外部DEM建立一個二階多項式來去除殘余地形誤差以及噪聲引起的偏移Ores，公式如下

(5)

Huang Jilei等[31]在2017年對該方法進行改進，提出一種更優的影像配準方法：采用互相關系數重新定義信噪比(SNR)，并利用估算的最大沉降值確定搜索半徑和搜索范圍，通過最大化SNR來自動提取最優的局部搜索窗大小，最后確定像元偏移的整數部分和小數部分，轉換得到地表形變值。選取陜西大柳塔礦區的模擬影像與雷達影像進行試驗，并將該方法得到的測量結果與實測水準數據進行對比，結果表明該方法得到的測量結果比固定搜索窗大小得到的結果更精確，但該方法需要確定搜索半徑及范圍、重新定義SNR、優化SNR提取最優搜索窗，所需的算法更復雜，花費的時間更多。

綜上所述，該方法將SAR影像的相位信息與幅度信息相結合，可有效地彌補相位信息在失相干區域的不足，且幅度信息不受干涉失相干影響，無需進行相位解纏，可有效提高大梯度形變信息的測量精度。但像元追蹤法仍存在下述幾點不足：①受地表覆蓋物變化的影響，其幅度值會因地表變化而引起較明顯的變化，影響影像配準過程；②對影像的配準精度要求較高，在實際操作中很難滿足其精度要求；③與SAR影像的分辨率有關，目前的SAR影像分辨率大多為米級，且影像配準精度很難達到0.01個像元，因此地表形變測量精度大多為分米級或厘米級，難以提取更高精度的地表形變信息。

2 結合礦山開采沉陷模型與InSAR技術

該技術是結合礦山開采沉陷模型，利用InSAR技術測量得到的微小形變區的點的時序形變值，采用遺傳算法或最小二乘算法擬合估算各沉陷模型的未知參數，利用得到的開采沉陷模型來預測估算礦區大梯度的地表形變信息。

礦山開采沉陷模型中最常用的是基于隨機介質理論的概率積分法[32]與基于地表單點的下沉時間函數模型(如：Knothe模型[33-34]、雙曲線模型[35-36]、Gompertz模型[37-38]、Logistic模型[39-40]、Weibull模型[41-42]與改進的Knothe模型[43-44])。許多學者已針對上述幾種沉陷模型，結合InSAR測量結果來預測和提取礦區的大梯度形變信息。

概率積分法是因其所用的移動和變形預計公式中含有概率積分而得名，且該方法的基礎是隨機介質理論，又被稱為隨機介質理論法[32]。隨機介質理論最初是由波蘭學者李特威尼申于20世紀50年代引入巖層及地表移動的研究，后來中國學者劉寶琛等將其發展成為概率積分地表沉陷預計模型[45]。目前該方法已經成為我國較為成熟且應用較為廣泛的礦區沉降預計方法之一[46-49]。

Wang Xunchun等[50]在2011年針對目前開采沉陷預測和監測方法中存在的問題，首次提出采用概率積分法與InSAR相結合的方法預測和監測礦區開采沉陷值與沉陷時間。Fan Hongdong等[51-52]在2014年和2015年兩次采用該方法提取礦區大梯度形變信息。Yang Zefa等[53-54]在2016年和2017年兩次對該技術進行改進，首次建立概率積分函數模型W與DInSAR測量得到的視線向形變結果dLOS間的函數關系

(6)

式中，θ為入射角；α為方位角；UN為地表水平南北向形變；UE為地表水平東西向形變。分析概率積分函數確定待計算的未知參數集P與采空區的地質參數集WP分別為

(7)

式中，q為下沉系數；b為水平移動系數；tanβ為主要影響角正切值；θ0為煤層開采影響傳播角；s3、s1與s2分別為煤層走向、下山向與上山向的拐點偏移距；m為開采高度；H、H1與H2分別為煤層走向、下山向與上山向的開采深度；αm為煤層傾角；D3與D1分別為采空區的寬度。由此，可建立采空區三維形變值與未知參數和地質參數間的函數關系

(8)

將式(8)代入式(6)中，并將其轉換為矩陣如下

dLOS(x,y)=

(9)

選取工作面沉降邊緣區點位置的視線向形變值，結合遺傳算法來估算概率積分模型的各參數值；最后利用得到的概率積分模型預測估算試驗區的大梯度沉降值。

Diao Xinpeng等[55]在2016年對概率積分法與InSAR相結合的技術進行簡化，將形變區分為大變形區域和小變形區域兩類，大變形區域采用概率積分法，根據工作面的實測數據信息以及工作經驗，確定概率積分模型的各參數，首次建立工作面任意點(x,y)的相位整周數φ(x,y)與概率積分值W(x,y)之間的函數關系

(10)

式中，λ為波長；θ為入射角；Int為取整符號。根據式(10)估算工作面任意點位置的整周相位，小變形區域采用DInSAR技術，獲取干涉影像對的纏繞相位，與概率積分模型得到的對應干涉對的整周相位結合得到任意點的總相位值，將其轉換為總形變值。

Logistic模型是一種生長模型，該模型又被稱為Verhulst-Pearl模型，它在地表沉降監測方面得到廣泛應用[56-57]。Logistic時間函數為[58]

(11)

式中，W0為工作面上某點的最大可能沉降值；a為模型參數；b為時間影響系數。

Lee等[40]在2006年首次利用Logistic模型對韓國某廢棄礦區進行形變監測，并分析模型參數與地表沉降值間的關系。Yang Zefa等[59]在2017年首次將InSAR技術與Logistic模型相結合獲取礦區大梯度形變信息，并首次建立工作面每個點位置的視線向測量結果dLOS與Logistic模型間的函數關系

(12)

為計算式(12)中W0、a和b這3個未知參數，需建立影像干涉對DInSAR技術得到的視線向形變值dobs與Logistic模型中各參數間的函數關系

(13)

綜上所述，概率積分法可同時預測估算整個工作面的形變值，Logistic模型可預測估算工作面中任意點在不同時刻的形變值。概率積分法可采用較少的SAR影像，利用InSAR技術獲取工作面范圍內充足的沉降邊緣點的精確的形變值，建立模型參數的多余觀測方程，估算模型中的未知參數，并利用得到的概率積分模型對礦區大梯度沉降信息進行模擬與預計；Logistic模型可利用InSAR技術提取的工作面某點不同時間段的形變值，建立該點模型的未知參數的多余觀測方程，估算得到該點的開采沉陷模型。但該方法中仍存在下述幾點不足：①概率積分法中一些參數的確定需依據實際的工作面地質條件及開采情況，不同的地質開采條件會對參數的結果造成較大的估算誤差；且當估算的參數與實際情況存在差異時，會對模型后續的擬合估算帶來誤差影響。②Logistic模型在開采沉陷初始時刻的沉降值與沉降速率都不為0，造成在開采初期的估算預測精度較低；且該模型在開采預測過程中為保證精度還需確定幾個特征點，當特征點缺失時對該模型擬合結果的精度影響較大[60]。

3 結論與展望

本文對目前基于InSAR技術的礦區大梯度形變監測方法進行總結，通過分析國內外研究現狀，對兩種方法所采用的模型公式進行了詳細分析，總結了目前各方法存在的不足，可為后續相關內容的研究提供參考。

當前,隨著越來越多雷達衛星的發射(Sentinel-1A/B、ALOS-2、TerraSAR-X等)，國內外學者可獲取影像的時間分辨率更高、空間分辨率更大，使得干涉測量的可檢測形變范圍增大，給礦區短時間大梯度的變形監測帶來幫助，可在一定程度上提高測量結果的準確性。此外，針對不同的礦山開采沉陷模型，依據不同地質環境因素及地表覆蓋物的影響，研究選取更加合理的開采沉陷模型，結合InSAR技術獲取的視線向形變值來完成礦區大梯度形變信息的提取，將是未來研究的主要方向。