基于時(shí)序InSAR的改進(jìn)型PIM礦山采空區(qū)沉降模型構(gòu)建與應(yīng)用

毛嘉騏,李素敏,成 睿,張龍宇,沈顯名

(昆明理工大學(xué) 國(guó)土資源工程學(xué)院,云南 昆明 650093)

0 引言

地下礦體開(kāi)采后易引發(fā)地表塌陷等地質(zhì)災(zāi)害,獲取沉降區(qū)具體位置對(duì)減少災(zāi)害事故至關(guān)重要[1]。傳統(tǒng)的獲取方式以物探為主,雖然精度較高,但工作范圍小,且需投入大量人力物力[2]。合成孔徑雷達(dá)干涉測(cè)量(Interferometric Synthetic Aperture Radar,InSAR)具有全天時(shí)、全天候和覆蓋范圍廣等優(yōu)點(diǎn),為采空區(qū)地表沉降監(jiān)測(cè)提供了有力的技術(shù)支持[3]。

利用InSAR對(duì)采空區(qū)地表進(jìn)行監(jiān)測(cè)時(shí),由于礦山開(kāi)采引起的地表沉降速度快、形變量大,且SAR衛(wèi)星存在一定的重復(fù)周期,導(dǎo)致沉降區(qū)中心附近出現(xiàn)低相干、失相干現(xiàn)象,形變結(jié)果缺失,無(wú)法獲取準(zhǔn)確的沉降區(qū)模型[4]。為解決上述問(wèn)題,當(dāng)前采用的主要方法是概率積分法(PIM)[5],PIM是由劉寶琛基于隨機(jī)介質(zhì)理論而提出的,其主要應(yīng)用于煤礦規(guī)則采空區(qū)。PIM利用高精度監(jiān)測(cè)點(diǎn)對(duì)相關(guān)參數(shù)進(jìn)行反演,以解決采空區(qū)地表沉降中心相干性較低的問(wèn)題。為得到準(zhǔn)確度更高的PIM參數(shù),目前應(yīng)用最廣泛的參數(shù)優(yōu)化反演方法有最小二乘法[6]、遺傳算法[7]、經(jīng)驗(yàn)設(shè)計(jì)值法等。然而不同的反演算法均有其固有缺陷,最小二乘法和經(jīng)驗(yàn)設(shè)計(jì)值法的抗粗差能力較弱,反演結(jié)果波動(dòng)性較大;遺傳算法全局搜索能力不足,易導(dǎo)致局部最優(yōu)解而使收斂速度過(guò)快,其在采空區(qū)沉降監(jiān)測(cè)中的精度較差[8]。

本文提出了一種將時(shí)序InSAR與PIM相融合的方法,旨在獲取準(zhǔn)確的采空區(qū)地表沉降信息。該方法首先利用時(shí)序InSAR技術(shù)獲取礦區(qū)工作面沉降結(jié)果,并對(duì)相干性進(jìn)行分析;之后選取InSAR沉降區(qū)邊界相干性大于0.3的點(diǎn)以及沉降中心線性插值后的少量水準(zhǔn)數(shù)據(jù)作為特征點(diǎn),通過(guò)改進(jìn)型PSO算法對(duì)PIM參數(shù)進(jìn)行反演,以跳出局部最優(yōu)解,增強(qiáng)其抗粗差能力,得到可靠的PIM預(yù)計(jì)參數(shù)以及沉降區(qū);提取時(shí)序InSAR沉降區(qū)邊界形變值,與由PIM得到的沉降中心進(jìn)行插值融合,完成工作面地表下沉的信息提取,在采用少量地表監(jiān)測(cè)數(shù)據(jù)的情況下構(gòu)建完整的地表沉降模型。

1 研究區(qū)概況及前期研究存在的問(wèn)題

1.1 研究區(qū)概況

以云南省玉溪市大紅山鐵礦地下采區(qū)某工作面為研究對(duì)象,其地理位置和范圍見(jiàn)圖1。礦區(qū)地屬滇中臺(tái)南坳南端以及紅河斷裂與綠汁江斷裂夾持的三角地帶,西側(cè)出露有變質(zhì)較深、混合巖化較強(qiáng)的太古代哀牢山群。圖1c表示該礦區(qū)某地下開(kāi)采工作面的某巷道開(kāi)采位置,該工作面對(duì)應(yīng)的地表植被覆蓋情況一般[9]。工作面礦層平均厚度72.58 m,主要采用無(wú)底柱分段崩落法開(kāi)采,礦體的分布標(biāo)高為25.72～945.00 m。該工作面走向開(kāi)采長(zhǎng)度為719 m,傾向開(kāi)采長(zhǎng)度為505 m,平均采深500 m,礦層傾角10°,開(kāi)采厚度5 m。圖1c中:黑色線表示工作面上方走向線,走向線自東南向西北方向布設(shè)有H1-H5點(diǎn)位;白色線表示傾向線,自西南向東北方向布設(shè)有J1、J2、J3點(diǎn)位。測(cè)量工作自2021年6月15日開(kāi)始,至2022年7月18日結(jié)束,共7期。

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

