基于SBAS-InSAR的山東濟(jì)陽礦區(qū)沉降監(jiān)測(cè)與分析

潘光永，陶秋香，陳 洋，王 珂

(山東科技大學(xué)測(cè)繪科學(xué)與工程學(xué)院，山東 青島 266590)

0 引言

隨著國家經(jīng)濟(jì)和科技水平的飛速發(fā)展，礦產(chǎn)資源的需求量越來越大。我國煤炭資源不僅存儲(chǔ)量大，而且分布較為廣泛，是重要的礦產(chǎn)資源之一，其廣泛應(yīng)用于生產(chǎn)生活的各個(gè)領(lǐng)域[1-2]。然而煤炭資源的開采必然會(huì)帶來礦區(qū)及周圍地表的沉降與塌陷等一系列地質(zhì)災(zāi)害，為了降低地面沉降帶來的安全隱患，預(yù)防礦產(chǎn)資源的過度開采，對(duì)礦區(qū)地面沉降的監(jiān)測(cè)工作尤為重要[3-7]。

現(xiàn)有的地面沉降監(jiān)測(cè)方法主要有水準(zhǔn)測(cè)量、GPS(Global Position System)測(cè)量及合成孔徑雷達(dá)干涉測(cè)量(Interferometric Synthetic Aperture Radar，InSAR)等[8]。孟憲綱[9]曾經(jīng)采用水準(zhǔn)測(cè)量的方法對(duì)京滬高速鐵路進(jìn)行了沉降分析，測(cè)量結(jié)果精確可靠，但是水準(zhǔn)測(cè)量的方法在布設(shè)點(diǎn)位時(shí)難以做到分布均勻，更新周期較長(zhǎng)，并且需要很多人力物力支撐[10-11]。張訓(xùn)虎等[12]采用GPS測(cè)量的方法對(duì)高速鐵路進(jìn)行了沉降監(jiān)測(cè)，得出該方法可快速、高效的完成沉降監(jiān)測(cè)工作的結(jié)論。GPS的方法雖然具有較高的定位精度，但對(duì)高程信息不敏感，也難以滿足地表微小沉降量監(jiān)測(cè)的需求[13-14]。InSAR技術(shù)為新型遙感技術(shù)，其差分干涉模式D-InSAR技術(shù)具有全天時(shí)、全天候、大范圍等監(jiān)測(cè)地表沉降的特點(diǎn)，其形變監(jiān)測(cè)精度理論上可以達(dá)到厘米級(jí)甚至毫米級(jí)[15]。1989年Gabriel首次采用D-InSAR成功獲取了農(nóng)田區(qū)的形變信息，盡管學(xué)者們隨后對(duì)D-InSAR技術(shù)做了一系列改進(jìn)，但其精度易受時(shí)空失相干及大氣等因素的影響[16]。2001年Berardino等人在D-InSAR基礎(chǔ)上克服時(shí)間、空間失相干和大氣效應(yīng)的限制，發(fā)展了小基線集技術(shù)(Small Baseline Subset，SBAS)對(duì)長(zhǎng)時(shí)間緩慢沉降的監(jiān)測(cè)也有很好的效果[17]，因此成為了監(jiān)測(cè)地表沉降的主要手段[18-21]。本文基于SBAS-InSAR技術(shù)對(duì)覆蓋濟(jì)陽礦區(qū)2017年5月20日至2018年10月18日期間40景C波段Sentinel-1A升軌數(shù)據(jù)進(jìn)行干涉對(duì)的選取、去除大氣效應(yīng)、濾波、相位解纏、去平地效應(yīng)等處理，得到期間礦區(qū)的沉降量及平均沉降速率，結(jié)果表明，基于SBAS-InSAR技術(shù)可以形象直觀地再現(xiàn)礦區(qū)沉降變化過程，其沉降監(jiān)測(cè)結(jié)果與礦區(qū)實(shí)際開采進(jìn)程相符，可作為預(yù)防地質(zhì)災(zāi)害的有力依據(jù)。

