基于SAR 技術(shù)的開采沉陷全盆地分區(qū)形變提取方法

鄒英杰，范洪冬，孫 葉，許耀宗

（中國礦業(yè)大學(xué) 自然資源部國土環(huán)境與災(zāi)害監(jiān)測重點實驗室，江蘇 徐州 221116）

煤炭工業(yè)作為我國重要的基礎(chǔ)工業(yè)，在帶來巨大經(jīng)濟(jì)效益的同時，地下煤層開采活動也會引發(fā)地表人工建構(gòu)筑物損壞和一系列的生態(tài)環(huán)境問題[1-2]。因此，對礦區(qū)地表進(jìn)行精細(xì)化形變監(jiān)測十分必要。常規(guī)地表形變監(jiān)測手段，如水準(zhǔn)測量和GNSS 等方法觀測精度雖高，但存在工作量大、觀測點位密度低且易被破壞等缺點。

合成孔徑雷達(dá)差分干涉測量（DInSAR）技術(shù)具有全天時、全天候、空間分辨率高、覆蓋區(qū)域廣等優(yōu)點[3]，其利用雷達(dá)視線向（LOS）相位信息能夠監(jiān)測到地表厘米級甚至亞厘米級的形變[4]；但是該技術(shù)受失相干影響較為嚴(yán)重，且受限于相位解纏的理論形變梯度閾值[5]，只能監(jiān)測地表小量級形變。偏移跟蹤（OT）方法基于SAR 影像的強(qiáng)度而非相位信息，因此其不受時空失相干和形變梯度的限制，多用于地震、冰川、滑坡、礦區(qū)等大量級形變監(jiān)測領(lǐng)域[6-10]；OT方法的監(jiān)測精度會受到匹配窗口大小的影響[11]，實現(xiàn)匹配窗口的自適應(yīng)選取可有效提升OT 方法的監(jiān)測精度。

近年來，DInSAR 和OT 相結(jié)合的方法在礦區(qū)地表形變監(jiān)測方面已經(jīng)取得了相應(yīng)成果，通過設(shè)定相干系數(shù)閾值或掩膜大形變區(qū)是將地表形變分為大小2 個量級的主流方法[12-16]。在此基礎(chǔ)上，相關(guān)學(xué)者根據(jù)先驗知識和目視解譯等方法劃分出中等量級形變區(qū)域[11,17-19]，并對該區(qū)域進(jìn)行不同的融合處理。上述劃分形變區(qū)域方法的適用性和范圍會受到現(xiàn)場實測數(shù)據(jù)和人為主觀因素的影響，所采用的OT 方法均未考慮到匹配窗口的自適應(yīng)選取。因此，通過相干系數(shù)和地表形變梯度劃分出中等量級形變區(qū)域，聯(lián)合該區(qū)域累積DInSAR 方法和自適應(yīng)OT 方法的監(jiān)測結(jié)果進(jìn)行加權(quán)融合，實現(xiàn)對礦區(qū)沉陷盆地大、中、小不同量級形變的分級與精細(xì)監(jiān)測。

1 開采沉陷盆地分區(qū)形變提取方法

1.1 累積DInSAR 方法原理

SAR 衛(wèi)星通過向地面發(fā)射電磁波信號，接收機(jī)接收到地面的反射信號后，經(jīng)處理得到SAR 單視復(fù)數(shù)影像，其中包含地物的輻射強(qiáng)度信息和相位信息。DInSAR 技術(shù)利用形變前后觀測的相位進(jìn)行干涉測量，其沿LOS 向路徑差△r 可以用干涉相位的相位差表示，即：

式中：Φ 為干涉相位的相位差；Φflat、Φtopo分別為平地相位和地形相位，可以利用SAR 衛(wèi)星成像參數(shù)和外部高精度DEM 進(jìn)行去除；Φdef為地表形變引起的形變相位；Φnoise為其他因素引起的噪聲相位。

SAR 衛(wèi)星可以以一定的周期重復(fù)獲取同一地區(qū)的時間序列影像，隨著SAR 數(shù)據(jù)的積累，一系列時序InSAR 方法得到快速發(fā)展，如SBAS-InSAR、PS-InSAR、StaMPS-InSAR 等。累積DInSAR 方法作為1 種最簡單的時序InSAR 技術(shù)，對一系列SAR 數(shù)據(jù)進(jìn)行差分干涉處理，將差分處理后的解纏相位進(jìn)行累積，表達(dá)式如下：