1.2 前期研究存在的問(wèn)題

前期采用時(shí)序InSAR方法對(duì)地表形變進(jìn)行了監(jiān)測(cè),其原理是利用同一地區(qū)的多幅時(shí)間基線較短的SAR影像形成的干涉對(duì),通過(guò)解纏、濾波等方式去除軌道誤差、噪音以及地形的殘余相位,然后采用奇異值分解(SVD)的方法,將多個(gè)基線集聯(lián)合求解,并對(duì)時(shí)間域和空間域的濾波進(jìn)行分析,分離出殘余相位中的大氣相位和非線性形變誤差,得到目標(biāo)區(qū)內(nèi)覆蓋整個(gè)觀測(cè)時(shí)間的地表形變信息[10]。本文選取IW模式下的2021年6月10日至2022年7月23日的34景VV極化方式的Sentinel-1A降軌數(shù)據(jù),重訪周期為12 d,波長(zhǎng)5.6 cm,中心入射角39.5°,使用NASA提供的30 m分辨率SRTM DEM去除干涉計(jì)算地形相位對(duì)結(jié)果的影響。

干涉對(duì)參數(shù)信息見(jiàn)表1。設(shè)置最大臨界空間基線占比10%,最大時(shí)間基線為60 d,按照基線配對(duì)的方式兩兩差分,按照相干系數(shù)法設(shè)置0.3為平均相干系數(shù)及最小相干系數(shù)閾值選取高相干點(diǎn),之后利用自適應(yīng)濾波的方法去除噪聲對(duì)干涉相位的影響,得到的沉降結(jié)果見(jiàn)圖2。

表1 干涉對(duì)參數(shù)信息

圖2 時(shí)序InSAR視線向沉降結(jié)果

由圖2可知:研究區(qū)出現(xiàn)了明顯的失相干現(xiàn)象,導(dǎo)致形變結(jié)果缺失;邊緣部分的形變結(jié)果顯示,2022年1月24日采區(qū)沉降面初步形成,說(shuō)明此時(shí)該工作面的開(kāi)采活動(dòng)對(duì)上覆巖層穩(wěn)定性產(chǎn)生了影響,但由于該工作面開(kāi)采時(shí)間較短,沉降量較小,截至2022年7月23日,利用時(shí)序InSAR監(jiān)測(cè)到的最大累積沉降量為48 mm。隨著開(kāi)采沿東南向西北方向移動(dòng),采區(qū)地表逐漸形成了較大的沉降區(qū)。

由研究區(qū)相干系數(shù)圖(見(jiàn)圖3)可知,在靠近采空區(qū)地表中心附近的少部分區(qū)域相干系數(shù)低于0.3,相干性較差[11],除失相干區(qū)域外,其余大部分區(qū)域相干性良好。結(jié)合InSAR監(jiān)測(cè)結(jié)果對(duì)實(shí)地進(jìn)行了踏勘,發(fā)現(xiàn)地表出現(xiàn)了明顯的局部沉降現(xiàn)象(見(jiàn)圖4a),且周邊存在較大的貫穿裂縫(見(jiàn)圖4b、圖4c),說(shuō)明該區(qū)域地表產(chǎn)生了較大形變。

圖3 研究區(qū)相干系數(shù)圖

圖4 采空區(qū)地表現(xiàn)場(chǎng)踏勘

2 研究方法

2.1 PIM模型

PIM是一種非連續(xù)的隨機(jī)介質(zhì)模型,由于巖層在移動(dòng)過(guò)程中會(huì)破壞介質(zhì)的連續(xù)性,介質(zhì)單元發(fā)生變化,導(dǎo)致沉積環(huán)境也發(fā)生變化,從而產(chǎn)生不同的新巖層,在地質(zhì)構(gòu)造作用下,新巖層中形成節(jié)理、斷層等弱面,弱面即為非連續(xù)的隨機(jī)介質(zhì)[11-12]。礦層開(kāi)采后巖層與地表移動(dòng)是一種隨機(jī)的過(guò)程,形成的地表塌陷區(qū)見(jiàn)圖5。

