許耀宗，李 濤，范洪冬

(1.中國礦業(yè)大學(xué) 自然資源部國土環(huán)境與災(zāi)害監(jiān)測(cè)重點(diǎn)實(shí)驗(yàn)室，江蘇 徐州 221116；2.山東省煤田地質(zhì)局物探測(cè)量隊(duì)，山東 濟(jì)南 250104)

0 引言

煤炭資源在我國能源結(jié)構(gòu)中占據(jù)主導(dǎo)地位，為國民經(jīng)濟(jì)快速增長提供重要力量。同時(shí)，因煤炭資源大規(guī)模開發(fā)利用，造成地表塌陷、山體滑坡、水土流失以及建構(gòu)筑物毀壞等，給生態(tài)環(huán)境和人民生產(chǎn)生活帶來負(fù)面影響。因此，對(duì)礦區(qū)開采造成的地表形變進(jìn)行高精度監(jiān)測(cè)具有重要的現(xiàn)實(shí)意義[1-2]。

差分合成孔徑雷達(dá)干涉測(cè)量(Different Interferometric Synthetic Aperture Radar,D-InSAR)技術(shù)理論上可以擺脫對(duì)地面特征點(diǎn)的依賴，獲取亞厘米級(jí)的地表形變，但該技術(shù)容易受時(shí)空失相干和大氣延遲等因素制約。因此，時(shí)序干涉SAR技術(shù)被提出，例如永久散射體干涉測(cè)量技術(shù)(PS-InSAR)[3]、斯坦福PS方法(StaMPS)[4]、小基線集(SBAS-InSAR)[5]、臨時(shí)相干點(diǎn)(TCP-InSAR)[6]，該技術(shù)克服傳統(tǒng)D-InSAR技術(shù)不足，在實(shí)現(xiàn)低成本、大規(guī)模、高精度地表沉降監(jiān)測(cè)方面具有較大優(yōu)勢(shì)。但時(shí)序InSAR方法在非城鎮(zhèn)、建構(gòu)筑物較少的野外，難以獲取所需監(jiān)測(cè)點(diǎn)密度，監(jiān)測(cè)結(jié)果不能完整反映整個(gè)形變場(chǎng)。針對(duì)這一問題，在常規(guī)時(shí)序InSAR基礎(chǔ)上，學(xué)者提出SqueeSAR[7]、CAEPSInSAR[8]等融合分布式散射目標(biāo)(Distributed Scatterer,DS)的時(shí)序地表形變監(jiān)測(cè)方法，顯著提高非城鎮(zhèn)地區(qū)監(jiān)測(cè)點(diǎn)密度，成為時(shí)序InSAR領(lǐng)域研究熱點(diǎn)之一[9-11]。

同質(zhì)點(diǎn)選取是DS-InSAR需要完成的首要任務(wù)，其結(jié)果直接影響相位優(yōu)化和最終地表形變精度。為提高DS-InSAR在礦區(qū)應(yīng)用中的精度、點(diǎn)位密度及處理效率，本文提出融合LRT的Bootstrap區(qū)間估計(jì)同質(zhì)點(diǎn)選取方法[12-17](LRT_Bootstrap)，并以鄂爾多斯某煤礦為例，驗(yàn)證DS-InSAR方法的可靠性。

1 DS-InSAR方法原理

(1)

當(dāng)δ的分布已知，可得式(2)：

P(δ1-α/2≤δ≤δα/2)=1-α

(2)

式中：P(*)表示發(fā)生*的概率；1-α為置信度；δ1-α/2為(1-α/2)的概率置信下限位置；δα/2為上α/2分位點(diǎn)。

結(jié)合式(1)～式(2)可得式(3)：

(3)

(4)

由Bootstrap原理以及大數(shù)定理可知，當(dāng)Bootstrap樣本個(gè)數(shù)足夠多時(shí)，可用δ*分布代替δ分布，由式(3)可得式(5)：

(5)

(6)

(7)

圖1 DS-InSAR數(shù)據(jù)處理流程Fig.1 Data processing flow chart of DS-InSAR

2 研究區(qū)域和數(shù)據(jù)

2.1 研究區(qū)域

