基于MT-InSAR和灰色預(yù)測模型的金屬礦區(qū)地表變形研究*

王 華 ，李素敏,2,3，袁利偉，張 瑋，馬相松，管慶丹

(1.昆明理工大學(xué) 國土資源工程學(xué)院，云南 昆明650093；2.云南省高校高原山區(qū)空間信息測繪技術(shù)應(yīng)用工程研究中心，云南 昆明650093；3.中國有色金屬工業(yè)協(xié)會智慧礦山地理空間信息集成創(chuàng)新重點(diǎn)實(shí)驗(yàn)室，云南 昆明650093；4．昆明理工大學(xué) 公共安全與應(yīng)急管理學(xué)院，云南 昆明650093；5.玉溪大紅山礦業(yè)有限公司，云南 玉溪 653400；6.昆明有色冶金設(shè)計研究院股份公司，云南 昆明 650000)

0 引言

金屬礦山由于高強(qiáng)度和大面積的地下開采，導(dǎo)致巖層結(jié)構(gòu)遭到破壞,地表沉陷現(xiàn)象頻發(fā)，易對礦區(qū)基礎(chǔ)設(shè)施和周邊生態(tài)環(huán)境造成不良影響[1]。掌握金屬礦山重大風(fēng)險災(zāi)害總體空間分布特征，明確其發(fā)展動態(tài)，并進(jìn)行遙感監(jiān)測和模型預(yù)測分析，對礦山安全監(jiān)管與預(yù)測預(yù)警研究具有重要的現(xiàn)實(shí)意義。因此，近年來學(xué)者針對礦區(qū)地表變形監(jiān)測和預(yù)測模型做了大量研究[2-3]。

礦區(qū)變形監(jiān)測傳統(tǒng)方法有水準(zhǔn)測量、GPS測量等，雖然這些方法精度高，但是由于技術(shù)的局限性，僅能監(jiān)測地面“離散點(diǎn)”，很難反映地表大面積的“面域”變形信息和時空演化特征，且需耗費(fèi)大量的人力、物力和財力[4]。近些年，在微波遙感領(lǐng)域，合成孔徑雷達(dá)干涉測量(Interferometric Synthetic Aperture Radar, InSAR)已成為對地觀測的研究熱點(diǎn)之一，InSAR技術(shù)具有全天時、全天候的工作能力和高精度、高分辨率、覆蓋范圍廣等特點(diǎn)，并能對地連續(xù)觀測獲取地表“面域”的變形信息[5]。時間序列InSAR技術(shù)(Multiple Temporal InSAR, MT-InSAR)是在D-InSAR的基礎(chǔ)上發(fā)展而來的。小基線集SABS技術(shù)最先由BERARDINO等[6]在2002年提出。李達(dá)等[7]利用SBAS和D-InSAR兩種方法獲取了某礦區(qū)的累計沉降值，并對其進(jìn)行差值分析和量化分析，為礦區(qū)地表的變形監(jiān)測與分析提供了新手段。目前對于礦區(qū)地表的監(jiān)測研究成果已有很多，如：胡海斌等[8]運(yùn)用灰色模型對沉陷區(qū)的沉降量進(jìn)行了預(yù)測，經(jīng)實(shí)例驗(yàn)證,灰色預(yù)測模型能滿足工程預(yù)測精度要求；李路等[9]加入權(quán)值思想獲得了改進(jìn)后的灰色系統(tǒng)模型，并將其運(yùn)用于實(shí)際沉降預(yù)測中；李金超等[10]提出了一種基于灰色支持向量機(jī)的組合預(yù)測模型并對礦區(qū)變形進(jìn)行了預(yù)測，結(jié)果表明，該模型可實(shí)現(xiàn)對礦區(qū)變形的快速監(jiān)測。

如何利用現(xiàn)有的貧信息、小樣本監(jiān)測數(shù)據(jù)[11]對礦區(qū)開展沉降預(yù)測是目前亟待解決的問題。本文以大紅山鐵礦為研究對象，利用MT-InSAR技術(shù)獲取礦區(qū)時間序列沉降量及變形速率，根據(jù)礦區(qū)地質(zhì)結(jié)構(gòu)和降雨量分析礦區(qū)地表變形機(jī)理，結(jié)合地面水準(zhǔn)監(jiān)測數(shù)據(jù)驗(yàn)證MT-InSAR技術(shù)在礦區(qū)變形監(jiān)測中的可行性和可靠性，并引入灰色沉降預(yù)測模型，以實(shí)現(xiàn)對小樣本數(shù)據(jù)等劣勢條件下的變形趨勢短期預(yù)測，為礦區(qū)安全生產(chǎn)提供保障，預(yù)防重大災(zāi)害的發(fā)生。

