基于SBAS-InSAR技術(shù)和Logistic模型的礦區(qū)沉降動(dòng)態(tài)預(yù)測方法

2022-06-22 07:03:18徐子興季民張過陳振煒

自然資源遙感 2022年2期

關(guān)鍵詞：模型

徐子興，季民，張過，陳振煒

(1.山東科技大學(xué)測繪與空間信息學(xué)院，青島 266590； 2.武漢大學(xué)測繪遙感信息工程國家重點(diǎn)實(shí)驗(yàn)室，武漢 430079)

0 引言

地下采礦活動(dòng)會破壞開采工作面上覆巖層的原始應(yīng)力平衡，從而導(dǎo)致礦區(qū)地表發(fā)生水平移動(dòng)和下沉，不僅會對地面建筑和生產(chǎn)設(shè)施造成破壞，還可能誘發(fā)地質(zhì)災(zāi)害，給礦區(qū)居民的生產(chǎn)和生活留下安全隱患[1]。地下開采導(dǎo)致的地表沉降不僅與煤層的水文地質(zhì)條件有關(guān)，而且也與開采方法、開采速度等密切相關(guān)。利用礦山開采沉降規(guī)律對地下開采活動(dòng)引起的后續(xù)地表沉降進(jìn)行有效預(yù)測，是評估礦山開采風(fēng)險(xiǎn)并進(jìn)一步采取防護(hù)措施的關(guān)鍵。

以往基于水準(zhǔn)或GPS測量的沉降僅可對少量離散點(diǎn)進(jìn)行沉降預(yù)測，而合成孔徑雷達(dá)干涉測量(synthetic aperture Radar interferometry，InSAR)技術(shù)的發(fā)展和應(yīng)用為大范圍的沉降預(yù)測提供了數(shù)據(jù)基礎(chǔ)。InSAR技術(shù)憑借其大范圍監(jiān)測、高空間分辨率、穿云透霧等優(yōu)點(diǎn)被廣泛應(yīng)用于地表形變的監(jiān)測中[2]。礦區(qū)地表通常有植被覆蓋，并且容易發(fā)生大梯度形變，這使得礦區(qū)容易發(fā)生失相干現(xiàn)象[3]。而小基線集干涉測量技術(shù)(small baseline subset-InSAR，SBAS-InSAR)作為一種常用的InSAR形變監(jiān)測技術(shù)，可以有效減弱時(shí)空失相干和大氣效應(yīng)的影響，在礦區(qū)沉降監(jiān)測中得到了廣泛的研究與應(yīng)用[4-7]。

礦區(qū)開采沉降過程大致可分為3個(gè)階段：初始沉降期、主要沉降期和殘余沉降期，沉降增長呈“S”型[8]。描述該過程的時(shí)間函數(shù)模型主要有Knothe模型[9-11]、Richards模型[12-13]、正態(tài)分布時(shí)間函數(shù)[14-15]及Logistic模型[16-17]等。其中，Logistic模型作為典型的“S”型增長曲線，與礦山開采沉降的3個(gè)階段比較符合。張文志等[16]用Logistic模型曲線擬合開采沉陷-時(shí)間關(guān)系曲線，并進(jìn)行了單點(diǎn)開采沉陷預(yù)測，結(jié)果表明該模型應(yīng)用于礦區(qū)沉降擬合和預(yù)測具有較高的精度。楊澤發(fā)等[17]通過將合成孔徑雷達(dá)差分干涉測量(differential interferometry synthetic aperture Radar，D-InSAR)監(jiān)測到的礦區(qū)沉降累加，結(jié)合Logistic模型對礦區(qū)全盆地沉降規(guī)律進(jìn)行了分析，并對全盆地的后續(xù)沉降進(jìn)行了預(yù)測。但僅當(dāng)整個(gè)沉降盆地的沉降過程接近結(jié)束時(shí)，全盆地的所有點(diǎn)才滿足可對后續(xù)沉降進(jìn)行預(yù)測的條件。在煤炭開采過程中，整個(gè)下沉盆地會有部分點(diǎn)滿足預(yù)測條件、部分點(diǎn)不滿足預(yù)測條件，上述不考慮可用條件的預(yù)測方法并不適用。