圖5 沉降區(qū)在某一點(diǎn)處下沉示意圖

由PIM給出任意點(diǎn)T(x,y)的下沉預(yù)計(jì)模型:

(1)

(2)

(3)

式中,Wi為最大沉降值,Wi(x)、Wi(y)分別為某一點(diǎn)在主斷面上的走向、傾向沉降量,u為PIM開(kāi)采單元,H為工作面采深,m為礦層厚度,q為下沉系數(shù),α為礦層傾角,θ為開(kāi)采影響傳播角,tanβ為主要影響角的正切值,D1、D3分別為工作面的走向、傾向長(zhǎng),s1、s2分別為走向上的左右邊界拐點(diǎn)偏移距,s3、s4分別為上山下山方向的拐點(diǎn)偏移距[11]。

2.2 改進(jìn)型PSO算法

PSO(Particle Swarm Optimization)算法[13]是對(duì)鳥(niǎo)類集體覓食行為的模仿,不同個(gè)體之間通過(guò)集體協(xié)作以及信息共享,使種群的移動(dòng)方向保持一致,而個(gè)體間保持一定距離來(lái)尋找最優(yōu)位置,是一種基于粒子群的優(yōu)化算法。經(jīng)典的粒子群算法隨著種群進(jìn)化會(huì)使其多樣性呈遞減趨勢(shì),從而過(guò)早地快速收斂到局部最優(yōu)解。本文主要以較大的慣性權(quán)重和學(xué)習(xí)因子為開(kāi)局,有利于增強(qiáng)初始粒子群的全局迭代尋優(yōu)能力,最后使用較小的慣性權(quán)重和學(xué)習(xí)因子結(jié)束,以使粒子群跳出局部最優(yōu)解。

在對(duì)種群進(jìn)行初始化處理的基礎(chǔ)上,選擇合適的參數(shù),迭代計(jì)算每個(gè)粒子的適應(yīng)值,將個(gè)體極值和種群極值代入式(4)、式(5),不斷更新每個(gè)粒子的位置和移動(dòng)方向,最終得到所需的尋優(yōu)解[13]。

Vi=ω×Vi+c1×rand()×(pbesti-Xi)
+c2×rand()×(gbesti-Xi),

(4)

Xi=Xi+Vi,

(5)

式中,Vi為第i個(gè)粒子的當(dāng)前速度,Xi為第i個(gè)粒子的當(dāng)前位置,c1、c2為學(xué)習(xí)因子,rand()是(0,1)之間的隨機(jī)數(shù),ω為整體慣性權(quán)重。

2.3 工作流程

a.對(duì)覆蓋礦區(qū)的SAR影像進(jìn)行時(shí)序InSAR處理,獲取礦區(qū)的時(shí)序形變信息。收集礦區(qū)與SAR影像同時(shí)段內(nèi)的走向和傾向線上的水準(zhǔn)數(shù)據(jù);

b.在工作面坐標(biāo)系下,選取時(shí)序InSAR沉降區(qū)邊緣附近高相干點(diǎn)和水準(zhǔn)數(shù)據(jù)作為特征點(diǎn),利用改進(jìn)型PSO算法對(duì)PIM參數(shù)進(jìn)行尋優(yōu)反演;

c.根據(jù)改進(jìn)型PSO算法反演出的PIM參數(shù),建立礦區(qū)的PIM沉降模型,并將結(jié)果統(tǒng)一至WGS-84坐標(biāo)系,利用克里金插值法,將InSAR沉降區(qū)邊緣形變值與由PIM獲得的中心形變值結(jié)合,獲取連續(xù)且完整的沉降結(jié)果。具體的工作流程見(jiàn)圖6。

圖6 工作流程

3 基于融合方法的礦區(qū)沉降結(jié)果提取及分析

3.1 基于時(shí)序InSAR的特征點(diǎn)選取