1 研究區(qū)概況及數(shù)據(jù)來源

1.1 研究區(qū)概況

大紅山鐵礦是我國主要的鐵礦石生產(chǎn)基地之一，位于云南省玉溪市新平縣城西119 km的新化鄉(xiāng)尺莫村，研究區(qū)地理位置與地形如圖1所示。礦區(qū)位于滇中臺坳南端以及紅河斷裂與綠汁江斷裂夾持的三角地帶，礦區(qū)西側(cè)出露有變質(zhì)較深、混合巖化較強(qiáng)的太古代哀牢山群，該區(qū)為北西、東西和南北向3組構(gòu)造線交會地帶。礦區(qū)屬亞熱帶氣候，干濕季分明，雨量較為充沛，年降雨量700～1 200 mm，多年平均降雨量930.8 mm，60%以上的雨量集中在6-9月。礦區(qū)目前有排土場、露天開采和地下開采區(qū)域，其中露天采場基本工作面高程在800～1 000 m，地下開采強(qiáng)度大且開采模式復(fù)雜，地表裂縫和滑坡不斷出現(xiàn)，傳統(tǒng)測量方法實(shí)施難度大、成本高。

圖1 研究區(qū)地理位置與地形

1.2 數(shù)據(jù)來源

采用歐洲航天局擁有12 d重訪周期的C波段傳感器的Sentinel-1A衛(wèi)星數(shù)據(jù)，選擇干涉寬幅模式下的斜距單視復(fù)數(shù)影像，極化方式為垂直極化。采用30 m分辨率的SRTM DEM描述地形。影像數(shù)據(jù)參數(shù)見表1。

表1 Sentinel-1A影像數(shù)據(jù)參數(shù)

2 研究方法及原理

2.1 SBAS-InSAR數(shù)據(jù)處理

InSAR技術(shù)可以提供基于時間序列雷達(dá)數(shù)據(jù)的開采沉陷信息和地表變形的時空演化信息[12]。SBAS-InSAR技術(shù)是通過獲取N景雷達(dá)數(shù)據(jù)設(shè)置短時間基線和空間基線閾值(本文為80%和10%)再將其組合，用Deluanary MCF最小費(fèi)用流解纏方法得到662對短基線差分干涉圖，根據(jù)相干系數(shù)選擇穩(wěn)定的高相干點(diǎn)目標(biāo)，利用奇異值分解法(SVD)[13]解算估計線性變形速率和高程誤差方程，再利用殘余相位估計大氣效應(yīng)等相位，最終通過地理編碼將相位轉(zhuǎn)變?yōu)樽冃螐亩@得相干點(diǎn)的變形速率及累計時間序列沉降。本文通過如圖2所示的SBAS-InSAR數(shù)據(jù)處理流程依次去除平地相位、地形相位、軌道誤差、大氣相位等，最終提取地面變形相位。通過獲取累計InSAR時間序列變形并將其作為灰色模型的輸入值，以實(shí)現(xiàn)對礦區(qū)變形趨勢的準(zhǔn)確預(yù)測。

圖2 數(shù)據(jù)處理流程

2.2 預(yù)測模型原理

灰色預(yù)測模型是通過對原始數(shù)據(jù)進(jìn)行累加處理生成一組灰色序列，減小其原始序列的隨機(jī)性，以此凸顯其指數(shù)增長規(guī)律，然后建立一階灰色微分方程，利用最小二乘法原理求取灰參數(shù)，最后經(jīng)累減還原得到預(yù)測值，其原理[14]敘述如下。

1)數(shù)據(jù)的檢驗(yàn)和處理

原始數(shù)據(jù)序列：

x(0)=[x(0)(1),x(0)(2),…,x(0)(n)]。

(1)

2)建立灰色模型

生成x(0)的累計序列：

x(1)=[x(1)(1),x(1)(2),…,x(1)(n)]

=[x(0)(1),x(0)(1)+x(0)(2),…,

x(0)(1)+…+x(0)(n)] 。

(2)

生成x(1)的均值序列：

z(1)(k)=0.5x(1)(k)+0.5x(1)(k-1),k=2,3,…,n。

(3)

建立白化微分方程：

(4)

建立灰色微分方程：

x(0)(k)+az(1)(k)=b。

(5)

最小二乘法計算發(fā)展系數(shù)a和灰作用量b，得到預(yù)測方程：

