基于SBAS-InSAR 技術(shù)的楊伙盤礦區(qū)地表形變監(jiān)測(cè)

張 帆，常 樂，荀張媛

（1.西安煤航遙感信息有限公司，陜西 西安 710100；2.沈陽城市建設(shè)學(xué)院，遼寧 沈陽 110167）

0 引言

我國是礦產(chǎn)資源大國，煤炭等資源一直在能源開發(fā)利用中占據(jù)重要地位。煤炭資源的長(zhǎng)期開采不僅引起地面塌陷等一系列問題，還會(huì)造成建筑物受損、農(nóng)田被破壞[1]，以及誘發(fā)崩塌或滑坡等次生災(zāi)害[2]，對(duì)人民的生產(chǎn)生活造成嚴(yán)重的經(jīng)濟(jì)損失，甚至?xí)＜吧踩玔3-4]。因此，重視礦區(qū)地表形變，實(shí)現(xiàn)礦區(qū)生態(tài)環(huán)境可持續(xù)發(fā)展是科學(xué)開采煤炭等礦產(chǎn)資源的重要任務(wù)。

礦區(qū)地表形變是一個(gè)基于面區(qū)域的大范圍的地表運(yùn)動(dòng)，但是在實(shí)際生產(chǎn)過程中，更多的還是依賴于傳統(tǒng)的水準(zhǔn)測(cè)量、GPS 測(cè)量等點(diǎn)監(jiān)測(cè)技術(shù)，在生產(chǎn)工作面上每隔一段距離布設(shè)一個(gè)水準(zhǔn)點(diǎn)或GPS 點(diǎn)。一方面，布設(shè)較多的點(diǎn)，人工消耗成本高、周期長(zhǎng)；另一方面，水準(zhǔn)測(cè)量和GPS 測(cè)量都是基于點(diǎn)觀測(cè)，對(duì)于礦區(qū)大范圍大區(qū)域監(jiān)測(cè)，需要大量加密點(diǎn)，通過擬合獲取最終結(jié)果，無法保證精度[5-7]。

合成孔徑雷達(dá)干涉測(cè)量（InSAR）技術(shù)是一種高效對(duì)地觀測(cè)技術(shù)，具有全天時(shí)、全天候?qū)Φ爻上?，監(jiān)測(cè)結(jié)果精確度高，測(cè)量范圍廣的優(yōu)勢(shì)，大大降低人力、物力、財(cái)力的消耗[8-10]。特別是時(shí)間序列InSAR 技術(shù)能夠更快速獲取地表大范圍高精度時(shí)間序列形變信息，為礦區(qū)開采沉陷監(jiān)測(cè)提供了新的方向[11-13]。

越來越多的學(xué)者已經(jīng)將InSAR 技術(shù)應(yīng)用于煤礦地表監(jiān)測(cè)研究。尹宏杰等[14]利用時(shí)間序列InSAR 技術(shù)對(duì)湖南冷水江錫煤礦地面形變進(jìn)行監(jiān)測(cè)，證實(shí)了煤礦地表時(shí)間序列形變漏斗的發(fā)展過程；趙偉穎等[15]利用SBAS-InSAR 技術(shù)獲取了徐州礦區(qū)地表的時(shí)間序列形變及累積形變量，并與水準(zhǔn)測(cè)量結(jié)果進(jìn)行了外部精度符合驗(yàn)證，且兩者的相關(guān)性較高；張香凝等[16]采用SBAS-InSAR 技術(shù)獲取了煤礦地表形變的時(shí)空分布特征，并利用數(shù)值模擬技術(shù)對(duì)觀測(cè)形變結(jié)果進(jìn)行模擬分析，獲取了地面沉降在時(shí)間和空間上的變形規(guī)律和機(jī)制；姚佳明等[17]選用升軌、降軌L波段PALSAR-2 影像數(shù)據(jù)，利用InSAR 技術(shù)對(duì)煤礦采空區(qū)開展了短期動(dòng)態(tài)地表沉降監(jiān)測(cè)，結(jié)合研究區(qū)開采信息對(duì)煤礦采空范圍及開采時(shí)間進(jìn)行反演，驗(yàn)證了InSAR 技術(shù)對(duì)煤礦采區(qū)反演的合理性與可靠性；何榮興等[18]介紹了采空區(qū)災(zāi)害類型及不同類型的相應(yīng)案例，分析了采空區(qū)災(zāi)害發(fā)生的一般規(guī)律和特征，提出了盡量選擇不產(chǎn)生采空區(qū)的采礦方法。