為構(gòu)建概率積分沉降模型,需先將InSAR的LOS向形變投影至垂直方向,之后提取InSAR沉降結(jié)果中大于等于0.5的點(diǎn)位進(jìn)行反演分析[14]。選取InSAR沉降盆地邊緣的30個(gè)相干性較高的點(diǎn)位,并對(duì)該工作面上方的走向和傾向線的7個(gè)水準(zhǔn)點(diǎn)進(jìn)行線性插值,使InSAR沉降結(jié)果與水準(zhǔn)數(shù)據(jù)在時(shí)間上統(tǒng)一,共選取37個(gè)特征點(diǎn),特征點(diǎn)位置見(jiàn)圖7(白色圓點(diǎn)表示InSAR高相干性點(diǎn),黑色三角點(diǎn)表示走向和傾向線水準(zhǔn)監(jiān)測(cè)點(diǎn),黃色圓點(diǎn)表示邊界GPS點(diǎn),黃色三角點(diǎn)表示其余水準(zhǔn)點(diǎn))。

圖7 特征點(diǎn)分布

3.2 基于改進(jìn)型PSO算法的PIM參數(shù)反演

PIM預(yù)計(jì)的參數(shù)分別為:下沉系數(shù)q,主要影響角正切值tanβ,開(kāi)采影響傳播角θ,拐點(diǎn)偏移距s1、s2、s3、s4。設(shè)M=[q,tanβ,s,θ],W為M的搜索空間,若地面點(diǎn)的測(cè)量沉降量為wi,PIM計(jì)算沉降量為wi′,以二者的誤差平方和最小為準(zhǔn)則,將PIM參數(shù)的反演計(jì)算過(guò)程表示為式(6)的約束優(yōu)化問(wèn)題,即在一個(gè)確定的空間W中找到一個(gè)向量M0,使得目標(biāo)函數(shù)的值達(dá)到最小[15],即

(6)

式中,fp為目標(biāo)函數(shù),m為引入的觀測(cè)站點(diǎn)數(shù)量。

在算法設(shè)計(jì)過(guò)程中,需反復(fù)變換參數(shù)進(jìn)行調(diào)試分析,以確定改進(jìn)的PSO算法參數(shù)。設(shè)置初始的種群規(guī)模為300,最大迭代次數(shù)為2 000,學(xué)習(xí)因子c2=c3=0.5,慣性系數(shù)ω=0.9,粒子飛行速度v=1.5,尋優(yōu)的空間維度為4。利用InSAR計(jì)算出沉降區(qū)邊緣高相干點(diǎn)的形變值,并用位于走向和傾向線上的水準(zhǔn)數(shù)據(jù)對(duì)PIM模型參數(shù)進(jìn)行反演。參數(shù)反演結(jié)果見(jiàn)表2。

表2 改進(jìn)PSO算法反演結(jié)果準(zhǔn)確性驗(yàn)證

地下采區(qū)經(jīng)驗(yàn)設(shè)計(jì)值根據(jù)同類型的大紅山銅礦開(kāi)采經(jīng)驗(yàn)、礦山實(shí)際開(kāi)采方式和礦山地質(zhì)巖性等資料獲得,相比于經(jīng)驗(yàn)設(shè)計(jì)值法,利用改進(jìn)型PSO算法反演得到的開(kāi)采影響傳播角相對(duì)誤差最小,主要影響角正切值相對(duì)誤差最大,反演得到的參數(shù)值與相關(guān)經(jīng)驗(yàn)設(shè)計(jì)值存在一定誤差,還需進(jìn)一步分析。

3.3 融合結(jié)果分析

由于礦區(qū)在開(kāi)采后地表形變規(guī)律易受地質(zhì)、地形等因素影響,為表示該工作面開(kāi)采后的地表形變情況,將SBAS-InSAR與PIM進(jìn)行融合,融合后的下沉曲面見(jiàn)圖7。由圖7可知,擬合效果較好,能夠反映出地表的沉降情況,沉降規(guī)律符合PIM模型[16]。將得到的沉降結(jié)果轉(zhuǎn)至WGS-84坐標(biāo)系,采用克里金插值法獲得具有連續(xù)性的結(jié)合PIM的礦區(qū)沉降結(jié)果(見(jiàn)圖8)。

從圖7、圖8中可以看出,在研究時(shí)段內(nèi),探測(cè)到該工作面上方存在漏斗狀沉降區(qū),最大沉降值為175 mm,周邊沉降值遞減,最終在沉降區(qū)邊緣逐漸收斂,其收斂速度在走向線上沿西北方向較慢、東南方向較快,在傾向線上均以較快的速度收斂。其原因是礦區(qū)沉降受地形、斷層及其他外部因素的影響,而在PIM中僅考慮開(kāi)采引起的沉降,故監(jiān)測(cè)到的地面沉降集中在開(kāi)采工作面的上方[3,11,17]。