(6)

3 結(jié)果分析

3.1 精度驗(yàn)證

為驗(yàn)證MT-InSAR技術(shù)的監(jiān)測精度，將水準(zhǔn)點(diǎn)數(shù)據(jù)與InSAR數(shù)據(jù)作對比分析。InSAR測量結(jié)果與真實(shí)三維空間位移的函數(shù)關(guān)系[15]為

DLOS=DV·cosθ+DN·sinφsinθ-DE·cosφsinθ，

(7)

式中，DLOS表示沿衛(wèi)星雷達(dá)視線方向的位移；DV表示地面點(diǎn)沿垂直方向的位移；DN表示地面點(diǎn)沿南北方向的位移；DE表示地面點(diǎn)沿東西方向的位移；θ表示雷達(dá)衛(wèi)星入射角度；φ表示雷達(dá)衛(wèi)星飛行的方位角。

由于地表的水平位移不能忽略，且只有降軌數(shù)據(jù)，所以只能將獲取的12個地表觀測點(diǎn)實(shí)測水準(zhǔn)數(shù)據(jù)三向位移(2020年4月8日至2020年10月8日)代入式(7)獲得歸一化的D′，并與由MT-InSAR觀測得到的DLOS進(jìn)行對比，以驗(yàn)證InSAR監(jiān)測的可靠性，結(jié)果見表2、表3。

根據(jù)地面觀測點(diǎn)的實(shí)測三向數(shù)據(jù)(北向、東向、垂向)經(jīng)歸一化后得到的D′與InSAR取得的DLOS比較，求得絕對誤差和相對誤差(見表3)。由表3可知，兩者的絕對誤差值均小于2 mm，相對誤差在15%以內(nèi)。這表明MT-InSAR技術(shù)不僅能獲取地面變形情況，且其觀測精度比傳統(tǒng)測量技術(shù)的高。

表2 地面觀測點(diǎn)實(shí)測三向位移 單位：mm

表3 歸一化后的D′與DLOS的誤差

3.2 InSAR結(jié)果分析

通過獲取研究區(qū)138景(2014年10月至2020年10月)降軌數(shù)據(jù)，利用SBAS-InSAR技術(shù)獲得研究區(qū)的變形場，從而掌握礦區(qū)的變形空間分布和變形特征。研究發(fā)現(xiàn)礦區(qū)出現(xiàn)了3個明顯的沉降中心，即A、B、C區(qū)域，沉降面積分別為0.37、0.20、0.10 km2。其中，A沉降區(qū)域已經(jīng)發(fā)育成較成熟的沉降集中帶。另外，礦區(qū)還存在局部沉降區(qū)域D、E、F(見圖3)。A、B、E沉降區(qū)域主要與堆排土體有關(guān)，C、D沉降區(qū)域主要是地下開采所致，F(xiàn)沉降區(qū)域主要受露天開采影響。為進(jìn)一步探討礦區(qū)的變形機(jī)理和變形特征，分別在6個沉降區(qū)域選取9個特征點(diǎn)，并截取A、B、C沉降區(qū)域中心的剖面圖。

p-特征取樣點(diǎn)；JC-水準(zhǔn)點(diǎn)。

A、B區(qū)域是礦區(qū)的南部排土場，此處已經(jīng)發(fā)育形成2個典型的沉降集中帶，累計變形量在-555～-70 mm，平均變形速率在-90～-12 mm/a。對A區(qū)域經(jīng)過兩個沉降最大中心處沿a-a′的剖面(見圖4b)分析可知，這兩個沉降中心呈西北至東南方向的走勢，且該排土場的海拔在800～1 200 m，落差較大，沉降邊緣不斷向兩個方向擴(kuò)散，且存在明顯的不均勻沉降。

圖4 A、B區(qū)域特征點(diǎn)時序曲線及a-a′剖面圖

獲取p1、p2、p3、p4四個點(diǎn)的時間序列變形，研究發(fā)現(xiàn)4個取樣點(diǎn)呈近乎一致的沉降規(guī)律，沉降曲線逐漸減緩，表明變形速率不斷減?。辉诿磕甑?-8月雨季期間，沉降曲線會發(fā)生小幅波動，這是因?yàn)樵撆磐翀鲭m然不存在地下開采，但堆排土質(zhì)疏松(見圖5)，且土壤中存在大量積水，弱化了排土場的穩(wěn)定性，從而導(dǎo)致沉降的發(fā)生。

圖5 排土場現(xiàn)狀

