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

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

0 引言

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

利用InSAR對采空區(qū)地表進(jìn)行監(jiān)測時(shí),由于礦山開采引起的地表沉降速度快、形變量大,且SAR衛(wèi)星存在一定的重復(fù)周期,導(dǎo)致沉降區(qū)中心附近出現(xiàn)低相干、失相干現(xiàn)象,形變結(jié)果缺失,無法獲取準(zhǔn)確的沉降區(qū)模型[4]。為解決上述問題,當(dāng)前采用的主要方法是概率積分法(PIM)[5],PIM是由劉寶琛基于隨機(jī)介質(zhì)理論而提出的,其主要應(yīng)用于煤礦規(guī)則采空區(qū)。PIM利用高精度監(jiān)測點(diǎn)對相關(guān)參數(shù)進(jìn)行反演,以解決采空區(qū)地表沉降中心相干性較低的問題。為得到準(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ǎo)致局部最優(yōu)解而使收斂速度過快,其在采空區(qū)沉降監(jiān)測中的精度較差[8]。

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

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

1.1 研究區(qū)概況

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

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

1.2 前期研究存在的問題

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

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

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

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

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

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

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

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

2 研究方法

2.1 PIM模型

PIM是一種非連續(xù)的隨機(jī)介質(zhì)模型,由于巖層在移動過程中會破壞介質(zhì)的連續(xù)性,介質(zhì)單元發(fā)生變化,導(dǎo)致沉積環(huán)境也發(fā)生變化,從而產(chǎn)生不同的新巖層,在地質(zhì)構(gòu)造作用下,新巖層中形成節(jié)理、斷層等弱面,弱面即為非連續(xù)的隨機(jī)介質(zhì)[11-12]。礦層開采后巖層與地表移動是一種隨機(jī)的過程,形成的地表塌陷區(qū)見圖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開采單元,H為工作面采深,m為礦層厚度,q為下沉系數(shù),α為礦層傾角,θ為開采影響傳播角,tanβ為主要影響角的正切值,D1、D3分別為工作面的走向、傾向長,s1、s2分別為走向上的左右邊界拐點(diǎn)偏移距,s3、s4分別為上山下山方向的拐點(diǎn)偏移距[11]。

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

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

在對種群進(jìn)行初始化處理的基礎(chǔ)上,選擇合適的參數(shù),迭代計(jì)算每個(gè)粒子的適應(yīng)值,將個(gè)體極值和種群極值代入式(4)、式(5),不斷更新每個(gè)粒子的位置和移動方向,最終得到所需的尋優(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.對覆蓋礦區(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算法對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é)果。具體的工作流程見圖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)位,并對該工作面上方的走向和傾向線的7個(gè)水準(zhǔn)點(diǎn)進(jìn)行線性插值,使InSAR沉降結(jié)果與水準(zhǔn)數(shù)據(jù)在時(shí)間上統(tǒng)一,共選取37個(gè)特征點(diǎn),特征點(diǎn)位置見圖7(白色圓點(diǎn)表示InSAR高相干性點(diǎn),黑色三角點(diǎn)表示走向和傾向線水準(zhǔn)監(jiān)測點(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β,開采影響傳播角θ,拐點(diǎn)偏移距s1、s2、s3、s4。設(shè)M=[q,tanβ,s,θ],W為M的搜索空間,若地面點(diǎn)的測量沉降量為wi,PIM計(jì)算沉降量為wi′,以二者的誤差平方和最小為準(zhǔn)則,將PIM參數(shù)的反演計(jì)算過程表示為式(6)的約束優(yōu)化問題,即在一個(gè)確定的空間W中找到一個(gè)向量M0,使得目標(biāo)函數(shù)的值達(dá)到最小[15],即

(6)

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

在算法設(shè)計(jì)過程中,需反復(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ù)對PIM模型參數(shù)進(jìn)行反演。參數(shù)反演結(jié)果見表2。

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

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

3.3 融合結(jié)果分析

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

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

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

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

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

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

圖10 走向線沉降量對比圖

圖11 傾向線沉降量對比圖

圖12 現(xiàn)場監(jiān)測點(diǎn)沉降量對比圖

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

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

4 結(jié)論

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

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

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