礦山開采后地表沉降邊界的D-InSAR 監(jiān)測與分析*

趙 恒，李亞星，楊可明*，竇 智

（1.淮北雙龍礦業(yè)有限責(zé)任公司地測科，安徽 淮北235000；2.中國礦業(yè)大學(xué)（北京）地球科學(xué)與測繪工程學(xué)院，北京100083）

地下煤層開采會對其上覆的巖層結(jié)構(gòu)產(chǎn)生破壞，進(jìn)而引起地面塌陷與沉降變形，煤礦區(qū)的采動嚴(yán)重影響威脅了采區(qū)地表的土地、水體、建（構(gòu)）筑物等的安全。為了保障煤礦區(qū)生產(chǎn)和生活安全的同時盡可能多回收煤炭資源，需要科學(xué)地確定開采沉降范圍，以便明確資源開采位置和界定相關(guān)責(zé)任，降低煤礦企業(yè)的開采風(fēng)險。

礦山開采引發(fā)地表沉降影響范圍的常規(guī)分析主要是依據(jù)開采工作面的實(shí)際地質(zhì)采礦條件、地表巖移觀測數(shù)據(jù)以及巖移規(guī)律等資料，運(yùn)用實(shí)測數(shù)據(jù)分析、巖層移動角幾何圖解、概率積分法（Probability Integral Method，PIM）預(yù)計和數(shù)值模擬計算等方法[1-4]，但這些常規(guī)方法存在較多且明顯的不足，如：巖移觀測站實(shí)測數(shù)據(jù)分析與幾何圖解主要體現(xiàn)在開采工作面的走向與傾向方向上，觀測點(diǎn)數(shù)量與圖解范圍相當(dāng)有限；PIM 與數(shù)值模擬方法雖然能大范圍確定地表沉降邊界，但不能反映不同的開采工作面存在不同的上覆巖層巖移與物理參數(shù)而引發(fā)實(shí)際沉降邊界的不同變化，導(dǎo)致所確定的邊界與實(shí)際位置存在一定甚至較大的偏差。而微波遙感的合成孔徑雷達(dá)干涉測量（Interferometric Synthetic Aperture Radar，InSAR）技術(shù)并不需要礦山開采的地質(zhì)采礦資料，而是基于側(cè)視雷達(dá)電磁波的地面反射回波相位干涉處理技術(shù)來獲取地表的沉降信息，常用技術(shù)有1989 年Gabriel 等首次提出的差分干涉測量（Differential InSAR，D-InSAR）[5]、2000 年Ferretti正式提出的永久散射體干涉測量（Permanent Scatterer InSAR，PS-InSAR）[6-7]和2002 年Berardino 提出的小基線集干涉測量（Small Baseline Subset InSAR，SBAS-InSAR）[8]等技術(shù)。InSAR 技術(shù)由于可以克服因地下開采條件復(fù)雜、地表移動變形，以及易受地下開采、覆巖巖性及結(jié)構(gòu)、地下水等多種因素影響而難以準(zhǔn)確界定地下開采所形成影響邊界的局限，所以現(xiàn)今已成為礦山開采形變監(jiān)測[9-10]、地質(zhì)災(zāi)害分析與環(huán)境影響評價[11-13]等開采沉陷研究的前沿技術(shù)。

D-InSAR 是目前較為先進(jìn)的一種地表形變監(jiān)測方法，與常規(guī)方法相比在監(jiān)測地表微小形變方面有著顯著的優(yōu)勢，它不僅能以厘米量級甚至毫米量級監(jiān)測地表形變，且能對大面積測區(qū)進(jìn)行快速監(jiān)測，從而得到采區(qū)主要形變區(qū)域與范圍[14-15]。本文將結(jié)合安徽省淮北礦區(qū)某煤礦井下311 工作面煤層開采下保護(hù)地面重要建（構(gòu)）筑物運(yùn)營安全等實(shí)際需求，利用D-InSAR 干涉測量技術(shù)監(jiān)測煤層工作面采動影響的地表沉降邊界，以下沉10mm 的位置為確定邊界的依據(jù)，并基于現(xiàn)場實(shí)測數(shù)據(jù)、SBAS-InSAR 技術(shù)分析結(jié)果和概率積分法變形預(yù)測等對所監(jiān)測到的地表沉降邊界進(jìn)行驗(yàn)證，進(jìn)而能為煤礦安全開采和最大程度回收地下煤炭資源提供科學(xué)的理論依據(jù)與技術(shù)保障。