C、D區(qū)域位于礦區(qū)地下開采的上覆地表，研究區(qū)變形場和變形速率圖顯示，該區(qū)域地表也形成了不同程度的沉陷，這兩處區(qū)域累計變形量在-348～-70 mm，平均變形速率在-55～12 mm/a。該區(qū)域長期存在地下開采，開采深度達(dá)1 200 m，規(guī)模為400萬t/a，地下高強(qiáng)度和大范圍的采動引起地表大面積沉降。分析C、D區(qū)域的p5、p6點(diǎn)的時間序列沉降(見圖6a)發(fā)現(xiàn)，雨季(6、7、8月)前后，兩處的變形曲線有明顯的波動，尤其是在2016年和2018年，且有加速下沉的趨勢，這可能與該時期的降雨量較大有關(guān)。分析沿兩個下沉中心的d-e剖面圖(見圖6b)發(fā)現(xiàn)，C、D區(qū)域依舊存在不均勻沉降，表明降雨也是造成地面塌陷的原因，且降雨量與地面下沉速率呈正相關(guān)，但此變形主要由地下開采所致。

圖6 p5、p6點(diǎn)時間序列沉降曲線及d-e剖面圖

由E區(qū)域的p7點(diǎn)的變形時間序列曲線(見圖7a)可知，該區(qū)域的變形速率逐漸減小，與實(shí)際情況相符。F區(qū)域是露天采場，邊坡上部坡度較大，由MT-InSAR獲取的變形(見圖7b)情況來看，下沉量較小，曲線呈周期性波動，這與露天邊坡長期的礦體開采有關(guān)。

圖7 p7、p和s點(diǎn)時間序列沉降曲線

3.3 預(yù)測分析

選取礦區(qū)水準(zhǔn)點(diǎn)作為建模數(shù)據(jù)，以CD01、JC04兩個實(shí)測點(diǎn)為研究對象，分別獲取12期的水準(zhǔn)監(jiān)測數(shù)據(jù)，提取同時間段SBAS-InSAR值，將其作為輸入值獲取沉降預(yù)測方程，繪制預(yù)測曲線，通過已獲取的小樣本數(shù)據(jù)預(yù)測礦區(qū)變形趨勢。將CD01、JC04點(diǎn)的第10、11、12期的預(yù)測值分別與InSAR值和水準(zhǔn)數(shù)據(jù)對比，以驗(yàn)證模型的預(yù)測精度和可靠性，結(jié)果見圖8、表4。

圖8 CD01、JC04點(diǎn)的沉降曲線對比

表4 CD01、JC04點(diǎn)的沉降量相對誤差

由圖8、表4可知：監(jiān)測點(diǎn)CD01、JC04的灰色模型預(yù)測曲線與水準(zhǔn)及InSAR曲線較為吻合，模型擬合效果較好；CD01、JC04點(diǎn)在第10、第11期的監(jiān)測相對誤差在5%以內(nèi)、預(yù)測相對誤差在9%以內(nèi)，而第12期的監(jiān)測相對誤差在3%以內(nèi)、預(yù)測相對誤差在12%以內(nèi)。但隨著監(jiān)測周期的拉長，CD01、JC04點(diǎn)的監(jiān)測相對誤差和預(yù)測相對誤差會增大，這與文獻(xiàn)[16]的結(jié)論一致。該灰色模型計算簡單快捷。

4 結(jié)論

a.利用時序SBAS-InSAR技術(shù)獲取研究區(qū)的變形場，研究發(fā)現(xiàn)礦區(qū)內(nèi)存在4個發(fā)育成熟的沉降集中帶A、B、C、D及局部沉降區(qū)域E、F。區(qū)域A和B主要是由排土場局部的不穩(wěn)定和雨季降水綜合影響所致。區(qū)域C和D主要是由長期頻繁的地下開采及雨季強(qiáng)降雨綜合影響所致。區(qū)域E是廢棄的排土場，其變形情況與實(shí)際相符。區(qū)域F是由南向北開采的露天采場，沉降量較小，曲線呈周期性波動，這與露天邊坡長期的礦體開采有關(guān)。

b.比較地面觀測點(diǎn)實(shí)測三向位移經(jīng)歸一化后得到的位移與InSAR取得的LOS向位移后發(fā)現(xiàn)，絕大部分點(diǎn)的相對誤差在15%以內(nèi)，表明MT-InSAR在觀測精度上是可靠的。

c.灰色模型獲得的預(yù)測曲線與水準(zhǔn)及InSAR曲線較為吻合，預(yù)測模型的相對誤差在12%以內(nèi)。