結(jié)合InSAR的礦區(qū)大梯度變形監(jiān)測技術(shù)進展

陳 磊，趙學(xué)勝

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

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

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

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

1 結(jié)合像元追蹤法與InSAR技術(shù)

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

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

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

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

(1)

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

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

Dmax=dmax+0.002(γ-1)

(2)

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

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

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

(3)

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

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

(4)

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

(5)

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

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

2 結(jié)合礦山開采沉陷模型與InSAR技術(shù)

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

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

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

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

(6)

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

(7)

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

(8)

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

dLOS(x,y)=

(9)

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

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

(10)

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

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

(11)

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

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

(12)

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

(13)

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

3 結(jié)論與展望

本文對目前基于InSAR技術(shù)的礦區(qū)大梯度形變監(jiān)測方法進行總結(jié)，通過分析國內(nèi)外研究現(xiàn)狀，對兩種方法所采用的模型公式進行了詳細(xì)分析，總結(jié)了目前各方法存在的不足，可為后續(xù)相關(guān)內(nèi)容的研究提供參考。

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