本文以楊伙盤煤礦為例，采用SBAS-InSAR 技術(shù)，利用中分辨率Sentinel-1A 數(shù)據(jù)對(duì)煤礦地表進(jìn)行時(shí)間序列形變監(jiān)測(cè)分析，提取礦區(qū)開采沉陷形變場(chǎng)，并結(jié)合光學(xué)影像數(shù)據(jù)，研究其沉陷變化規(guī)律。

1 研究區(qū)概況

陜北煤炭生產(chǎn)基地是全國煤炭生產(chǎn)基地之一，尤其是神木、府谷等區(qū)域，煤炭資源豐富。但是長(zhǎng)期高強(qiáng)度的資源開采，加之陜北地區(qū)大多為黃土地貌，最終導(dǎo)致礦區(qū)的土地大面積沉陷，并引發(fā)了滑坡、崩塌等多種次生災(zāi)害。因此，為保護(hù)礦區(qū)生態(tài)環(huán)境，減少人民的財(cái)產(chǎn)損失，需要對(duì)礦區(qū)的采空區(qū)形變進(jìn)行長(zhǎng)時(shí)間的有效監(jiān)測(cè)；結(jié)合多種技術(shù)，加大對(duì)采煤塌陷區(qū)域的監(jiān)測(cè)監(jiān)管力度，形成一套行之有效的陜北榆神府礦區(qū)采空區(qū)沉陷監(jiān)測(cè)體系。

楊伙盤煤礦位于陜西省神木市城北約30 km 的府店一級(jí)公路北側(cè)，行政區(qū)劃屬神木市店塔鎮(zhèn)，靠近府谷縣（圖1），工作區(qū)面積約26.92 km2。礦區(qū)地表廣泛覆蓋著現(xiàn)代風(fēng)積沙及第四系黃土層、新近系紅土層等，僅部分地區(qū)有基巖出露，基巖以砂巖、泥巖為主。根據(jù)某一鉆孔資料，3-1煤層主要位于侏羅系中統(tǒng)延安組第三段，煤層距地表約160 m。根據(jù)現(xiàn)有巷道布置情況，結(jié)合各煤層開采范圍及煤層賦存特點(diǎn)，一水平3-1煤層劃分4 個(gè)盤區(qū)。圖2 展示了3-1煤層開采工作面分布情況。

圖1 楊伙盤煤礦地理位置及高程Fig.1 Geographical location and elevation of Yanghuopan Coal Mine

圖2 工作面分布圖Fig.2 Distribution map of working face

2 InSAR 數(shù)據(jù)及技術(shù)方法

2.1 InSAR 技術(shù)方法

BERARDINO 等[19]提出的SBAS-InSAR 技術(shù)，是目前較為常用的一種時(shí)間序列InSAR 技術(shù)，技術(shù)流程如圖3 所示。首先，根據(jù)SAR 數(shù)據(jù)量，通過設(shè)置時(shí)間基線閾值和空間基線閾值組合干涉對(duì)，一般根據(jù)相干系數(shù)圖和干涉圖挑選質(zhì)量較高的干涉對(duì)；其次，通過多視和濾波等操作提高所選干涉圖的信噪比，基于具有穩(wěn)定后向散射特性的高相干點(diǎn)去除地形相位和平地相位，并進(jìn)行解纏；之后，根據(jù)干涉相位及時(shí)空基線去除DEM 誤差，并分離線性形變相位，解纏殘余相位，通過時(shí)間域高通濾波和空間域低通濾波分離出殘余相位中的大氣相位，剩下即為非線性相位，補(bǔ)償線性形變相位即可得到完整的形變相位；最后，采用奇異值分解，連接多個(gè)子集，求得影像序列地表形變速率的最小范數(shù)、最小二乘解[20]。該方法能有效去除或減弱時(shí)空失相干等誤差的影響，并在一定程度消除大氣相位和限制地形誤差的影響。

圖3 SBAS 技術(shù)流程圖Fig.3 Flow chart of SBAS technology

獲取研究 區(qū)t0······tn時(shí)間段內(nèi)N+1幅SAR 影像，將影像配準(zhǔn)到同一個(gè)坐標(biāo)系，設(shè)定時(shí)間基線閾值及空間基線閾值，得到M個(gè)差分干涉對(duì)。假設(shè)其中第i幅干涉圖由tA、tB（tA＜tB）兩時(shí)期影像生成，此干涉圖上雷達(dá)坐標(biāo)系下（a,r）處的干涉相位可表示為式（1）。