東勝煤田某煤礦位于內(nèi)蒙古自治區(qū)鄂爾多斯市烏審旗與伊金霍洛旗之間，距伊金霍洛旗新街鎮(zhèn)較近，礦區(qū)內(nèi)人口稀少，以牧業(yè)為主，地形整體呈東高西低，最高點(diǎn)位海拔標(biāo)高為1 418 m，位于井田北部，最低點(diǎn)海拔標(biāo)高為1 316 m，位于井田西北部，最大地形高差為102 m，平均高差約15 m。井田南北方向平均寬約7.4 km，東西方向平均長約9.4 km，井田具有典型的高原堆積型沙丘地貌特征，地表大部被第四系風(fēng)積沙覆蓋，局部有白堊下統(tǒng)地層出露，植被稀疏，為沙漠-半沙漠地區(qū)。

2.2 數(shù)據(jù)集

本文選用研究區(qū)域2018年6月至2019年6月共31景Sentinel-1A衛(wèi)星影像數(shù)據(jù)，該影像數(shù)據(jù)波長約56 mm，方位向像元尺寸約13.9 m，斜距向像元尺寸約2.3 m，重訪周期為12 d，入射角約39°。由于獲取初始SAR影像覆蓋范圍較大，考慮計(jì)算機(jī)運(yùn)算時(shí)間以及研究區(qū)域范圍，最終將SAR影像進(jìn)行裁剪和過采樣，處理后的影像距離向和方位向尺寸均為600個(gè)×600個(gè)像元。

3 結(jié)果分析

3.1 模擬數(shù)據(jù)驗(yàn)證

本文采用蒙特卡羅隨機(jī)實(shí)驗(yàn)，檢驗(yàn)BWS、LRT和LRT_Bootstrap 3種同質(zhì)點(diǎn)選取方法。在復(fù)圓高斯分布假設(shè)前提下，設(shè)置格網(wǎng)尺寸大小為[15×15]，設(shè)置[15×8]格網(wǎng)中像元幅度真值為μ1，剩余[15×7]像元幅度真值為μ2；以格網(wǎng)中心位置(8,8)為參考像元點(diǎn)，剩余像元作為待探測(cè)像元；設(shè)置樣本N=30，分別從參數(shù)σ1，σ2的瑞利分布總體中模擬加噪聲幅度序列；用3種同質(zhì)點(diǎn)選取方法進(jìn)行實(shí)驗(yàn)，該過程重復(fù)5 000次，得到3種方法的拒絕率均值和拒絕率標(biāo)準(zhǔn)差如圖2所示。通過變換σ1/σ2比值，重復(fù)實(shí)驗(yàn)還可得到不同對(duì)比度下拒絕率均值和標(biāo)準(zhǔn)差[18]。

圖2 3種同質(zhì)點(diǎn)選取方法的拒絕率均值和拒絕率標(biāo)準(zhǔn)差Fig.2 Mean rejection rate and standard deviation of rejection rate by three homogeneous pixels selection methods

理論上，真實(shí)拒絕率為異質(zhì)像元個(gè)數(shù)與置信水平為α?xí)r同質(zhì)像元的第1類誤差像元個(gè)數(shù)之和，當(dāng)α取0.05時(shí)，真實(shí)拒絕率為49.33%。由圖2(a)可知，隨對(duì)比度增加，3種方法均收斂，其中，LRT_Bootstrap的拒絕率均值相對(duì)較高，收斂速度更快。在圖2(b)中，LRT_Bootstrap的拒絕率標(biāo)準(zhǔn)差相對(duì)較小且穩(wěn)定，當(dāng)σ1/σ2的比值小于1.5時(shí)，也可保證拒絕率標(biāo)準(zhǔn)差較小。同質(zhì)點(diǎn)選取主要用于相干性矩陣的精確估計(jì)，當(dāng)選取的獨(dú)立樣本個(gè)數(shù)大于30時(shí)，得到的估計(jì)結(jié)果近似無偏[21]。綜上，LRT_Bootstrap擁有較高拒絕率均值和較低且穩(wěn)定的拒絕率標(biāo)準(zhǔn)差，表明其適用于沒有明顯地物區(qū)分度的野外地區(qū)同質(zhì)點(diǎn)選取。

圖3 不同算法在研究區(qū)域內(nèi)同質(zhì)點(diǎn)的選取結(jié)果Fig.3 Selection results of homogeneous pixels in study area by different algorithms

3.2 真實(shí)數(shù)據(jù)同質(zhì)點(diǎn)選取結(jié)果