式中：Φsum為累積后的形變相位；Φiumw為第i 個差分干涉處理后的解纏相位；n 為SAR 序列影像的數(shù)量。

最后，可根據(jù)式（2）計算出成像期間地面點沿雷達(dá)LOS 向的累積形變量。

1.2 自適應(yīng)OT 方法原理

傳統(tǒng)OT 方法利用光學(xué)影像匹配中的互相關(guān)系匹配算法（NCC），對主影像和從影像按照一定窗口進(jìn)行精確匹配，通過式（3）計算互相關(guān)系數(shù)和尋找互相關(guān)系數(shù)矩陣峰值的位置來確定2 幅影像在方位向和距離向上的像素級偏移量。為達(dá)到亞像素級的配準(zhǔn)精度，可通過插值、曲面擬合等方法對相關(guān)系數(shù)進(jìn)行擬合，將極值點的位置作為匹配的最佳位置。

式中：ρ（x，y）為像元點（x，y）的互相關(guān)系數(shù)；I1、I2分別為模板窗口和搜索窗口內(nèi)像素點的像元強(qiáng)度；Iˉ1、Iˉ2分別為模板窗口和搜索窗口內(nèi)所有像素點的像元強(qiáng)度平均值；（u，v）為像素點在從影像中相對于主影像的距離向和方位向偏移量。

自適應(yīng)OT 算法旨在實現(xiàn)匹配窗口的自適應(yīng)選取，提高OT 方法的監(jiān)測精度。因此，需要尋找1 個最合適的匹配窗口，為了衡量不同窗口間的匹配結(jié)果的可信度，采用信噪比SNR 作為衡量標(biāo)準(zhǔn)，其定義如式（4）：

式中：ρmax為相關(guān)系數(shù)峰值；ρˉ為3×3 窗口內(nèi)相關(guān)系數(shù)平均值。

SNR 值越大，代表相干系數(shù)峰值所在位置越突出，越有利于進(jìn)行圖像的精確匹配，將SNR 最大值的窗口大小作為OT 方法匹配的最佳窗口。自適應(yīng)OT 方法通過尋找以初始步長為b（通常為8）的一定窗口區(qū)間［a，d］內(nèi)的SNR 最大值的窗口大小a+bi，并繼續(xù)在［a+b（i-1）a+b（i+1）］窗口區(qū)間內(nèi)以步長為b/2 尋找新的SNR 最大值的窗口大小，以此類推，直到找到最終的最佳匹配窗口；在窗口大小自適應(yīng)選取的基礎(chǔ)上還需完成OT 方法后續(xù)相關(guān)系數(shù)插值以及距離向與方位向亞像素級偏移量解算等步驟。自適應(yīng)OT 方法的實現(xiàn)過程如圖1，

圖1 自適應(yīng)OT 方法流程圖Fig.1 Processing flow of adaptive OT

1.3 累積DInSAR 和自適應(yīng)OT 聯(lián)合加權(quán)方法原理

地下煤層開采過程中，上覆地表形變量級通常是從下沉中心到邊界逐漸減小的，從前面的分析可以看出，累積DInSAR 方法可以測量下沉盆地邊緣小尺度地表形變，自適應(yīng)OT 方法可以測量下沉盆地中心部分大尺度地表形變。對累積DInSAR 和自適應(yīng)OT 方法監(jiān)測精度較低的中等量級形變可以利用二者得到的結(jié)果進(jìn)行加權(quán)融合，實現(xiàn)對礦區(qū)多量級地表形變精細(xì)監(jiān)測。

對于小量級形變區(qū)域，采用累積DInSAR 方法，由于相干系數(shù)越高代表地表形變的解算結(jié)果越可靠，多數(shù)文獻(xiàn)及GAMMA 軟件說明書在進(jìn)行差分干涉處理時，一般取相干系數(shù)0.3 作為執(zhí)行相位解纏的閾值[13，20]。對于不同分辨率SAR 影像的處理閾值并無明確說明，考慮到大柳塔礦52303 工作面地表形變大而導(dǎo)致利用相位解算地表形變誤差大，在進(jìn)行試驗時，曾選取多個相干性閾值進(jìn)行試驗，當(dāng)閾值為0.45 時，所劃分的小量級形變區(qū)域效果最佳。