1 理論與方法

1.1 D-InSAR 技術(shù)原理

D-InSAR 技術(shù)以其大范圍、低成本、速度快、高精度等優(yōu)勢，已被廣泛應(yīng)用于各種地表形變監(jiān)測中。如圖1 所示，D-InSAR 的基本思想從兩幅SAR 影像的干涉相位中分離出地表形變所引起的相位分量，從而對影像區(qū)域的地表沿雷達(dá)視線向地表形變進(jìn)行制圖。干涉相位是兩幅配準(zhǔn)過后的SAR 影像對應(yīng)像元的相位差，其組成可以表示如下[16]：

圖1 平地-高程相位函數(shù)模型

式中，φint表示干涉相位；φtop、φflat、φorb、φatm和φnoise分別表示高程相位、平地相位、軌道誤差相位、大氣相位和相位噪聲；φdef表示形變相位。軌道誤差可以通過精密定軌和模型擬合去除，噪聲一般通過濾波來抑制。礦區(qū)沉降空間范圍較小，大氣相位常常可以忽略不計。高程相位和平地相位可通過函數(shù)模型來去除。

圖1 中S1、S2分別表示主輔影像成像時的天線位置。B 為空間基線，B||和B⊥分別表示平行和垂直于雷達(dá)入射方向的基線分量，稱為平行基線和垂直基線。α 為空間基線與水平方向的夾角，θ 為雷達(dá)入射角，Δθ 為高程引起的雷達(dá)入射角變化。R1、R2表示天線與目標(biāo)點(diǎn)P 之間的斜距。R 表示P 點(diǎn)等斜距投影到參考橢球面的R1距離（R=R1）。h 為目標(biāo)點(diǎn)P 的高程。根據(jù)圖中幾何關(guān)系和泰勒公式，有：

式中，λ 為波長。因此可以得出平地相位和高程相位的表達(dá)式：

形變相位φdef與地表視線向形變dlos之間的關(guān)系如下：

1.2 地表沉降邊界分析方法

D-InSAR 是在InSAR 的基礎(chǔ)上發(fā)展起來的監(jiān)測地表微小形變的技術(shù)，在監(jiān)測工礦區(qū)地表形變方面具有連續(xù)空間覆蓋，高度自動化和高精度的優(yōu)勢。一般常用二軌法D-InSAR 處理流程，即采用外部的數(shù)字高程模型（Digital Elevation Model，DEM）數(shù)據(jù)模擬并去除地形相位。二軌法D-InSAR 的基本處理流程如圖2 所示。

圖2 二軌法D-InSAR 基本處理流程

1.3 概率積分法地表移動變形預(yù)計思想

地表移動與變形計算的概率積分法（PIM）是指把巖層移動看作服從統(tǒng)計規(guī)律的隨機(jī)過程，從而將開采引起的地表下沉剖面表示成概率密度函數(shù)積分公式的預(yù)計方法。已有大量理論推導(dǎo)及實(shí)例驗(yàn)證表明，當(dāng)?shù)叵麻_采工作面的頂板管理方法采用自然垮落時，概率積分法得到的地表沉降盆地穩(wěn)定后的預(yù)計結(jié)果與實(shí)測數(shù)據(jù)具有較高的擬合優(yōu)度和精度，是用于礦區(qū)開采沉陷預(yù)計的一種行之有效且廣泛使用的方法。

2 結(jié)果與分析

2.1 研究區(qū)工作面概況與遙感數(shù)據(jù)源

2.1.1 311 工作面概述

