融合LRT的Bootstrap區間估計同質點識別方法研究*

許耀宗，李 濤，范洪冬

(1.中國礦業大學 自然資源部國土環境與災害監測重點實驗室，江蘇 徐州 221116；2.山東省煤田地質局物探測量隊，山東 濟南 250104)

0 引言

煤炭資源在我國能源結構中占據主導地位，為國民經濟快速增長提供重要力量。同時，因煤炭資源大規模開發利用，造成地表塌陷、山體滑坡、水土流失以及建構筑物毀壞等，給生態環境和人民生產生活帶來負面影響。因此，對礦區開采造成的地表形變進行高精度監測具有重要的現實意義[1-2]。

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

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

1 DS-InSAR方法原理

(1)

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

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

(2)

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

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

(3)

(4)

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

(5)

(6)

(7)

圖1 DS-InSAR數據處理流程Fig.1 Data processing flow chart of DS-InSAR

2 研究區域和數據

2.1 研究區域

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

2.2 數據集

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

3 結果分析

3.1 模擬數據驗證

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

圖2 3種同質點選取方法的拒絕率均值和拒絕率標準差Fig.2 Mean rejection rate and standard deviation of rejection rate by three homogeneous pixels selection methods

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

圖3 不同算法在研究區域內同質點的選取結果Fig.3 Selection results of homogeneous pixels in study area by different algorithms

3.2 真實數據同質點選取結果

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

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

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

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

3.3 時序監測結果與驗證

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

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

圖6 基于LRT_Bootstrap的DS-InSAR監測結果與實測水準數據對比Fig.6 Comparison of DS-InSAR monitoring results and measured level data

為進一步驗證基于LRT_Bootstrap的DS-InSAR形變結果有效性，本文選取4種方法形變結果的16個點進行比較，如圖7所示，4種方法形變結果具有高度一致性，基于LRT_Bootstrap的形變結果與實測水準數據的平均絕對誤差為11.5 mm，均方根誤差為12.9 mm；基于LRT和BWS的形變結果與實測水準數據的平均絕對誤差分別為12.2，13.0 mm，均方根誤差分別為14.5，14.5 mm，形變精度略低于基于LRT_Bootstrap的DS-InSAR。

圖7 基于BWS、LRT和LRT_Bootstrap的DS-InSAR與水準形變結果對比Fig.7 Comparison of deformation results of DS-InSAR based on BWS,LRT,LRT_Bootstrap and level

4 結論

1)對比3種同質點選取算法，發現LRT_Bootstrap更適合于區分度不明顯的野外地區的同質點探測。

2)基于LRT_Bootstrap的DS-InSAR監測結果和水準數據的監測結果具有高度一致性，并且基于LRT_Bootstrap的DS-InSAR形變結果精度更高。

3)利用DS-InSAR方法對鄂爾多斯某煤礦開采過程的沉降情況進行成功監測，DS-InSAR監測結果可以反映實際自然地表下煤礦沉降位置、范圍以及大小，可為煤礦開采沉陷邊界確定、沉降區域地質災害預測以及地表動態沉陷規律獲取提供參考。