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

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

0 引言

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

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

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

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

1.1 研究區(qū)概況

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

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

1.2 數(shù)據(jù)來(lái)源

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

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

2 研究方法及原理

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

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

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

2.2 預(yù)測(cè)模型原理

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

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

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

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

(1)

2)建立灰色模型

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

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)

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

(6)

3 結(jié)果分析

3.1 精度驗(yàn)證

為驗(yàn)證MT-InSAR技術(shù)的監(jiān)測(cè)精度，將水準(zhǔn)點(diǎn)數(shù)據(jù)與InSAR數(shù)據(jù)作對(duì)比分析。InSAR測(cè)量結(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個(gè)地表觀測(cè)點(diǎn)實(shí)測(cè)水準(zhǔn)數(shù)據(jù)三向位移(2020年4月8日至2020年10月8日)代入式(7)獲得歸一化的D′，并與由MT-InSAR觀測(cè)得到的DLOS進(jìn)行對(duì)比，以驗(yàn)證InSAR監(jiān)測(cè)的可靠性，結(jié)果見(jiàn)表2、表3。

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

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

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

3.2 InSAR結(jié)果分析

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

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

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

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

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

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

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

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

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

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

3.3 預(yù)測(cè)分析

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

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

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

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

4 結(jié)論

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

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

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