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

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

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

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

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

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

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

1.1 累積DInSAR 方法原理

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

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

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

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

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

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

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

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

自適應(yīng)OT 算法旨在實(shí)現(xiàn)匹配窗口的自適應(yīng)選取，提高OT 方法的監(jiān)測(cè)精度。因此，需要尋找1 個(gè)最合適的匹配窗口，為了衡量不同窗口間的匹配結(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 方法通過(guò)尋找以初始步長(zhǎng)為b（通常為8）的一定窗口區(qū)間［a，d］內(nèi)的SNR 最大值的窗口大小a+bi，并繼續(xù)在［a+b（i-1）a+b（i+1）］窗口區(qū)間內(nèi)以步長(zhǎng)為b/2 尋找新的SNR 最大值的窗口大小，以此類推，直到找到最終的最佳匹配窗口；在窗口大小自適應(yīng)選取的基礎(chǔ)上還需完成OT 方法后續(xù)相關(guān)系數(shù)插值以及距離向與方位向亞像素級(jí)偏移量解算等步驟。自適應(yīng)OT 方法的實(shí)現(xiàn)過(guò)程如圖1，

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

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

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

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

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

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

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

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

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

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

2.1 礦區(qū)概況

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

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

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

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

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

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

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

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

圖3 累積DInSAR 方法形變監(jiān)測(cè)結(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］，初始步長(zhǎng)b 為8。同時(shí)，根據(jù)現(xiàn)場(chǎng)實(shí)測(cè)資料以及礦區(qū)采煤地質(zhì)條件，判斷出該工作面地表LOS 向形變值應(yīng)處于-0.5～4.5 m 的范圍內(nèi)，因此需要對(duì)自適應(yīng)OT 方法監(jiān)測(cè)結(jié)果內(nèi)大于4.5 m 和小于-0.5 m 的形變值進(jìn)行濾除處理。研究使用均值濾波對(duì)濾除后空洞進(jìn)行填充，得到自適應(yīng)OT方法在礦區(qū)沉陷盆地中心大量級(jí)形變亞像素級(jí)監(jiān)測(cè)結(jié)果如圖4。

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

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

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

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

3.2 試驗(yàn)結(jié)果分析

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

圖6 地表觀測(cè)站LOS 向形變對(duì)比圖Fig.6 Comparison of LOS deformation at the surface observation station

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

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

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

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

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

4 結(jié) 語(yǔ)

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

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

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

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