本文對使用Logistic時(shí)間函數(shù)模型進(jìn)行礦區(qū)沉降預(yù)測中的可用條件進(jìn)行了分析和模擬實(shí)驗(yàn)，并提出了一種基于小基線集InSAR(SBAS-InSAR)技術(shù)和Logistic時(shí)間函數(shù)模型的礦區(qū)沉降動(dòng)態(tài)預(yù)測方法。首先，采用SBAS-InSAR技術(shù)對礦山開采引起的沉降進(jìn)行監(jiān)測，獲取礦區(qū)沉降時(shí)間序列; 然后，以時(shí)序沉降數(shù)據(jù)作為擬合數(shù)據(jù)，采用信賴域算法逐像元計(jì)算Logistic模型參數(shù)，并根據(jù)Logistic模型的可用條件判斷出可對后續(xù)沉降進(jìn)行預(yù)測的像元范圍; 最后，根據(jù)Logistic模型計(jì)算對可預(yù)測范圍內(nèi)的所有像元的后續(xù)沉降進(jìn)行預(yù)測。

1 方法與原理

1.1 Logistics模型

選擇Logistic模型以確定礦區(qū)沉降動(dòng)態(tài)預(yù)測中的可預(yù)測條件和預(yù)測礦區(qū)沉降，這是因?yàn)槠鋵τ诘V區(qū)沉降具有良好的擬合精度，且具有最大加速度點(diǎn)、最大速度點(diǎn)和最小加速度點(diǎn)3個(gè)特征點(diǎn)，利于沉降階段的劃分。Logistic模型是經(jīng)典的“S”型增長曲線，其公式為：

，

(1)

式中：W(t)為某一點(diǎn)t時(shí)刻的下沉值；W0為該點(diǎn)最大下沉值；a和b為形狀參數(shù)；待求參數(shù)為W0，a和b。對式(1)分別求一階導(dǎo)數(shù)和二階導(dǎo)數(shù)，可得到Logistic模型下沉速度和下沉加速度的表達(dá)式，分別為：

，

(2)

(3)

Logistic模型具有3個(gè)特征點(diǎn)，分別為最大加速度點(diǎn)、最大速度點(diǎn)和最小加速度點(diǎn)。根據(jù)這3個(gè)特征點(diǎn)設(shè)置了擬合數(shù)據(jù)對沉降過程的不同覆蓋程度，并探究其對Logistic模型擬合精度的影響。假定一個(gè)W0，a和b參數(shù)分別為-1 m，120和0.04的Logistic模型，時(shí)間范圍設(shè)置為0～240 d，每12 d獲取一個(gè)沉降值得到該Logistic模型下的沉降時(shí)間序列。為了使模擬的數(shù)據(jù)符合實(shí)際情況，在沉降時(shí)間序列中加入了0.02 m的正態(tài)分布隨機(jī)誤差。根據(jù)擬合數(shù)據(jù)對沉降過程的不同覆蓋程度設(shè)置了4組模擬實(shí)驗(yàn)，采用信賴域算法進(jìn)行擬合得到Logistic曲線，如圖1所示。圖上3個(gè)紅色三角形自左到右分別代表最大加速度點(diǎn)、最大速度點(diǎn)和最小加速度點(diǎn)。

(a) 覆蓋范圍未超過最大加速度點(diǎn)(b) 覆蓋范圍超過最大加速度點(diǎn)

從圖1中可以看出，隨著擬合數(shù)據(jù)覆蓋范圍的增加，擬合曲線越來越逼近正確曲線。直到擬合數(shù)據(jù)覆蓋范圍超過最右邊特征點(diǎn)(最小加速度點(diǎn))時(shí)，擬合曲線才非常接近正確曲線。因此，將擬合數(shù)據(jù)覆蓋范圍是否超過最小加速度點(diǎn)作為Logistic模型的可用條件，并根據(jù)該可用條件逐像元進(jìn)行判斷得到可對后續(xù)沉降進(jìn)行預(yù)測的范圍。