3種算法在研究區(qū)域內(nèi)同質(zhì)點(diǎn)選取結(jié)果如圖3所示，廠房、公路以及部分零星分布的建筑物、巖石等為強(qiáng)散射體，多屬于PS目標(biāo)，該地點(diǎn)的同質(zhì)點(diǎn)數(shù)量相對(duì)較少；部分地點(diǎn)地表覆蓋不均勻，區(qū)域同質(zhì)樣本數(shù)量較少。BWS檢驗(yàn)和LRT方法選取的同質(zhì)點(diǎn)結(jié)果整體顏色偏黃，2種方法像元同質(zhì)點(diǎn)個(gè)數(shù)均值分別為145，132個(gè)，LRT_Bootstrap為122個(gè)，證實(shí)BWS檢驗(yàn)和LRT方法拒絕率相對(duì)較低。BWS檢驗(yàn)無法完全排除來自不同分布的樣本，得到的同質(zhì)點(diǎn)集合中包含異質(zhì)樣本，并最終導(dǎo)致同質(zhì)樣本數(shù)增多[18]。LRT估計(jì)是無偏的，可有效控制第1類誤差，但當(dāng)像元點(diǎn)差異不明顯時(shí)，似然比檢驗(yàn)的取偽概率較大，這也是似然比檢驗(yàn)選擇的同質(zhì)樣本數(shù)量較多的原因。同質(zhì)點(diǎn)選取主要是進(jìn)行相干性矩陣的精確估計(jì)，當(dāng)選取的獨(dú)立樣本個(gè)數(shù)大于30時(shí)，得到的估計(jì)結(jié)果近似無偏，LRT_Bootstrap在實(shí)驗(yàn)地區(qū)選取的同質(zhì)點(diǎn)個(gè)數(shù)大于30的像元個(gè)數(shù)為334 713個(gè)，占總體像元的93%，確保大部分地區(qū)有足夠的同質(zhì)點(diǎn)滿足相位優(yōu)化需求。此外，除使用LRT結(jié)果作為初值外，還可使用非參數(shù)檢驗(yàn)結(jié)果作為初始值，但若初始值不準(zhǔn)確，會(huì)導(dǎo)致選取的同質(zhì)點(diǎn)出現(xiàn)問題，甚至部分區(qū)域選不上點(diǎn)[22]。

在計(jì)算效率方面，選用CPU為Intel(R) Core(TM) i7-8700K，運(yùn)行內(nèi)存64 G，采用Matlab程序進(jìn)行同質(zhì)點(diǎn)選取，BWS檢驗(yàn)選取同質(zhì)點(diǎn)樣本所用時(shí)間為1 413 s，似然比檢驗(yàn)耗時(shí)20 s，耗時(shí)相對(duì)最短，LRT_Bootstrap方法耗時(shí)適中，為98 s。

圖4 研究區(qū)域參考像元及其同質(zhì)樣本Fig.4 Reference pixel and its homogeneous samples of study area

利用BWS檢驗(yàn)、LRT、LRT_Bootstrap和目視解譯的實(shí)驗(yàn)區(qū)域同質(zhì)點(diǎn)選取情況對(duì)比，如圖4所示。其中，紅色點(diǎn)為參考點(diǎn)，同質(zhì)區(qū)域和異質(zhì)區(qū)域參考點(diǎn)在影像中的坐標(biāo)分別為(552,344)，(67,104)(以圖像左上角為原點(diǎn))，綠色點(diǎn)為選出來的與參考點(diǎn)類似的點(diǎn)(即同質(zhì)點(diǎn))，沒有選出來的點(diǎn)即與參考點(diǎn)不一致的點(diǎn)(即非同質(zhì)點(diǎn))。由圖4可知，在異質(zhì)區(qū)域，BWS檢驗(yàn)拒絕率相對(duì)最低，同質(zhì)點(diǎn)樣本相對(duì)最多，選擇的同質(zhì)點(diǎn)樣本中明顯包含非同質(zhì)點(diǎn)；LRT具有較多異質(zhì)點(diǎn)；LRT_Bootstrap方法的高拒絕率使其完全排除異質(zhì)點(diǎn)，雖然漏選7個(gè)同質(zhì)樣本點(diǎn)，但選取的同質(zhì)點(diǎn)數(shù)量達(dá)到32個(gè)。在同質(zhì)區(qū)域，BWS檢驗(yàn)得到的同質(zhì)點(diǎn)個(gè)數(shù)相對(duì)最多，其次為LRT，由于較高的拒絕率使得LRT_Bootstrap方法得到的同質(zhì)點(diǎn)個(gè)數(shù)略低，同質(zhì)樣本個(gè)數(shù)達(dá)196個(gè)。

