基于多軌道SBAS方法的露天礦高陡邊坡形變監(jiān)測

魏戀歡, 封秋月, 毛亞純, 劉善軍

(東北大學(xué) 資源與土木工程學(xué)院, 遼寧 沈陽 110819)

我國冶金礦山80%的礦石來源于露天開采,隨著開采深度的不斷增加,我國大多數(shù)露天礦已進(jìn)入深凹開采模式,導(dǎo)致礦坑邊坡高度及邊坡角不斷增大,邊坡穩(wěn)定性隨之降低.對于深凹露天礦而言,高陡邊坡的穩(wěn)定性是影響礦山安全生產(chǎn)的重要因素.露天礦邊坡滑移災(zāi)害的發(fā)生表面上是突發(fā)的,但其過程是逐步醞釀和發(fā)展的,在災(zāi)害發(fā)生前會存在長期的緩慢形變積累.因此,采取有效手段實時監(jiān)測高陡邊坡的穩(wěn)定性，研究其形變特點(diǎn)及成因、及時治理潛在的險情，是礦山安全生產(chǎn)的重要保障.

目前,露天礦高陡邊坡形變監(jiān)測一般采用全站儀或GPS等現(xiàn)場測量手段,存在效率低、人力成本高、僅能獲得少量點(diǎn)位形變、無法提取大范圍分布特征等不足.近年來,隨著合成孔徑雷達(dá)(synthetic aperture radar,SAR)衛(wèi)星與數(shù)據(jù)處理技術(shù)的不斷發(fā)展,時序InSAR(time series interferometry SAR)作為一種新型的形變測量手段,具有低成本、高精度等優(yōu)點(diǎn),已成功應(yīng)用于沉降、滑坡、地震、火山活動等地表形變監(jiān)測中[1].時序InSAR從一組時間序列SAR影像中提取不受時空去相干影響的目標(biāo)像元,每一個像元的干涉相位φ可分解為如下幾個部分:φflat，φtopo，φdef，φnoise.

平地效應(yīng)相位φflat和地形相位φtopo是由于SAR衛(wèi)星多次成像時所處位置不同、地球橢球和地表高程變化造成的;φdef為形變相位;φnoise主要是由軌道誤差、大氣效應(yīng)等造成的噪聲相位[1].針對這些像元開展相位分析,可達(dá)到分離各種分量并解算精細(xì)形變的目的,實現(xiàn)毫米級的形變監(jiān)測精度.目前已有的時序InSAR算法主要分為基于點(diǎn)目標(biāo)的永久散射體法(permanent scatterers InSAR,PS-InSAR)和基于分布式目標(biāo)的小基線子集法(small baseline subset,SBAS)兩大類[2-7].相比于PS-InSAR,SBAS可以有效利用所有時空基線較短的高質(zhì)量差分干涉圖,提高了數(shù)據(jù)利用率和點(diǎn)密度,更加適合進(jìn)行長時間的形變監(jiān)測.

受限于SAR衛(wèi)星固有的側(cè)視成像特點(diǎn),時序InSAR得出的是目標(biāo)像元沿衛(wèi)星視線(line of sight,LOS)方向的一維形變信息,是真實形變在LOS方向的投影,因此需要由LOS向形變求解真實形變.在大范圍沉降監(jiān)測中,通常認(rèn)為地物只存在垂直形變,不存在水平形變,簡單地用LOS方向的形變結(jié)果除以影像入射角的余弦值即可得到垂直形變結(jié)果[8-9].然而,在露天礦高陡邊坡形變監(jiān)測中,由于礦坑邊坡滑移通常都是沿坡面梯度方向的,而不同空間位置的坡度和坡向差別很大,無法直接根據(jù)入射角求解[10].

本文針對上述問題,介紹了SBAS法的基本原理,面向深凹露天礦高陡邊坡形變監(jiān)測需求,提出了一種基于地形的坡向形變求解方法,利用SBAS法分析了覆蓋鞍鋼集團(tuán)大孤山鐵礦的3組Sentinel-1數(shù)據(jù),提取了3組LOS方向形變參數(shù),將3組LOS方向形變分別轉(zhuǎn)換到坡向,與同時期的測量機(jī)器人監(jiān)測數(shù)據(jù)及降水量進(jìn)行了對比分析,揭示了大孤山鐵礦礦坑邊坡的穩(wěn)定性及影響因素.