由于相鄰像元形變梯度大于1/4 個波長的區(qū)域超出了單次相位解纏的極限要求，而礦區(qū)因開采導(dǎo)致的地表沉降超過2 m 的情況非常普遍，若要保證本研究累積DInSAR 方法能夠參與地表形變解算，只能在礦區(qū)地表形變梯度小于（n-1）/4 個波長的區(qū)域。概率積分法作為礦區(qū)開采沉陷適用性最廣的方法[21-22]，本研究根據(jù)地質(zhì)采礦參數(shù)和概率積分法預(yù)計并預(yù)計出相鄰像元梯度大于（n-1）/4 個波長的區(qū)域作為大量級形變區(qū)域，該區(qū)域采用自適應(yīng)OT 方法進(jìn)行處理。

對于其它區(qū)域則作為中等量級形變采用累積DInSAR 與自適應(yīng)OT 方法按照一定的權(quán)值進(jìn)行融合處理，表達(dá)式如下：

式中：DLOS為本方法最后得到的地表LOS 向形變信息；P1、P2分別為DInSAR 方法和自適應(yīng)OT 方法所占權(quán)值，一般依據(jù)先驗定權(quán)的方式或經(jīng)驗確定；C 為相干系數(shù)；Div 為概率積分法預(yù)計LOS 向形變梯度；ΣDDInSAR為累積DInSAR 監(jiān)測結(jié)果；DOT為自適應(yīng)OT 方法監(jiān)測結(jié)果；n 為SAR 影像數(shù)量；λ 為雷達(dá)波長。

根據(jù)現(xiàn)場實測資料先驗定權(quán)的方式，其權(quán)值確定的具體公式如下：

式中：σΣDInSAR、σOT分別為累積DInSAR 方法和自適應(yīng)OT 方法相對于實測水準(zhǔn)數(shù)據(jù)或GPS 數(shù)據(jù)的均方根誤差，并將其作為監(jiān)測精度。

2 試驗數(shù)據(jù)概況

2.1 礦區(qū)概況

試驗選取榆林大柳塔煤礦52303 工作面為研究區(qū)域。大柳塔煤礦地處陜西省榆林市神木市與內(nèi)蒙古鄂爾多斯市伊金霍洛旗交界處，研究區(qū)地理位置及GPS 點位信息如圖2，紅色的圓點為地表移動觀測GPS 點。該地平均海拔在1 250 m，屬于北溫帶大陸性半干旱氣候，冬季寒冷，夏季干熱，春季多風(fēng)，秋季涼爽。52303 工作面于2012 年12 月開始試采，工作面長度約301.5 m，推進(jìn)長度約為4 443 m，平均采深250 m，平均采高6.6 m，煤層傾角1°～3°，采用走向綜合機(jī)械化長壁采煤法，頂板管理方法為全部垮落法。作為1 座年產(chǎn)量超2 000 萬t 的超大型現(xiàn)代化煤礦，長久的地下采煤活動造成地表不同程度的地表變形，危及礦區(qū)生產(chǎn)安全和周圍人民生命財產(chǎn)安全。因此，對該礦區(qū)進(jìn)行大范圍、多量級的地表形變精細(xì)監(jiān)測具有重要的意義。

圖2 研究區(qū)地理位置及GPS 點位信息Fig.2 Location of the study area and the information of GPS points

2.2 SAR 數(shù)據(jù)概況

研究采用德國航天中心發(fā)射的TerraSAR-X 衛(wèi)星高空間分辨率雷達(dá)影像，該衛(wèi)星波長3.1 cm，衛(wèi)星重訪周期為11 d，其距離向和方位向空間分辨率分別為0.91 m 和1.82 m。為了便于試驗精度評定，選取TerraSAR-X 衛(wèi)星成像日期為2014-04-05—2014-05-08 日的4 景影像，現(xiàn)場實測GPS 數(shù)據(jù)日期為2014-04-06—2014-05-08 日，可認(rèn)為SAR 衛(wèi)星和GPS 二者為同期觀測。所選取的SAR 影像對數(shù)據(jù)具體參數(shù)見表1。

表1 TerraSAR 影像對參數(shù)Table 1 TerraSAR image pairs parameters

3 試驗結(jié)果及分析

3.1 數(shù)據(jù)處理結(jié)果