3.3 時(shí)序監(jiān)測(cè)結(jié)果與驗(yàn)證

基于Sentinel-1A數(shù)據(jù)和圖1，利用LRT_Bootstrap方法進(jìn)行同質(zhì)樣本點(diǎn)選取以及相位優(yōu)化，求取研究區(qū)域地表時(shí)序沉降。利用DS-InSAR技術(shù)，得到煤礦2018年6月至2019年6月監(jiān)測(cè)時(shí)段內(nèi)豎直向累積沉降值如圖5所示，最大豎直向累積沉降值為-128.9 mm。采用基于LRT_Bootstrap的DS-InSAR進(jìn)行處理，最終獲取形變點(diǎn)數(shù)量為70 784個(gè)，約占影像總像元的20%，可較好地反映煤礦開采引起的地表沉降情況。

本文選取與水準(zhǔn)點(diǎn)重合的DS-InSAR監(jiān)測(cè)點(diǎn)累計(jì)沉降結(jié)果，與水準(zhǔn)觀測(cè)數(shù)據(jù)進(jìn)行比較如圖6所示，若水準(zhǔn)點(diǎn)上沒有與之重合的InSAR監(jiān)測(cè)點(diǎn)，則取1個(gè)像元內(nèi)距離水準(zhǔn)點(diǎn)最近的InSAR監(jiān)測(cè)點(diǎn)，若距離超過1個(gè)像元，則舍棄該水準(zhǔn)點(diǎn)。由圖6可知，實(shí)測(cè)水準(zhǔn)觀測(cè)結(jié)果與基于LRT_Bootstrap的DS-InSAR觀測(cè)結(jié)果具有高度一致性，最大差值為27.4 mm，均方根誤差和平均絕對(duì)誤差分別為14.4，12.0 mm，一定程度上驗(yàn)證監(jiān)測(cè)結(jié)果的可靠性。

圖6 基于LRT_Bootstrap的DS-InSAR監(jiān)測(cè)結(jié)果與實(shí)測(cè)水準(zhǔn)數(shù)據(jù)對(duì)比Fig.6 Comparison of DS-InSAR monitoring results and measured level data

為進(jìn)一步驗(yàn)證基于LRT_Bootstrap的DS-InSAR形變結(jié)果有效性，本文選取4種方法形變結(jié)果的16個(gè)點(diǎn)進(jìn)行比較，如圖7所示，4種方法形變結(jié)果具有高度一致性，基于LRT_Bootstrap的形變結(jié)果與實(shí)測(cè)水準(zhǔn)數(shù)據(jù)的平均絕對(duì)誤差為11.5 mm，均方根誤差為12.9 mm；基于LRT和BWS的形變結(jié)果與實(shí)測(cè)水準(zhǔn)數(shù)據(jù)的平均絕對(duì)誤差分別為12.2，13.0 mm，均方根誤差分別為14.5，14.5 mm，形變精度略低于基于LRT_Bootstrap的DS-InSAR。

圖7 基于BWS、LRT和LRT_Bootstrap的DS-InSAR與水準(zhǔn)形變結(jié)果對(duì)比Fig.7 Comparison of deformation results of DS-InSAR based on BWS,LRT,LRT_Bootstrap and level

4 結(jié)論

1)對(duì)比3種同質(zhì)點(diǎn)選取算法，發(fā)現(xiàn)LRT_Bootstrap更適合于區(qū)分度不明顯的野外地區(qū)的同質(zhì)點(diǎn)探測(cè)。

2)基于LRT_Bootstrap的DS-InSAR監(jiān)測(cè)結(jié)果和水準(zhǔn)數(shù)據(jù)的監(jiān)測(cè)結(jié)果具有高度一致性，并且基于LRT_Bootstrap的DS-InSAR形變結(jié)果精度更高。

3)利用DS-InSAR方法對(duì)鄂爾多斯某煤礦開采過程的沉降情況進(jìn)行成功監(jiān)測(cè)，DS-InSAR監(jiān)測(cè)結(jié)果可以反映實(shí)際自然地表下煤礦沉降位置、范圍以及大小，可為煤礦開采沉陷邊界確定、沉降區(qū)域地質(zhì)災(zāi)害預(yù)測(cè)以及地表動(dòng)態(tài)沉陷規(guī)律獲取提供參考。