1 小基線子集法

作為時序InSAR方法中的一個典型代表,SBAS是在PS-InSAR的基礎(chǔ)上發(fā)展而來.PS-InSAR從一組SAR影像集中選取一個公共主影像,將所有從影像與主影像干涉得到時序干涉圖集.公共主影像的使用使一些時空基線均較大的低質(zhì)量干涉圖參與了PS點(diǎn)目標(biāo)提取與形變解算,導(dǎo)致人工地物稀少的區(qū)域點(diǎn)目標(biāo)密度偏低,難以精確估算形變分布信息[2-3].SBAS針對這一問題,采用了多主影像干涉模式,根據(jù)時空基線分布特征將SAR數(shù)據(jù)集劃分為多個子集,每個子集內(nèi)時空基線均較短,各子集間通過少量長時空基線的干涉對相連,可以保證數(shù)據(jù)的高時間采樣率和高點(diǎn)密度,降低時空去相干的影響.SBAS的多主影像干涉模式增加了干涉圖數(shù)量,避免了因干涉圖數(shù)量少導(dǎo)致的方程秩虧問題,短基線干涉組合可以降低數(shù)字高程模型(digital elevation model,DEM)誤差對形變量測精度的影響,有利于直接解算研究區(qū)域的形變[6].因此,本文采用SBAS法對大孤山鐵礦的高陡邊坡進(jìn)行時序InSAR分析.

在SBAS法中,首先根據(jù)數(shù)據(jù)集的時空基線分布特征選擇超級主影像,將所有從影像都以超級主影像為參考進(jìn)行粗配準(zhǔn)和精配準(zhǔn).其次,根據(jù)多主影像策略生成小基線差分干涉圖集.再次,因為干涉圖中獲得的相位是真實相位對2π的余數(shù)值,其取值范圍是在-π到π之間的，稱為“纏繞相位”，不能直接轉(zhuǎn)換為形變值。通過最小費(fèi)用流算法對序列差分干涉圖進(jìn)行解纏，將“纏繞相位”加上2π的整數(shù)倍，即可恢復(fù)為“真實相位”，稱為“解纏相位”。用奇異值分解(singular value decomposition,SVD)法將解纏后的差分干涉圖結(jié)合起來生成一個包含大氣、地形和形變信號的時間序列.然后,通過二次解纏和濾波消除地形、大氣信號,最終得到LOS方向形變速率和形變量[4,11].

2 邊坡坡向形變求解方法

以圖1a中的A點(diǎn)為例，假設(shè)其沿坡度方向的位移為dA，可按如下公式分解為水平分量和垂直分量:

(1)

式中：β為坡度；dA,v與dA,h分別為dA的垂直與水平分量。SBAS提取的LOS方向形變值dA,los實際上是dA,v和dA,h在LOS向的投影之和。其中dA,v在LOS向的投影可直接根據(jù)傳感器的入射角θ計算。為了求解dA,h在LOS向的投影，需要先將dA,h投影到LOS向所在的鉛垂面[10].由于升軌數(shù)據(jù)和降軌數(shù)據(jù)的LOS向不同,需要分別討論.這里“升軌”是指衛(wèi)星由南向北飛行,“降軌”是指衛(wèi)星由北向南飛行.受到地球自轉(zhuǎn)的影響,升軌時衛(wèi)星實際前進(jìn)方向為北偏西,相應(yīng)的LOS方向在水平面的投影為東偏北；降軌時衛(wèi)星實際前進(jìn)方向為南偏西,相應(yīng)的LOS方向在水平面的投影為西偏北.

根據(jù)降軌數(shù)據(jù)的幾何關(guān)系(圖1a)可知，dA,h在LOS所在鉛垂面的投影dA,h-los可表示為