1 研究區(qū)概況

濟(jì)陽井田礦區(qū)位于華北平原的南緣，黃河北岸，地勢(shì)平坦，局部地勢(shì)低洼，雨季易積水。黃河從井田南部1 km處流過，井田中東部有邢家渡引黃總干渠，中西部有大寺河通過。礦區(qū)東西向約9.9 km，南北向約5 km，面積約50 km2。本井田為全隱蔽的華北型石炭二迭系海陸交互相煤田。區(qū)內(nèi)地層自下而上為奧陶系石灰?guī)r、石炭系本溪組和太原組，二疊系山西組和石盒子組，第三系明化鎮(zhèn)組，第四系平原組。主要含煤地層為太原組和山西組。具體位置及范圍如圖1所示(紅色矩形框?yàn)榈V區(qū)范圍)。研究區(qū)主要為四采區(qū)7煤層，北靠DF15斷層(H=19 m∠47°)，南鄰李家斷層(H=70～260 m ∠70°)，西至煤7風(fēng)氧化帶邊界，東至DF16斷層(H=80 m ∠70°)，平均走向長(zhǎng)約2 378 m，傾斜長(zhǎng)約954 m，面積2.04×106m2。四采區(qū)內(nèi)7煤層位于太原組上部，上距一灰11.50～26.70 m，平均17.97 m，下距二灰30.50～40.30 m，平均34.72 m。下距第10層煤40.00～58.85 m，平均50.87 m。開采上限標(biāo)高-375.0 m左右，下限標(biāo)高-581.1 m左右，埋深398.2～604.3 m之間，平均506.3 m，煤厚0～1.81 m，平均1.1 m，傾角5°～12°，平均6°，工業(yè)儲(chǔ)量為2.596×106t，可采儲(chǔ)量為1.37×106t。其結(jié)構(gòu)較簡(jiǎn)單，屬穩(wěn)定型薄煤層。主要布設(shè)有4701、4701東、4701西、4702、4702西、4703、4704上、4704下和4705工作面，其分布如圖2所示，開采時(shí)間情況如表1所示。

圖1 研究區(qū)地理位置Fig.1 Geographical location of research area

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

2 SBAS-InSAR技術(shù)原理

常規(guī)的D-InSAR的測(cè)量方法受制于SAR影像的數(shù)量，且易受時(shí)空去相關(guān)的影響，時(shí)間、空間基線過長(zhǎng)等都會(huì)對(duì)差分的結(jié)果造成較大影響。由于礦區(qū)的開采過程較為緩慢，開采時(shí)間較長(zhǎng)，因此應(yīng)用D-InSAR的測(cè)量方法監(jiān)測(cè)長(zhǎng)時(shí)間段內(nèi)礦區(qū)形變結(jié)果具有一定局限性。小基線集技術(shù)則可通過利用多景SAR數(shù)據(jù)對(duì)長(zhǎng)時(shí)間的微小形變進(jìn)行高精度觀測(cè)，并通過產(chǎn)生的具有時(shí)間基線和空間基線特征相對(duì)較小的干涉圖序列，更好的克服去相關(guān)的影響，提高形變監(jiān)測(cè)質(zhì)量。

假設(shè)研究區(qū)域在(t0,…,tN)的時(shí)間段里共有N+1景雷達(dá)影像，根據(jù)干涉相對(duì)之間的自由組合，便可得到共M幅干涉圖，且滿足：

(1)

先將tA,tB兩景不同時(shí)刻的影像去除地形相位，假設(shè)tA

δφi(x,y)=φ(tB,x,y)-φ(tA,x,y)

(2)

式中φ(tA,x,y)和φ(tB,x,y)分別表示tA和tB時(shí)刻相對(duì)應(yīng)的相位，λ表示雷達(dá)波長(zhǎng)，d(tA,x,y)和d(tB,x,y)分別表示在tA和tB時(shí)刻任意像元(x,y)相對(duì)于基準(zhǔn)時(shí)刻t0在LOS方向上的地表形變。