1.2 SBAS-InSAR礦區(qū)沉降監(jiān)測

采用SBAS-InSAR對礦山開采引起的沉降進(jìn)行監(jiān)測，SBAS-InSAR可以有效減弱常規(guī)InSAR方法中易受大氣延遲誤差干擾、時(shí)空失相干等問題。假設(shè)有N+1幅覆蓋研究區(qū)域的SAR影像時(shí)，選取其中一幅作為超級主影像，將其余N幅SAR影像都配準(zhǔn)到超級主影像上，根據(jù)小基線原則和礦區(qū)形變速度快、形變量級大的特點(diǎn)設(shè)定時(shí)空基線閾值進(jìn)行基線組合，得到M個(gè)干涉對。M滿足:

(N+1)/2

，

(4)

其中任一干涉對主影像和從影像獲取時(shí)間分別設(shè)為tA和tB(tB>tA)，對干涉對進(jìn)行干涉、去平地效應(yīng)和去地形相位后，得到的第j幅干涉圖上高相干點(diǎn)(x,y)上的相位值可以表示為[18]:

，

(5)

式中：λ為波長；d(tB,x,y)和d(tA,x,y)為相對于d(t0,x,y)=0時(shí)的雷達(dá)視線方向的累積形變量；Δφjtopo(x,y)，Δφjatm(x,y)，Δφjorb(x,y)和Δφjnoise(x,y)分別表示殘余地形、軌道信息不精確、大氣延遲和噪聲導(dǎo)致的相位誤差，去除這些誤差相位后，式(5)可簡化為：

(6)

d(tB,x,y)-d(tA,x,y)=vi(tB-tA)

，

(7)

式中：vi為tA至tB時(shí)間段內(nèi)雷達(dá)視線方向的平均形變速率。將式(6)和式(7)轉(zhuǎn)化為矩陣形式為：

BV=δφ

，

(8)

式中：B為參數(shù)矩陣；V為平均形變速率矩陣；δφ為相位矩陣。采用最小二乘法和奇異值分解法對式(8)求解，得到各時(shí)間段內(nèi)平均形變速率的最小范數(shù)解，并對其進(jìn)行積分得到視線(line of sight，LOS)方向的時(shí)間序列形變量[19]。

利用式(9)將LOS方向形變轉(zhuǎn)換為沉降值，公式為：

W=dLOS/cosθ

，

(9)

式中：W為沉降值；dLOS為LOS方向形變；θ為入射角。

1.3 信賴域算法求解Logistic模型參數(shù)

采用信賴域算法逐像元將SBAS-InSAR監(jiān)測的時(shí)序沉降進(jìn)行Logistic模型擬合，計(jì)算其最優(yōu)Logistic模型參數(shù)，即求解式(10)：

(10)

設(shè)xk為第k次迭代點(diǎn)，fk=f(xk),gk=g(xk) =f(xk)，Bk為2f(xk)的第k次近似，則第k次迭代步信賴域二次模型子問題如式(11)所示[20]：

，

(11)

式中： Δk>0為信賴域半徑。設(shè)子問題的解為dk，通過式(12)反映二次函數(shù)模型與目標(biāo)函數(shù)之間的近似程度，公式為：

，

(12)

若rk< 0，則xk+dk不能作為下一個(gè)迭代點(diǎn)，此時(shí)需要縮小信賴域半徑；若rk接近于1，可以認(rèn)為模型函數(shù)與目標(biāo)函數(shù)在信賴域范圍內(nèi)具有良好的一致性，可以擴(kuò)大信賴域半徑；反之，若rk離1 較遠(yuǎn)，則需要縮小信賴域半徑。不斷迭代直至得到收斂解，流程圖如圖2所示。