dA,h-los=dA,h·sin(αA-α0).

(2)

其中：αA是dA,h與北向的夾角；α0是雷達(dá)的航向角。設(shè)傳感器入射角為θ，A點(diǎn)的LOS向形變dA,los由dA,v和dA,h-los在LOS向投影的矢量和組成，可表示為

dA,los=dA,h-los·sin(θ)-dA,v·cos(θ).

(3)

將式(1)、式(2)代入式(3)可得：

dA,los=dA·cos(β)·sin(αA-α0)·sin(θ)-dA·sin(β)·cos(θ).

(4)

因此，在降軌數(shù)據(jù)條件下，A點(diǎn)的LOS向形變可轉(zhuǎn)換為坡向形變：

(5)

類似地，根據(jù)升軌數(shù)據(jù)的幾何成像(圖1b)可知，dA,h在LOS所在鉛垂面的投影dA,h-los可表示為

dA,h-los=dA,h·sin(αA-α0-π)=-dA,h·sin(αA-α0).

(6)

A點(diǎn)的LOS向形變dA,los由dA,v和dA,h-los在LOS向投影的矢量和組成，表示為

dA,los=-dA,h-los·sin(θ)-dA,v·cos(θ).

(7)

將式(1)、式(6)代入式(7)后可以得到dA和dA,los的關(guān)系，與式(4)相同。因此，在升軌數(shù)據(jù)條件下，A點(diǎn)的LOS向形變轉(zhuǎn)為坡向形變的公式與式(5)相同[10]。在本文實驗中，三組數(shù)據(jù)的航向角(以北向起始順時針為正)α0和入射角θ如表1所示。由于大孤山鐵礦礦坑西北幫巖體破碎，工程地質(zhì)條件復(fù)雜，為了保證安全生產(chǎn)，開采坡度角β約為37°，坡度方向與北向的水平夾角與點(diǎn)目標(biāo)所在位置有關(guān).

表1 時序SAR數(shù)據(jù)信息

3 研究區(qū)及數(shù)據(jù)概況

3.1 研究區(qū)概況

大孤山鐵礦位于我國遼寧省鞍山市,是鞍山鋼鐵集團(tuán)鞍山礦業(yè)公司5大露天鐵礦之一,是我國典型的短深型深凹露天礦,也是亞洲最深的露天鐵礦之一[12].大孤山鐵礦于1970年轉(zhuǎn)入深凹露天縱向開采,形成了東西長1 700 m、南北寬1 520 m、垂直深度達(dá)348—398 m左右的橢圓形露天采坑[12].礦區(qū)年降雨量平均為720 mm,一年中以6—8月降雨量最多,合計438 mm,9—11月降雨量次之,為151 mm.

如圖2所示,大孤山鐵礦的采坑(黃框區(qū)域)位于礦區(qū)西部,開采技術(shù)為高陡邊坡開采,露采工作呈階梯形向地下延伸,其東側(cè)及東北側(cè)的“倒L”形區(qū)域為排土場.該排土場大部分區(qū)域為巖石松散裸露堆積.排土場的東南側(cè)是大球場尾礦庫.大球場尾礦庫屬于山谷型尾礦庫,溝口向北,其余方向均為山脊,尾礦庫壩的筑壩方式是排土場渣土分層碾壓筑壩.對大孤山鐵礦的高陡邊坡的穩(wěn)定性進(jìn)行持續(xù)監(jiān)測至關(guān)重要,是確保礦山安全生產(chǎn)的關(guān)鍵.

3.2 數(shù)據(jù)概況

本文收集了覆蓋大孤山鐵礦的3組Sentinel-1數(shù)據(jù),空間覆蓋范圍如圖3所示,詳細(xì)參數(shù)如表1所示.本文采用的DEM數(shù)據(jù)為德國TanDEM-X計劃獲取的DEM,空間分辨率為3 rad·s,絕對高程精度約1 m[13-14].由于SBAS法使用的差分干涉圖時空基線都相對較短,DEM誤差對差分干涉質(zhì)量的影響相對較小.