將公式(2)中的相位φ表示成任意兩個(gè)獲取時(shí)間之間的平均相位速度v與時(shí)間t的乘積，以此來獲得具有物理意義的沉降序列，表達(dá)形式為：

(3)

于是，第i幅干涉圖的相位可表示為：

(4)

將其改寫成矩陣形式：

Bv=δφ

(5)

式中：B——M×N矩陣。最后，利用奇異值分解的方法得到速度矢量v的廣義逆解，再根據(jù)對(duì)各時(shí)間段速度在時(shí)間域上進(jìn)行積分便可得到各個(gè)時(shí)間段的形變量。

3 數(shù)據(jù)處理

研究采用2017年5月20日至2018年10月18日期間共計(jì)40景C波段Sentinel-1A升軌的雷達(dá)影像數(shù)據(jù)，高程數(shù)據(jù)(Digital Elevation Model，DEM)采用美國地質(zhì)調(diào)查局所提供的SRTM(Shuttle Radar Topography Mission，即航天飛機(jī)雷達(dá)地形測(cè)繪使命)90 m 分辨率DEM，實(shí)驗(yàn)通過SARscape軟件設(shè)置合適的時(shí)間閾值和空間閾值，共生成68個(gè)差分干涉對(duì)。差分干涉對(duì)的時(shí)空基線組合如圖3所示。

圖3 時(shí)空基線組合圖Fig.3 Time-space baseline combination map

3.1 差分干涉處理

利用獲取的DEM數(shù)據(jù)對(duì)生成的68個(gè)干涉對(duì)進(jìn)行差分處理，生成相干性干涉圖，通過Goldstein的濾波方法去除噪聲影響，進(jìn)行增強(qiáng)處理，濾波后的部分差分干涉圖如圖4所示。采用最小費(fèi)用流(Minimum Cost Flow)的方法進(jìn)行相位解纏，保證影像中相對(duì)孤立的高相干區(qū)域有更好的處理結(jié)果，獲取研究區(qū)域真實(shí)相位。

3.2 軌道精煉和重去平

為去除矯正相位偏移，選取位于沒有形變或形變較小的區(qū)域且具有高相干性的GCP點(diǎn)用于進(jìn)行修正干涉相位和解纏后的相位，進(jìn)而利用線性模型對(duì)初始位移進(jìn)行估計(jì)，去除殘余地形。

3.3 SBAS反演

為了更好的分析形變監(jiān)測(cè)的結(jié)果，首先估算形變速率與殘余地形，然后對(duì)差分干涉圖進(jìn)行二次解纏處理，進(jìn)而對(duì)干涉圖進(jìn)行優(yōu)化，得到第一次反演結(jié)果。在第一次反演得到的形變速率的基礎(chǔ)上，利用軌道精煉中所選的GCP點(diǎn)移除恒定相位或斜坡相位，根據(jù)大氣濾波的定制，估算并去除大氣效應(yīng)，最終得到時(shí)間序列上的位移變化。獲取該礦區(qū)年平均沉降速率圖(圖5)。

圖4 濾波后的干涉圖示例Fig.4 Example of filtered interferogram

圖5 年平均沉降速率Fig.5 Annual average sedimentation rate

4 結(jié)果分析

4.1 與礦區(qū)開采進(jìn)程的對(duì)比分析

由圖5沉降監(jiān)測(cè)結(jié)果分析可知：

(1)監(jiān)測(cè)期間該研究區(qū)出現(xiàn)三處沉降盆地，分別位于研究區(qū)西北、西南和東南方向的工作面附近，其中4705工作面與研究區(qū)西北方向4701和4702工作面附近沉降最為嚴(yán)重，其最大年平均沉降速率可達(dá)320 mm/a，礦區(qū)的沉降速率由沉降中心到邊緣逐漸減小。