累積DInSAR 試驗使用4 景時間序列SAR 影像做DInSAR 處理，在DInSAR 處理過程中引入外部STRM 分辨率為30 m 的DEM 數(shù)據(jù)去除地形相位。為了保證在進(jìn)行差分干涉處理過程中起始解纏相位不產(chǎn)生相對變化，在每次DInSAR 處理過程中選擇同1 個遠(yuǎn)離形變區(qū)域、保持高相干性的相位解纏起始點。將解纏后的形變相位進(jìn)行時序疊加后，根據(jù)式（2）轉(zhuǎn)換為LOS 向形變，累積DInSAR 方法得到的礦區(qū)沉陷盆地邊緣小形變區(qū)域監(jiān)測結(jié)果如圖3。

圖3 累積DInSAR 方法形變監(jiān)測結(jié)果Fig.3 Deformation monitoring results of cumulative DInSAR

利用2014-04-05 和2014-05-08 日2 景影像進(jìn)行自適應(yīng)OT 處理，本研究自適應(yīng)OT 方法設(shè)置初始窗口區(qū)間為［73-121］，初始步長b 為8。同時，根據(jù)現(xiàn)場實測資料以及礦區(qū)采煤地質(zhì)條件，判斷出該工作面地表LOS 向形變值應(yīng)處于-0.5～4.5 m 的范圍內(nèi)，因此需要對自適應(yīng)OT 方法監(jiān)測結(jié)果內(nèi)大于4.5 m 和小于-0.5 m 的形變值進(jìn)行濾除處理。研究使用均值濾波對濾除后空洞進(jìn)行填充，得到自適應(yīng)OT方法在礦區(qū)沉陷盆地中心大量級形變亞像素級監(jiān)測結(jié)果如圖4。

圖4 自適應(yīng)OT 方法形變監(jiān)測結(jié)果Fig.4 Deformation monitoring results of adaptive OT

研究提出的加權(quán)融合方法以相干系數(shù)0.45 和預(yù)計形變梯度0.023 為閾值，選擇出中等量級形變區(qū)域。結(jié)合該區(qū)域現(xiàn)場實測GPS 數(shù)據(jù)分別對自適應(yīng)OT 方法和累積DInSAR 方法進(jìn)行精度驗證，累積DInSAR 方法的監(jiān)測精度為8.5 cm，自適應(yīng)OT 方法的監(jiān)測精度為8.2 cm，根據(jù)式（6）進(jìn)行定權(quán)，確定中等量級形變區(qū)域累積DInSAR 方法和自適應(yīng)OT 方法所占權(quán)值分別為0.48 和0.52。根據(jù)式（5）對累積DInSAR 方法和自適應(yīng)OT 方法的監(jiān)測結(jié)果進(jìn)行加權(quán)融合，得到中等量級形變監(jiān)測結(jié)果。

對3 種量級形變進(jìn)行線性疊加后，得到礦區(qū)完整多量級地表形變場，其監(jiān)測結(jié)果如圖5。

圖5 本研究形變監(jiān)測結(jié)果Fig.5 Deformation monitoring results of the method in this paper

3.2 試驗結(jié)果分析

結(jié)合52303 工作面開采工作進(jìn)度和現(xiàn)場地表實測數(shù)據(jù)，此時52303 工作面開采工作已經(jīng)進(jìn)入最后階段直至停采了一段時間，因此沉陷盆地形變中心靠近5-2煤輔運大巷?？紤]到礦區(qū)開采過程中，地表移動具有一定的滯后性，此期間地表LOS 向形變最大值已超過3 m，所生成的地表沉陷盆地仍可以用本研究方法進(jìn)行形變監(jiān)測。通過圖3 和圖4 可以直觀發(fā)現(xiàn)，累積DInSAR 所監(jiān)測到的最大形變小于0.2 m，遠(yuǎn)小于礦區(qū)實際地表形變。自適應(yīng)OT 方法可以彌補(bǔ)DInSAR 方法無法監(jiān)測大梯度沉降礦區(qū)的缺陷，但是通過圖4 可以看出，自適應(yīng)OT 方法在礦區(qū)周邊非形變區(qū)域出現(xiàn)了錯誤的監(jiān)測結(jié)果，說明自適應(yīng)OT 方法在小量級形變區(qū)域解算精度并不高。對于礦區(qū)大量級沉陷盆地，DInSAR 和OT 方法進(jìn)行聯(lián)合監(jiān)測是非常有必要的，本研究所提出的融合累積DInSAR 和自適應(yīng)OT 的方法可以有效地提取礦區(qū)沉陷全盆地多量級形變場。利用礦區(qū)52303 工作面傾向19 個地表GPS 觀測點對本研究提出方法所監(jiān)測的多量級形變結(jié)果進(jìn)行精度分析，將二者結(jié)果進(jìn)行對比，19 個地表觀測站形變對比如圖6。