圖4中展示了3組Sentinel-1數(shù)據(jù)的獲取時間以及SBAS法生成的差分干涉組合.其中,升軌Sentinel-1A數(shù)據(jù)生成了145幅差分干涉圖；軌道編號為P105的降軌Sentinel-1B數(shù)據(jù)生成了250幅差分干涉圖；軌道編號為P3的降軌Sentinel-1B數(shù)據(jù)生成了229幅差分干涉圖.

4 結(jié)果分析

4.1 監(jiān)測結(jié)果分析

SBAS法提取的大孤山鐵礦采坑區(qū)域的LOS向平均形變速率如圖5所示,其中紅色色標(biāo)表示地物目標(biāo)向遠(yuǎn)離傳感器方向移動,藍(lán)色色標(biāo)表示地物目標(biāo)向靠近傳感器方向移動；黃色箭頭指示的Az方向表示衛(wèi)星飛行方向,LOS方向表示衛(wèi)星視線方向.交叉驗證表明,三組Sentinel-1結(jié)果具有較高一致性,監(jiān)測結(jié)果較可靠.大孤山鐵礦高陡邊坡的主要形變區(qū)域位于采坑西北幫,已在圖5中用綠色多邊形標(biāo)出,該區(qū)域的大部分地物目標(biāo)LOS向形變速率均超過了-60 mm/年,最高達(dá)到了-221 mm/年.

采坑西北幫的形變可大致分為A，B，C三個小區(qū)域.其中,A,B區(qū)域在升軌結(jié)果中為遠(yuǎn)離傳感器方向形變,在兩組降軌結(jié)果中均為朝向傳感器方向形變.由于邊坡的真實形變一般是沿坡向朝向坑底的,SBAS結(jié)果僅為真實形變在LOS向的投影.考慮到升降軌數(shù)據(jù)的LOS方向不同,升降軌結(jié)果存在差異是合理的.C區(qū)域在三組結(jié)果中都顯示為遠(yuǎn)離傳感器方向形變,主要是由于該區(qū)域位于礦坑邊緣部位、坡度相對較小且坡向接近南向,因此升降軌結(jié)果中的差異不明顯.

圖6中的Google Earth歷史影像顯示,2018年9月10日前A區(qū)域發(fā)生了大規(guī)模的巖體移動,巖體的移動方向為沿坡面指向坑底方向,符合露天開挖導(dǎo)致的滑移特征.根據(jù)本文提出的方法,求解了位于A區(qū)域滑坡體上緣P1點(diǎn)的坡向形變.圖7展示了P1點(diǎn)的坡向形變量曲線,為了方便對比,三組數(shù)據(jù)都以2017年6月5日為參考,截止2019年6月初累計坡向滑移量均超過了-600 mm.

大多數(shù)露天礦的滑坡均起因于降水引起的沿軟弱結(jié)構(gòu)面的滑動,每年的雨季都是滑坡災(zāi)害發(fā)生最頻繁的時期.降水對邊坡穩(wěn)定性的影響較大,水的浸泡、滲透作用可以降低軟弱面的抗剪強(qiáng)度,增加下滑力,從而導(dǎo)致邊坡破壞[15].圖7顯示,2017年和2018年夏季持續(xù)大量降水之后,P1點(diǎn)在三組結(jié)果中均表現(xiàn)出了加速形變,入冬后逐漸減速.不過,總體而言P1點(diǎn)的加速形變不大,未出現(xiàn)嚴(yán)重的邊坡失穩(wěn)問題.

在現(xiàn)場布設(shè)了2個監(jiān)測點(diǎn)(圖6中P2，P3),位于B區(qū)域,于2017年8月13日至2018年11月1日期間使用測量機(jī)器人進(jìn)行了連續(xù)監(jiān)測.圖8所示的是P2和P3點(diǎn)的坡向形變時間序列,為了方便對比分析,3組SBAS結(jié)果都以測量機(jī)器人第一次觀測時刻(2017年8月13日)為參考.監(jiān)測結(jié)果表明,3組SBAS結(jié)果與測量機(jī)器人監(jiān)測結(jié)果存在高度一致性,P2點(diǎn)和P3點(diǎn)形變規(guī)律相似,2017年8月13日至2019年6月初累計坡向形變均超過了-500 mm,均在每年夏季大量降水后產(chǎn)生了一定的加速形變,入冬時均有減速,未發(fā)生明顯失穩(wěn)問題.