(2)為進(jìn)一步獲取監(jiān)測(cè)期間礦區(qū)的累積沉降結(jié)果，對(duì)各時(shí)間段沉降速率進(jìn)行時(shí)間域上的積分，獲取2017年5月20日至2018年10月18日期間共39幅累積沉降量圖。圖6給出了研究區(qū)部分成像時(shí)刻相對(duì)于2017年5月20日的累積沉降量圖。

圖6 累積沉降量Fig.6 Cumulative settlement

由圖6結(jié)合煤層各個(gè)工作面實(shí)際開采情況分析得到：

①2017年5月至2017年12月期間，隨著開采工作的進(jìn)行，4702工作面上方沉降持續(xù)增大，最大沉降量達(dá)到130 mm；4701西工作面隨開采工作的開始，上方出現(xiàn)下沉，并且該工作面南鄰剛剛結(jié)束開采的4701工作面，其上覆巖層尚未達(dá)到力學(xué)平衡狀態(tài)，仍在緩慢下沉，從而導(dǎo)致地下巖層結(jié)構(gòu)更趨于不穩(wěn)定狀態(tài)，加劇了沉降現(xiàn)象，最大沉降量達(dá)220 mm。

②2017年12月至2018年5月期間，4701西、4703及4701東工作面開采工作正在進(jìn)行。根據(jù)圖6(c)與圖6(d)可看出：隨著4701西工作面開采工作繼續(xù)進(jìn)行，西南處沉降盆地持續(xù)下沉，沉降范圍有所擴(kuò)大，最大沉降量達(dá)到390 mm；4703工作面北鄰剛剛結(jié)束開采的4702工作面，隨著其開采工作的進(jìn)行，加劇了該區(qū)域的地面沉降，最大沉降量達(dá)300 mm。

③2018年5月至2018年11月期間，4701東工作面因開采工作繼續(xù)進(jìn)行，致使該區(qū)域沉降量及沉降區(qū)域逐漸增大，至2018年10月份開采工作結(jié)束，該區(qū)域沉降量最大達(dá)到300 mm，并且因該區(qū)域地下煤層開采導(dǎo)致東北方向地面也發(fā)生了沉降現(xiàn)象；4702西工作面從2018年8月份開始開采后便加劇了該地區(qū)的沉降，最大沉降量達(dá)到447 mm；由于4701西、4702和4703三個(gè)工作面的開采，使采區(qū)附近巖層結(jié)構(gòu)遭到破壞，原有的力學(xué)平衡結(jié)構(gòu)逐步變得不穩(wěn)固，并使開采區(qū)上覆巖層發(fā)生移動(dòng)和變形，從而影響周圍地區(qū)的巖層結(jié)構(gòu)，導(dǎo)致4705工作面附近在正式開采前就發(fā)生沉降現(xiàn)象，加上附近的4704工作面原本就進(jìn)行過開采工作，地下巖層結(jié)構(gòu)早已破壞，因此受附近三面開采工作的影響使該地區(qū)沉降更為嚴(yán)重，在2018年8月4705工作面開采后又加劇了沉降，至2018年10月18日沉降量最大達(dá)到447 mm。

(3)分別在沉降較小的4701東工作面附近選取A點(diǎn)，沉降較為嚴(yán)重的4705工作面附近選取B點(diǎn)進(jìn)行時(shí)間序列分析，具體位置如圖6(f)所示，結(jié)果如圖7(a)和圖7(b)所示。

圖7 時(shí)間序列分析Fig.7 Time series analysis

由圖7分析可知，兩點(diǎn)處隨著開采工作的進(jìn)行，總體上沉降量越來越大且沒有減緩的趨勢(shì)。其中圖7(a)可看出，監(jiān)測(cè)初期，由于4701西和4702工作面的開采，A點(diǎn)處開始出現(xiàn)沉降，至2017年8月左右，沉降量達(dá)45 mm，之后因?yàn)樵搮^(qū)域巖層硬度及4701西工作面開采強(qiáng)度減弱等因素致使該處沉降量明顯減少。至2018年3月，4701東工作面開采活動(dòng)的進(jìn)行導(dǎo)致A點(diǎn)處再次出現(xiàn)較大下沉，隨著開采活動(dòng)的持續(xù)進(jìn)行，沉降量逐漸增大，至2018年10月18日，最大沉降量大于140 mm。

