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

陳 磊，趙學(xué)勝

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

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

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

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

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

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

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

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

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

(1)

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

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

Dmax=dmax+0.002(γ-1)

(2)

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

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

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

(3)

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

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

(4)

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

(5)

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

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

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

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

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

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

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

(6)

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

(7)

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

(8)

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

dLOS(x,y)=

(9)

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

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

(10)

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

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

(11)

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

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

(12)

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

(13)

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

3 結(jié)論與展望

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

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