P2點(diǎn)和P3點(diǎn)的SBAS監(jiān)測結(jié)果與測量機(jī)器人監(jiān)測結(jié)果之間的均方根誤差如表2所示,其中S1B_P105數(shù)據(jù)的均方根誤差最小,分別為17.4 mm和19.2 mm.時間序列結(jié)果表明,持續(xù)大量降水會導(dǎo)致邊坡形變量和形變速率的增大,這也是在雨季提高測量機(jī)器人的監(jiān)測頻率、縮短監(jiān)測周期的主要原因.為了保障安全生產(chǎn),應(yīng)在每年雨季對西北幫高陡邊坡實施穩(wěn)定性監(jiān)測.

表2 P2，P3點(diǎn)均方根誤差

4.2 工程地質(zhì)條件分析

工程地質(zhì)資料(圖9)顯示,大孤山鐵礦露天采坑的東端幫與南幫以花崗巖為主,巖性及構(gòu)造條件良好,總體較為穩(wěn)定.西幫(包括西北幫和西南幫)可分為典型的4個區(qū)域[16].西南幫為混合巖區(qū)域,也是典型的塊裂體工程地質(zhì)模型,但節(jié)理發(fā)育情況好于西北幫,且斷層較少,穩(wěn)定性較好.西北幫的巖體較為破碎,工程地質(zhì)條件復(fù)雜,其穩(wěn)定性直接影響礦山的安全生產(chǎn).

西北幫可大致分為A,B,C三個區(qū)域.A區(qū)域為品位較低的楔形礦體,其西南側(cè)以F15斷層為界,東北側(cè)以F14斷層與B區(qū)域接觸,F14斷層周邊的破碎帶穩(wěn)定條件較差,威脅著采場北幫的鐵路運(yùn)輸；B區(qū)域的巖體以混合巖為主,其穩(wěn)定性受幾條較大的斷層控制,屬于典型的塊裂體工程地質(zhì)模型[17]；C區(qū)域的巖體以千枚巖為主,為典型的碎裂和散體工程地質(zhì)模型,節(jié)理裂隙非常發(fā)育,巖石強(qiáng)度較低,在雨水的作用下容易發(fā)生滑坡現(xiàn)象[17].雖然B,C區(qū)域在Google Earth歷史影像中并未出現(xiàn)明顯的邊坡滑動,但從巖土結(jié)構(gòu)上看,易發(fā)生形變破壞的地段多為千枚巖地段和混合巖地段(B,C區(qū)域),這主要是由于該類巖體整體完整度低,堅固程度差所造成的.考慮到C區(qū)域位于礦坑邊緣且所占范圍逐漸減小,盡管千枚巖的抗壓、抗剪強(qiáng)度均較低,易產(chǎn)生形變破壞,其穩(wěn)定性對礦山安全生產(chǎn)的威脅不大[17].因此,B區(qū)域是威脅大孤山鐵礦安全生產(chǎn)的潛在滑坡區(qū),應(yīng)持續(xù)進(jìn)行形變監(jiān)測.

5 結(jié) 論

1) 利用SBAS法處理了多軌道時序SAR數(shù)據(jù),對大孤山露天礦的高陡邊坡進(jìn)行了穩(wěn)定性監(jiān)測.

2) 提出了一種基于地形特征的LOS向形變到坡向形變的解算方法,得到了采坑西北幫的坡向形變結(jié)果,與實測數(shù)據(jù)存在高度一致性.

3) 通過研究同期的降水?dāng)?shù)據(jù)及大孤山工程地質(zhì)條件,對變形原因進(jìn)行分析,表明露天礦的穩(wěn)定性主要受巖體結(jié)構(gòu)、巖性與降水的影響.