式中：δφi（disp）為形變相位；δφi（topo）為殘余地形相位；δφi（aps）為大氣延遲相位；δφi（noise）為噪聲相位。忽略殘余地形相位、大氣延遲相位及噪聲相位，將形變相位以LOS向形變表示，式（1）可簡(jiǎn)化為式（2）。

式中：λ為雷達(dá)波長(zhǎng)；dtA（a,r）、dtB（a,r）為（a,r）處在tA時(shí)期、tB時(shí)期相對(duì)于初始t0時(shí)期的形變。針對(duì)（a,r）這一像元處，假設(shè)IE為主影像，IS為從影像，對(duì)應(yīng)全部的M個(gè)干涉對(duì)，見式（3）。

式（3）滿足IEi＞ISi，其中，i=1,···,M。則對(duì)于任意干涉圖，對(duì)應(yīng)的有式（4）。

式（4）為n個(gè)未知數(shù)M個(gè)方程所組成的方程組，以矩陣形式表示 δφ=Aφ，其中，A為M*n矩陣。對(duì)于小基線子集內(nèi)部，M≥n，方程可由最小二乘法求解。各子集之間M接近或者小于n時(shí)方程會(huì)出現(xiàn)多解，采用奇異值分解以及最小范數(shù)進(jìn)行解算獲得唯一解[19]。

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

本次研究使用的數(shù)據(jù)為歐洲航天局Sentinel-1A衛(wèi)星在TOPS 成像模式下的寬幅干涉SLC 數(shù)據(jù)，為目前常用的免費(fèi)中分辨率SAR 數(shù)據(jù)，數(shù)據(jù)產(chǎn)品詳細(xì)參數(shù)見表1。根據(jù)監(jiān)測(cè)周期，共收集了2020 年8 月—2022 年11 月共58 景升軌數(shù)據(jù)，軌道號(hào)為113/126。

表1 Sentinel-1A IW 產(chǎn)品基本參數(shù)Table 1 Basic parameters of Sentinel-1A IW product

本次采用的DEM 數(shù)據(jù)為ALOS Global Digital Surface Model“ALOS World 3D-30 m”（AW3D30），是由日本宇宙航空研究開發(fā)機(jī)構(gòu)（JAXA）2015 年5 月免費(fèi)提供的高精度全球數(shù)字地表模型數(shù)據(jù)，水平分辨率為30 m（1"），高程精度5 m，是目前世界上最精確的3D 地圖，覆蓋全球所有的土地尺度，工作區(qū)的DEM 數(shù)據(jù)如圖1 所示。精密軌道數(shù)據(jù)采用成像21 d之后發(fā)布的POD 精密軌道數(shù)據(jù)。

3 結(jié)果與分析

3.1 InSAR 形變分析

工作區(qū)年平均形變速率圖如圖4 所示。形變部分主要分布在①、②、③三個(gè)區(qū)域，沿著工作面均表現(xiàn)為長(zhǎng)條狀分布，主要覆蓋的工作面有30112 工作面、30114 工作面、30116 工作面、30115 工作面、30117工作面、30301 工作面和30302 工作面。三個(gè)形變中心可探測(cè)的最大形變速率分別約為-52 mm/a、-39 mm/a 和-39 mm/a，形變最大的區(qū)域位于區(qū)域①，分別對(duì)三個(gè)形變區(qū)進(jìn)行時(shí)間序列和剖線形變分析。

圖4 工作區(qū)年平均形變速率圖Fig.4 Annual average deformation rate of the working area

圖5 為區(qū)域①的時(shí)間序列形變圖。由圖5 可知，從2020 年8 月開始，30112 工作面處于開采狀態(tài)，圖5（a）為基準(zhǔn)面，在圖5（b）中地表開始出現(xiàn)形變。隨著開采面的推進(jìn)，地表形變范圍逐漸擴(kuò)大，直至監(jiān)測(cè)期結(jié)束，最終形變主要覆蓋的工作面為30112 工作面、30114 工作面和30116 工作面。其中，30114 工作面地表形變量級(jí)及范圍都較大，最大累積形變量約為-130 mm，為整個(gè)工作區(qū)的可探測(cè)到的最大形變量，在30114 工作面上（圖4）畫一條剖線AA′以及AA′剖線上的形變最值特征點(diǎn)a 進(jìn)行時(shí)間序列分析。

圖5 區(qū)域①時(shí)間序列形變圖Fig.5 Time series deformation of region ①