1.4 礦區(qū)沉降動(dòng)態(tài)預(yù)測方法

本預(yù)測方法首先采用SBAS-InSAR對開采沉降進(jìn)行監(jiān)測，獲取沉降時(shí)間序列；然后以沉降時(shí)間序列作為擬合數(shù)據(jù)，采用信賴域算法逐像元進(jìn)行Logistic模型求參，并根據(jù)可預(yù)測條件判斷可預(yù)測范圍；最后利用Logistic模型對可預(yù)測范圍內(nèi)的像元進(jìn)行沉降預(yù)測。流程圖如圖3所示。

圖3 礦區(qū)沉降動(dòng)態(tài)預(yù)測方法流程圖

2 礦區(qū)沉降預(yù)測實(shí)驗(yàn)與結(jié)果

2.1 研究區(qū)與數(shù)據(jù)源

所選研究礦區(qū)為內(nèi)蒙古某礦區(qū)，位于鄂爾多斯市東南部，地處毛烏素沙漠的東北邊緣，與陜西省榆林市毗鄰，如圖4所示。研究礦區(qū)表面植被覆蓋率較低，受時(shí)間去相干影響較小。所處位置地勢較為平坦，平均海拔為1 310 m。利用40景IW模式、VV極化的Sentinel-1A影像進(jìn)行SBAS-InSAR處理，時(shí)間跨度為2017年9月20日—2019年2月18日，影像參數(shù)如表1所示。其中2017年9月20日—2018年11月02日的時(shí)序沉降用于作為擬合數(shù)據(jù)計(jì)算Logistic模型參數(shù)，并分別對2018年12月08日和2019年02月18日的累積沉降進(jìn)行預(yù)測，采用對應(yīng)日期的InSAR監(jiān)測累積沉降對預(yù)測結(jié)果進(jìn)行了驗(yàn)證。相鄰影像時(shí)間間隔大部分是12 d，少量為24 d，小的時(shí)間間隔可以最小化大部分去相干源的影響。采用30 m分辨率的航天飛機(jī)雷達(dá)地形測繪使命(Shuttle Radar Topography Mission，SRTM)數(shù)據(jù)用于模擬和去除地形相位。

圖4 研究礦區(qū)位置

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

2.2 SBAS-InSAR監(jiān)測

根據(jù)小基線集的原則和礦區(qū)大梯度形變的特點(diǎn)，設(shè)置空間基線和時(shí)間基線的閾值分別為150 m和36 d。經(jīng)過配準(zhǔn)、差分干涉、濾波、相位解纏、基線精化、估算線性形變速率、去除大氣延遲誤差、估算最終形變速率、相位轉(zhuǎn)形變和地理編碼，并將形變轉(zhuǎn)到垂直方向后，得到了該礦區(qū)2017年9月20日—2019年2月18日的39幅時(shí)序沉降圖，部分如圖5所示。

圖5-1 時(shí)序沉降圖

(e) 2017/09/20—2018/08/22(f) 2017/09/20—2018/10/21(g) 2017/09/20—2018/12/20(h) 2017/09/20—2019/02/18

圖5-2 時(shí)序沉降圖

圖5(h)中的A—D這4個(gè)點(diǎn)于下節(jié)中展示Logistic模型擬合情況。

2.3 Logistic模型擬合

以SBAS-InSAR處理得到的2017年9月20日—2018年11月02日的時(shí)序沉降數(shù)據(jù)作為擬合數(shù)據(jù)，對最終沉降大于0.05 m的像元逐像元采用信賴域算法計(jì)算Logistic模型參數(shù)。圖6展示圖5(h)中的A，B，C，D這4個(gè)點(diǎn)的Logistic模型擬合情況。

(a) A點(diǎn)(b) B點(diǎn)

2.4 沉降預(yù)測結(jié)果