圖3 為研究區(qū)煤礦的311 工作面及其地表走向觀測線的沉降監(jiān)測點(diǎn)位示意圖。311 工作面位于該礦井Ⅱ1 采區(qū)上山西翼，工作面可采走向長550m、面長240m。該工作面整體呈一單斜構(gòu)造，煤層傾角8°～18°，平均傾角為10°，煤層厚度為2.6m～4.2m，平均采厚為3.0m。該工作面標(biāo)高-160m～-280m，地面標(biāo)高為+31m～+32.5m，2018 年2 月9 日開始回采，于2019 年7 月27 日開采結(jié)束。

圖3 井下開采的311 工作面及地表走向觀測線點(diǎn)位圖

2.1.2 遙感數(shù)據(jù)源

所使用的SAR 數(shù)據(jù)為Sentinel-1A（哨兵-1A）IW 成像模式SLC 數(shù)據(jù)，此SAR 數(shù)據(jù)為C 波段數(shù)據(jù)，空間分辨率為5m×20m（距離向×方位向），其他詳細(xì)參數(shù)見表1 所示。根據(jù)311 工作面開采的時間節(jié)點(diǎn)，獲取了2018 年2月8 日至2019 年9 月7 日間的46 景哨兵1A 影像數(shù)據(jù)。

表1 所用SAR 數(shù)據(jù)主要參數(shù)

2.2 基于D-InSAR 技術(shù)的開采沉降邊界分析

D-InSAR 處理的具體技術(shù)路線如圖4 所示，所涉及的主要技術(shù)內(nèi)容包括：（1）影像預(yù)處理。包括數(shù)據(jù)格式轉(zhuǎn)換和目標(biāo)區(qū)裁剪。主要是將原始影像數(shù)據(jù)處理為標(biāo)準(zhǔn)的SLC 數(shù)據(jù)，然后根據(jù)研究區(qū)范圍的大小，對影像進(jìn)行目標(biāo)區(qū)裁剪。（2）基線估算。估算結(jié)果包括時間基線、空間基線、多普勒偏移、一個相位變化周期代表的高程變化等信息，主要目的是評估像對干涉質(zhì)量。（3）生成干涉圖。首先根據(jù)影像的精密星歷計算影像中心偏移值，并設(shè)置合適的窗口大小以進(jìn)行自動配準(zhǔn)，精度可以達(dá)到1 個像元的水平；然后將主輔影像與已有的外部基準(zhǔn)DEM 進(jìn)行配準(zhǔn)，再經(jīng)過采樣處理之后使配準(zhǔn)精度達(dá)到0.001 個像元的水平；在此基礎(chǔ)上根據(jù)配準(zhǔn)好的主輔影像的相位信息進(jìn)行干涉，生成干涉條紋圖。（4）濾波和相干性計算。濾波處理主要目的是去除大氣延遲和平地效應(yīng)所引起的誤差相位，同時得到相干系數(shù)圖。（5）相位解纏。使影像相位信息與地形信息相對應(yīng)。（6）控制點(diǎn)選擇。控制點(diǎn)主要用于軌道精煉和重去平，常在去平后的干涉圖上進(jìn)行選擇，選擇應(yīng)遵循以下原則：避免在有地形相位未去除和變化的區(qū)域選擇控制點(diǎn)；選擇相干性高的區(qū)域；控制點(diǎn)應(yīng)分布于整個范圍內(nèi)；避免解纏錯誤的區(qū)域，如相位孤島等。（7）軌道精煉和重去平。目的是估算和去除殘余的恒定相位和解纏后還存在的相位坡道。（8）相位轉(zhuǎn)形變和地理編碼。解纏后的相位經(jīng)過轉(zhuǎn)換就可得到相對于解纏起始點(diǎn)的視線向形變信息。地理編碼的目的是將所有的結(jié)果編碼到所需的坐標(biāo)系下。

圖4 D-InSAR 技術(shù)路線圖