圖7(b)可看出在4705工作面附近因受三個(gè)方向工作面開采的影響，也在持續(xù)發(fā)生沉降，并且隨著2018年5月4705工作面開采后，沉降量增大，最大大于400 mm。

4.2 與水準(zhǔn)監(jiān)測(cè)結(jié)果的對(duì)比分析

為了驗(yàn)證SBAS-InSAR監(jiān)測(cè)結(jié)果的可靠性，分別在沉降盆地邊緣及靠近沉降中心處選取了7個(gè)水準(zhǔn)點(diǎn)的累積沉降結(jié)果與SBAS-InSAR監(jiān)測(cè)的結(jié)果進(jìn)行對(duì)比分析。選取的水準(zhǔn)點(diǎn)布設(shè)情況如圖8所示，對(duì)比結(jié)果詳見表2。

圖8 水準(zhǔn)點(diǎn)位置Fig.8 Location of bench marks

表2 SBAS-InSAR形變監(jiān)測(cè)結(jié)果與水準(zhǔn)測(cè)量結(jié)果差異比較

根據(jù)表2的對(duì)比結(jié)果可看出，在本次實(shí)驗(yàn)中所得到的形變結(jié)果與水準(zhǔn)測(cè)量得到的形變結(jié)果吻合較好，兩種技術(shù)監(jiān)測(cè)結(jié)果最小差值為-7.71 mm，最大差值為28.12 mm，經(jīng)計(jì)算平均差值約為5.91 mm。由于SAR影像自身特性的限制、干涉過程中受到各種誤差因素的影響，SBAS-InSAR的監(jiān)測(cè)結(jié)果與實(shí)際水準(zhǔn)數(shù)據(jù)有所差別，但從表2可看出大部分點(diǎn)形變的差值相對(duì)較小，表明SBAS-InSAR監(jiān)測(cè)的結(jié)果是可靠的。

5 結(jié)論

本文基于C波段Sentinel-1A升軌的雷達(dá)影像數(shù)據(jù)，采用SBAS-InSAR技術(shù)方法對(duì)山東濟(jì)陽某礦區(qū)進(jìn)行了監(jiān)測(cè)，分別從研究區(qū)內(nèi)地面沉降的年平均沉降速率、累積沉降量以及沉降較為嚴(yán)重區(qū)域的時(shí)間序列進(jìn)行了分析。結(jié)果表明，2017年5月20日至2018年10月18日期間，研究區(qū)沉降最為嚴(yán)重的區(qū)域集中在4701西與4702西工作面和4705工作面，年平均沉降速率最大達(dá)到320 mm/a，累積沉降量最大為447 mm，且隨著時(shí)間推移，沉降量逐漸增大。此外，位于礦區(qū)周邊的104國道與濟(jì)南繞城高速受開采區(qū)影響，也發(fā)生了部分沉降，說明了在研究區(qū)工作面進(jìn)行開采工作后，不僅開采地區(qū)會(huì)有沉降現(xiàn)象發(fā)生，周圍地區(qū)也受其影響而出現(xiàn)下沉。通過與礦區(qū)實(shí)際開采情況的對(duì)比分析，表明SBAS-InSAR的方法在礦區(qū)沉降監(jiān)測(cè)上有很大的適用性，能夠快速準(zhǔn)確地定位礦區(qū)沉降的位置及范圍，且其監(jiān)測(cè)結(jié)果與礦區(qū)開采進(jìn)程相符，并對(duì)因礦區(qū)開采而導(dǎo)致的地質(zhì)災(zāi)害問題的預(yù)防有重要意義。