圖6 為剖線AA′的時(shí)間序列剖線圖，其中主要選擇了6 個(gè)時(shí)間段顯示，共形成了4 個(gè)形變漏斗。由圖6 可知，靠近A 點(diǎn)的第一個(gè)形變漏斗的變化相對(duì)較為平穩(wěn)；第二個(gè)形變漏斗在2021 年1 月后發(fā)生突變；第三個(gè)漏斗在2021 年5 月后發(fā)生突變；第四個(gè)形變漏斗在2021 年8 月發(fā)生突變，這三個(gè)突變情況與工作面的開采時(shí)間都較為吻合。其中，第四個(gè)形變漏斗的形變量遠(yuǎn)遠(yuǎn)超過其他區(qū)域的形變量。圖7 為剖線AA′上的形變最值點(diǎn)a 的時(shí)間序列形變折線圖。由圖7 可知，2021 年8 月，該點(diǎn)發(fā)生了突變，而后速率逐漸減小。2022 年5 月前后，再次出現(xiàn)較小幅度的突變，可能是受到30116 工作面開采的影響。

圖6 AA'時(shí)間序列剖線圖Fig.6 Time series deformation of AA' profile

圖7 a 點(diǎn)時(shí)間序列圖Fig.7 Time series deformation of point a

圖8 為區(qū)域②的時(shí)間序列形變圖。由圖8 可知，從2020 年8 月開始，30115 工作面處于開采狀態(tài)。由于初期的形變較小，絕對(duì)量值在20 mm 之內(nèi)，在圖8（d）中才顯示出來。最終形變主要覆蓋30115 工作面和30117 工作面。在30117 工作面上（圖4）畫一條剖線BB′以及BB′剖線上的形變最值特征點(diǎn)b 進(jìn)行時(shí)間序列分析。

圖8 區(qū)域②時(shí)間序列形變圖Fig.8 Time series deformation in region ②

圖9 為剖線BB′的時(shí)間序列剖線圖，其中，選擇了6 個(gè)時(shí)間段顯示，共形成多個(gè)形變漏斗。由圖9 可知，靠近B 點(diǎn)的形變突變基本都發(fā)生在2021 年1 月之后，后面時(shí)間段的形變變化相對(duì)都較為平穩(wěn)。圖10 為剖線BB′上的形變最值點(diǎn)b 的時(shí)間序列形變折線圖。由圖10 可知，在2021 年3 月前后，b 點(diǎn)發(fā)生了突變，與工作面開采的時(shí)間基本吻合，4 月之后速率逐漸減小。

圖9 BB'時(shí)間序列剖線圖Fig.9 Time series deformation of BB' profile

圖10 b 點(diǎn)時(shí)間序列圖Fig.10 Time series deformation of point b

圖11 為區(qū)域③的時(shí)間序列形變圖。由圖11 可知，2021 年此區(qū)域的30301 工作面才開始開采；2021年9 月，工作面的形變才開始顯示。最終，形變主要覆蓋的工作面為30301 工作面和30302 工作面。在30302 工作面上（圖4）畫一條剖線CC′以及CC′剖線上的形變最值特征點(diǎn)c 進(jìn)行時(shí)間序列分析。

圖11 區(qū)域③時(shí)間序列形變圖Fig.11 Time series deformation in region ③

圖12 為剖線CC′的時(shí)間序列剖線圖，其中，選擇了6 個(gè)時(shí)間段顯示，只形成1 個(gè)形變漏斗，漏斗中心基本靠近剖線CC′的中間位置。由于剖線的位置靠近30301 工作面，在2021 年8 月和2022 年1 月，剖線上的形變先后發(fā)生兩次突變，第一次突變是受到30301 工作面開采的影響，第二次突變是受到30302工作面開采的影響。圖13 為剖線CC′上的形變最值點(diǎn)c 的時(shí)間序列形變折線圖。由圖13 可知，在2021年9 月前后及2022 年3 月前后，c 點(diǎn)發(fā)生了兩次突變，與剖線整體的形變特征基本吻合。

圖12 CC'時(shí)間序列剖線圖Fig.12 Time series deformation of CC' profile

圖13 c 點(diǎn)時(shí)間序列圖Fig.13 Time series deformation of point c

隨著時(shí)間的推移，礦區(qū)內(nèi)主要形變區(qū)域的形變愈加明顯，時(shí)序InSAR 監(jiān)測(cè)結(jié)果數(shù)據(jù)對(duì)比反映了礦區(qū)內(nèi)局部地區(qū)的地表變形情況。從獲取的結(jié)果來看，礦區(qū)的形變區(qū)域是時(shí)刻發(fā)生變化的，這是因?yàn)榈V區(qū)開采造成的形變隨著時(shí)間的推移而逐漸演化，可以清楚地看到礦區(qū)形變的緩慢變化過程，礦區(qū)其他區(qū)域地表大部分處于穩(wěn)定狀態(tài)。