基于獲取的46 景哨兵1A 影像組成45 個干涉對，利用疊加二軌法D-InSAR 技術(shù)疊加了2018 年2 月8 日至2019 年9 月7 日的累積沉降，將每對D-InSAR 得到的沉降監(jiān)測結(jié)果按時序疊加即可獲得研究區(qū)內(nèi)的時序沉降，以此可確定311 工作面開采下沉影響邊界，如圖5 所示。311 工作面中心區(qū)域出現(xiàn)下沉梯度過大導(dǎo)致的地表形變失相關(guān)，因此監(jiān)測到的最大下沉為-302mm，其原因是因?yàn)檠退蛳鲁撂荻冗^大區(qū)域，InSAR 監(jiān)測失相干，無法監(jiān)測實(shí)際最大下沉量。

2.3 驗(yàn)證分析與下沉邊界的確定

以下沉值-10mm 作為開采沉陷盆地邊界，聯(lián)合D-InSAR、SBAS-InSAR 技術(shù)以及概率積分預(yù)計法，以多手段、多角度來綜合權(quán)衡確定311 工作面南側(cè)實(shí)際下沉邊界。

2.3.1 基于D-InSAR 技術(shù)確定的沉降邊界位置

如圖5 所示，利用D-InSAR 技術(shù)監(jiān)測到的下沉邊界大約在II311 工作面走向觀測線上的D013 點(diǎn)南側(cè)、距II311 工作面采空區(qū)邊緣的切眼位置平均距離為165m。

圖5 D-InSAR 技術(shù)確定的311 工作面實(shí)際下沉邊界

2.3.2 基于SBAS-InSAR 技術(shù)確定的沉降邊界位置

基于2018 年2 月8 日至2020 年9 月7 日的46 景哨兵1A 影像，再利用SBAS-InSAR 技術(shù)對311 工作面開采影響的下沉邊界進(jìn)行監(jiān)測，監(jiān)測到的下沉邊界大約在II311 工作面走向觀測線上的D013 點(diǎn)、距II311 工作面采空區(qū)邊緣的切眼位置平均距離為160m。

2.3.3 基于概率積分法確定的沉降邊界位置

如圖6 所示，可看出基于概率積分預(yù)計模型所獲得的下沉邊界在II311 工作面走向觀測線上的D013 點(diǎn)、距II311 工作面采空區(qū)邊緣的切眼位置平均距離為160m。

2.3.4 基于InSAR 及概率積分法預(yù)測的下沉邊界綜合確定

通過綜合上述D-InSAR、SBAS-InSAR 技術(shù)與概率積分預(yù)計三種方法分析311 工作面實(shí)際下沉邊界，得到如圖6 所示多方法確定的311 工作面實(shí)際下沉邊界位置，與實(shí)際觀測線上的點(diǎn)位相符合。

圖6 多方法綜合的311 工作面南側(cè)實(shí)際下沉邊界對比圖

3 結(jié)束語

工礦區(qū)地表形變監(jiān)測過程復(fù)雜，其破壞及危害性不可避免。傳統(tǒng)的變形監(jiān)測技術(shù)費(fèi)時費(fèi)力，難以得到大范圍內(nèi)的整體形變結(jié)果，隨著微波遙感的快速發(fā)展與應(yīng)用，D-InSAR 和SBAS-InSAR 等技術(shù)已廣泛應(yīng)用于地表沉降動態(tài)監(jiān)測。本研究主要是利用D-InSAR 技術(shù)對淮北礦區(qū)某煤礦311 工作面開采的采動影響邊界進(jìn)行確定研究，并通過概率積分法預(yù)計與SBAS-InSAR 干涉測量技術(shù)對D-InSAR 監(jiān)測結(jié)果進(jìn)行驗(yàn)證，提取到311 工作面開采后的沉降邊界，同時根據(jù)311 工作面走向觀測線的觀測點(diǎn)位，在走向觀測線上指定了沉降邊界的具體位置，即在D013 測點(diǎn)處，該測點(diǎn)至311 工作面采空區(qū)切眼巷道的平面距離為162m。