根據(jù)擬合得到的Logistic模型的參數(shù)，以參與擬合數(shù)據(jù)的覆蓋范圍是否超過最小加速度點(diǎn)作為判斷依據(jù)逐像元進(jìn)行判斷，得到可預(yù)測范圍。分別對2017年9月20日—2018年12月08日和2019年02月18日的累積沉降進(jìn)行預(yù)測，如圖7所示。

(a) 2017/09/20—2018/12/08(b) 2017/09/20—2019/02/18

圖7 預(yù)測沉降圖

2.5 預(yù)測精度驗(yàn)證

采用SBAS-InSAR獲得的2017年9月20日—2018年12月08日和2019年02月18日的累積沉降對預(yù)測結(jié)果進(jìn)行驗(yàn)證。將InSAR監(jiān)測值和預(yù)測值作一致性驗(yàn)證圖，如圖8所示，可以看出大部分點(diǎn)分布在y=x附近，兩者一致性較高。可預(yù)測范圍內(nèi)全部像元的預(yù)測沉降值與InSAR監(jiān)測沉降值的RMSE分別為0.010 1 m和0.023 6 m。統(tǒng)計(jì)預(yù)測誤差作柱狀圖如圖9所示，預(yù)測誤差小于0.03 m的比例分別達(dá)到98.9%和89.3%，表明該動(dòng)態(tài)預(yù)測模型預(yù)測精度較高。

(a) 2017/09/20—2018/12/08 (b) 2017/09/20—2019/02/18

2.6 預(yù)測精度影響因素分析

預(yù)測誤差的空間分布如圖10所示，從圖中可以看出：垂直于工作面開采方向，預(yù)測誤差較大的像元基本處于沉降盆地的中間位置，因?yàn)槌两蹬璧刂虚g位置下沉量大、下沉速度快，容易導(dǎo)致大的預(yù)測誤差；沿著工作面開采方向，預(yù)測誤差較大的像元大部分處于開采方向末端，分析原因是開采方向末端像元擬合數(shù)據(jù)對最小加速度點(diǎn)后下沉過程的覆蓋程度較低。將最小加速度點(diǎn)出現(xiàn)時(shí)間作為x軸，預(yù)測誤差作為y軸，繪制散點(diǎn)圖，如圖11所示。從圖11可以看出，最小加速度點(diǎn)出現(xiàn)時(shí)間越晚，出現(xiàn)較大預(yù)測誤差的概率也隨之增加。因?yàn)閿M合數(shù)據(jù)的時(shí)間范圍是一個(gè)定值，最小加速度點(diǎn)出現(xiàn)時(shí)間越晚，最小加速度點(diǎn)之后的擬合數(shù)據(jù)就越少，對最小加速度點(diǎn)后下沉過程的覆蓋程度就會越低，會對后續(xù)沉降的預(yù)測精度產(chǎn)生影響。

(a) 2017/09/20—2018/12/08(b) 2017/09/20—2019/02/18

3 結(jié)論

本文對Logistic模型在礦區(qū)沉降動(dòng)態(tài)預(yù)測中的可用條件進(jìn)行分析與模擬實(shí)驗(yàn)，并基于SBAS-InSAR技術(shù)和Logistic模型提出了一種新的礦區(qū)沉降動(dòng)態(tài)預(yù)測方法。以內(nèi)蒙古某礦區(qū)為實(shí)驗(yàn)區(qū)，利用2017年9月20日—2019年2月18日的40景Sentinel-1A數(shù)據(jù)進(jìn)行了礦區(qū)開采沉降預(yù)測實(shí)驗(yàn)，結(jié)果顯示36 d后和108 d后預(yù)測結(jié)果的RMSE分別為0.010 1 m和0.023 6 m，預(yù)測誤差小于0.03 m的比例分別達(dá)到98.9%和89.3%，表明該預(yù)測方法精度較高。該方法具有一定的實(shí)用性，對評估礦山開采風(fēng)險(xiǎn)、指導(dǎo)礦區(qū)開采規(guī)劃具有一定的理論借鑒價(jià)值和實(shí)際應(yīng)用意義。