3.2 精度評(píng)定

精度評(píng)定是驗(yàn)證獲取地表形變量精度的一種方法，一般分為內(nèi)符合精度評(píng)定和外符合精度評(píng)定。內(nèi)符合精度評(píng)定一般是獲取的形變的中誤差（標(biāo)準(zhǔn)差）、形變量分布直方圖或者不同數(shù)據(jù)源之間形變值的相互比較。外符合精度評(píng)定一般是借助外部的形變數(shù)據(jù)，如GPS 數(shù)據(jù)或水準(zhǔn)數(shù)據(jù)，與InSAR 獲取的形變數(shù)據(jù)進(jìn)行對(duì)比驗(yàn)證。

由于本次試驗(yàn)沒有收集到外部數(shù)據(jù)，且SAR 數(shù)據(jù)源只有Sentinel-1A。因此，通過形變量值的標(biāo)準(zhǔn)差的計(jì)算，進(jìn)行內(nèi)符合精度評(píng)定。圖14 為形變速率像元分布圖。由圖14 可知，85.26%的像元的形變速率的絕對(duì)值小于10 mm。圖15 為工作區(qū)形變速率標(biāo)準(zhǔn)差分布圖。由圖15 可知，大部分區(qū)域的像元的標(biāo)準(zhǔn)差均在6 mm 范圍內(nèi)，在形變量級(jí)大的區(qū)域的像元的標(biāo)準(zhǔn)差相對(duì)較大，在6～15 mm 范圍之間。

圖14 形變速率像元分布圖Fig.14 Pixel distribution of deformation rate

圖15 形變速率標(biāo)準(zhǔn)差分布圖Fig.15 Distribution map of standard deviation of deformation rate

3.3 光學(xué)影像解譯

通過無人機(jī)正射飛行獲取到區(qū)域①的高分辨率光學(xué)影像，將InSAR 形變與光學(xué)影像疊加到一起進(jìn)行光學(xué)解譯，如圖16 所示。礦區(qū)地面形變的表現(xiàn)特征一般為地面塌陷造成的地裂縫和塌陷坑等。

圖16 區(qū)域①光學(xué)影像與形變疊加圖Fig.16 Optical image and deformation superposition in region ①

由圖16 可知，區(qū)域①的形變主要為（a）（b）（c）三部分，地表形變最主要的特征表現(xiàn)為地裂縫展布，部分區(qū)域分布有多個(gè)塌陷坑，而形變量級(jí)較小的區(qū)域地表的光學(xué)特征不甚明顯。結(jié)果顯示，光學(xué)影像特征與InSAR 監(jiān)測(cè)的形變結(jié)果的分布較為吻合，但I(xiàn)nSAR 技術(shù)更適合小量級(jí)的形變監(jiān)測(cè)，尤其是光學(xué)影像上無明顯特征的區(qū)域。

4 結(jié)論

本文基于SBAS-InSAR 方法對(duì)楊伙盤煤礦在2020 年8 月至2022 年10 月間的地表進(jìn)行了形變監(jiān)測(cè)，并利用高分辨率光學(xué)影像進(jìn)行了特征解譯。

1）礦區(qū)范圍內(nèi)共分布了三處形變區(qū)域，均呈長(zhǎng)條狀分布，與工作面開采范圍高度吻合。整個(gè)區(qū)域可探測(cè)的最大形變值達(dá)到-130 mm，最大年平均沉陷速率約為-52 mm/a。

2）通過時(shí)間序列形變分析，每處形變區(qū)在對(duì)應(yīng)或相鄰的工作面開采時(shí)發(fā)生突變，持續(xù)時(shí)間在1～2個(gè)月之間，隨即開始緩慢的持續(xù)性形變。

3）通過高分辨率光學(xué)影像和InSAR 形變結(jié)果的綜合分析，在形變量級(jí)較大的區(qū)域，在光學(xué)影像上均展布多條地裂縫，部分區(qū)域還分布有塌陷坑，表明InSAR 監(jiān)測(cè)結(jié)果的可靠性。而形變量級(jí)較小的區(qū)域，光學(xué)影像上無明顯的特征，表明InSAR 技術(shù)在地表小量級(jí)形變監(jiān)測(cè)方面的優(yōu)勢(shì)。