圖6 地表觀測站LOS 向形變對比圖Fig.6 Comparison of LOS deformation at the surface observation station

在中等量級形變區(qū)域，自適應(yīng)OT 方法、累積DInSAR 方法和融合方法的監(jiān)測精度對比見表2。對研究區(qū)域內(nèi)3 種不同量級的形變監(jiān)測結(jié)果進(jìn)行精度分析，精度對比見表3。

表2 中等量級形變區(qū)方法精度分析表Table 2 Accuracy analysis table of medium order deformation region method

表3 研究區(qū)不同量級形變精度分析表Table 3 Analysis table of deformation accuracy of different orders in the study area

由表2 可知：采用累積DInSAR 方法，其均方根誤差為0.085 m，平均絕對偏差為0.065 m，自適應(yīng)OT 方法均方根誤差為0.082 m，平均絕對誤差為0.058 m。根據(jù)權(quán)值對2 種方法的監(jiān)測結(jié)果進(jìn)行融合后，將融合加權(quán)結(jié)果與現(xiàn)場實測結(jié)果相比，加權(quán)融合方法均方根誤差為0.046 m，平均絕對誤差為0.036 m，其監(jiān)測精度明顯高于單獨使用累積DInSAR 或自適應(yīng)OT 方法，證明本研究所提出的加權(quán)融合方法能夠有效提高礦區(qū)沉陷盆地邊緣中小量級形變監(jiān)測的精度。

由表3 得到本方法提取的GPS 觀測點形變值的均方根誤差為0.108 m，最大絕對偏差值點號為Q14，為0.319 m。其中，自適應(yīng)OT 方法均方根誤差為0.154 m，平均絕對誤差為0.124 m；累積DInSAR方法均方根誤差為0.048 m，平均絕對誤差為0.037 m。自適應(yīng)OT 方法均方根誤差最大，這是因為自適應(yīng)OT 監(jiān)測形變精度受到像元尺寸的影響，其均方根誤差與距離向像素尺寸比值為16.48%。采用融合方法獲得的監(jiān)測精度相比于累積DInSAR 方法并沒有顯著提高，但是其在礦區(qū)地表下沉盆地邊緣中小量級形變的監(jiān)測范圍要優(yōu)于累積DInSAR 方法。

4 結(jié) 語

針對DInSAR 方法和OT 方法在中等量級形變監(jiān)測精度較低的問題，提出1 種融合累積DInSAR和自適應(yīng)OT 方法，實現(xiàn)對礦區(qū)多量級地表形變精細(xì)監(jiān)測。利用4 景TerraSAR 影像分別采用累積DInSAR、自適應(yīng)OT 以及二者加權(quán)融合方法對大柳塔煤礦52303 工作面開采沉陷盆地大量級、小量級和中等量級形變進(jìn)行監(jiān)測，得到礦區(qū)沉陷全盆地形變信息，與同期現(xiàn)場實測GPS 數(shù)據(jù)對比評價。

1）提出的融合方法監(jiān)測結(jié)果與19 個GPS 點實測數(shù)據(jù)對比，證明了融合累積DInSAR 和自適應(yīng)OT 的方法能夠獲取礦區(qū)完整多量級地表形變的有效性和可靠性，其監(jiān)測精度優(yōu)于0.108 m。

2）針對礦區(qū)地表形變盆地邊緣中小量級形變，累積DInSAR 自適應(yīng)OT 方法加權(quán)融合方法的均方根誤差優(yōu)于0.046 m，其監(jiān)測精度明顯高于單獨使用累積DInSAR 或自適應(yīng)OT 方法，相較于二者均方根誤差，分別提高了46%和44%。

3）累積DInSAR 和OT 融合的方法能夠得到礦區(qū)地表完整的形變場，但是該方法監(jiān)測精度受到SAR影像分辨率和雷達(dá)波長的影響，待未來長波長、高空間分辨率SAR 衛(wèi)星投入使用后，該方法的監(jiān)測精度及中小層次形變監(jiān)測范圍將會進(jìn)一步提高。