圖8 基于改進(jìn)型PSO算法的預(yù)計(jì)下沉曲面

圖9 插值后連續(xù)沉降區(qū)

3.4 融合法監(jiān)測(cè)結(jié)果精度分析

為了對(duì)融合法得到的監(jiān)測(cè)精度進(jìn)行驗(yàn)證,將其與該工作面上走向及傾向線上8個(gè)水準(zhǔn)點(diǎn)的監(jiān)測(cè)數(shù)據(jù)進(jìn)行對(duì)比分析,結(jié)果見(jiàn)圖10、圖11。由圖10、圖11可知,采用融合方法獲取的形變量與水準(zhǔn)監(jiān)測(cè)得到的形變量較接近,而采用經(jīng)驗(yàn)設(shè)計(jì)值法得到的形變量與水準(zhǔn)監(jiān)測(cè)得到的形變量存在明顯差異。為進(jìn)一步判斷模型的準(zhǔn)確性,利用圖8中未參與模型反演的黃色點(diǎn)監(jiān)測(cè)值進(jìn)行精度分析,結(jié)果見(jiàn)圖12。圖12中,K1-K3表示邊緣部分GPS監(jiān)測(cè)點(diǎn),H6和G5表示走向和傾向線附近水準(zhǔn)點(diǎn),通過(guò)對(duì)比發(fā)現(xiàn),構(gòu)建出的沉降模型符合現(xiàn)場(chǎng)監(jiān)測(cè)值的變化趨勢(shì)。

圖10 走向線沉降量對(duì)比圖

圖11 傾向線沉降量對(duì)比圖

圖12 現(xiàn)場(chǎng)監(jiān)測(cè)點(diǎn)沉降量對(duì)比圖

在相同的條件下,對(duì)不同方法在走向和傾向線上的形變結(jié)果誤差進(jìn)行統(tǒng)計(jì),結(jié)果見(jiàn)表3。從表3中可以看出,融合方法的平均絕對(duì)誤差為43 mm,最大絕對(duì)誤差為70 mm,最大沉降量絕對(duì)誤差為52 mm,均方根誤差為46 mm,聯(lián)合方法的擬合效果明顯優(yōu)于經(jīng)驗(yàn)設(shè)計(jì)值法,說(shuō)明融合時(shí)序InSAR與改進(jìn)型PSO算法可以用于開(kāi)采沉陷預(yù)計(jì)參數(shù)的求取,并能獲得精度較高的沉降區(qū)模型。誤差產(chǎn)生的主要原因是PIM為一個(gè)理想化的模型[18],未考慮斷層、排土場(chǎng)和降雨等因素對(duì)地表沉降的影響。

表3 實(shí)測(cè)值與融合值、經(jīng)驗(yàn)設(shè)計(jì)值的總體誤差對(duì)比

4 結(jié)論

本文針對(duì)時(shí)序InSAR在礦山采空區(qū)地表沉降監(jiān)測(cè)中的不足,提出了一種融合時(shí)序InSAR與改進(jìn)型PIM的礦區(qū)地表沉降監(jiān)測(cè)方法。以大紅山礦區(qū)某采空區(qū)為研究對(duì)象,利用該方法獲取在研究時(shí)段內(nèi)的沉降區(qū)模型,并與時(shí)序InSAR在同一時(shí)段內(nèi)的水準(zhǔn)數(shù)據(jù)以及由經(jīng)驗(yàn)設(shè)計(jì)值法得到的沉降結(jié)果進(jìn)行對(duì)比,得到以下主要結(jié)論:

a.采用時(shí)序InSAR技術(shù)與改進(jìn)型PIM融合的方法,彌補(bǔ)了InSAR在沉降中心附近由于失相干而導(dǎo)致的形變結(jié)果缺失的問(wèn)題,且僅結(jié)合少量水準(zhǔn)數(shù)據(jù)便獲取了礦山地下開(kāi)采引起的地表沉降區(qū)完整模型。

b.融合改進(jìn)型PSO算法與時(shí)序InSAR技術(shù)反演出的PIM參數(shù)效果較好,在走向線和傾向線上誤差結(jié)果均優(yōu)于由經(jīng)驗(yàn)設(shè)計(jì)值法得到的沉降結(jié)果,且各項(xiàng)誤差均不超過(guò)100 mm,說(shuō)明該方法滿足實(shí)際工程需求,為時(shí)序InSAR技術(shù)在礦山開(kāi)采沉陷監(jiān)測(cè)中的應(yīng)用提供了新思路。