劉翠芝, 王興杰, 賀黎明, 唐永亮

(1. 東北大學 資源與土木工程學院, 遼寧 沈陽 110819; 2. 遼寧省地質(zhì)環(huán)境監(jiān)測總站, 遼寧 沈陽 110032)

煤礦開采帶來的采空沉陷是我國諸多資源開采型城市面臨的嚴重的地質(zhì)災害之一[1],由于煤采區(qū)開采的持續(xù)性和區(qū)域多煤層地質(zhì)結(jié)構(gòu)的復雜性,地表沉降主要表現(xiàn)為非線性和大梯度的特征,給當?shù)氐慕?jīng)濟發(fā)展和人民生活帶來了極大的危害,嚴重制約資源型城市的可持續(xù)性發(fā)展[2].煤采區(qū)地表沉降是地下煤炭資源被持續(xù)采出后,開采區(qū)域周圍巖體的原始應力平衡受到破壞,造成上覆巖層和地表失去支撐而發(fā)生強烈的非連續(xù)和非線性形變,且隨著時間推移應力重新分布達到新平衡的復雜過程,是開采煤層幾何形狀、位置等幾何狀態(tài)參數(shù)在地表的綜合直觀體現(xiàn).因此開展對煤田采空區(qū)的形變監(jiān)測和煤層參數(shù)反演的研究對區(qū)域地質(zhì)災害防范與煤礦安全生產(chǎn)有重要的意義[3-4].

傳統(tǒng)的全球?qū)Ш叫l(wèi)星系統(tǒng)(global navigation satellite system,GNSS)和水準等點式測量方法無法實現(xiàn)對沉陷區(qū)大范圍的有效監(jiān)測,而雷達差分干涉測量技術(shù)(interferometric synthetic aperture radar, InSAR)是一種典型的面式對地觀測手段,具有高精度、高分辨率、空間連續(xù)覆蓋的典型特征,特別適合大范圍、大梯度的地表沉降監(jiān)測[5-7].基于InSAR監(jiān)測獲得的覆蓋整個煤采區(qū)的地表形變場,結(jié)合地球物理非線性反演算法,可快速準確反演得到煤層幾何形狀、位置以及產(chǎn)狀等參數(shù)信息,這對于分析近距離多煤層開采引起地表形變的機理和演變規(guī)律至關(guān)重要.

近年,國內(nèi)外學者對煤層參數(shù)的反演研究逐漸增多,但大多集中于兩個方向.一種是基于FLAC3D有限元或UDEC離散元等數(shù)值模型來對近距離多煤層共同采動引起的地表位移進行分析[8-9],這類方法對近距離多煤層開采深入階段的實際狀況進行了較為周密的考慮,具有較高的實際應用指導價值,但由于其模擬過程中需要已知詳盡的多煤層位置、產(chǎn)狀等先驗參數(shù)信息,并且缺乏實際監(jiān)測的地表形變約束,因此嚴重制約了在很多未知先驗信息的煤層參數(shù)反演過程中的應用;另一種是依托InSAR監(jiān)測獲得的地面沉降量信息,基于Okada模型選取部分單個煤層巷道反演其幾何狀態(tài)參數(shù)[10],這類方法可以對單個煤層的幾何位置等參數(shù)進行合理的反演解釋,但單煤層巷道不能反映多煤層開采的實際狀況,因此缺乏對近距離多煤層實際開采狀態(tài)的合理表達.

為此,本文提出了一種基于多源Okada模型的受地表形變約束的近距離多煤層參數(shù)反演方法,有效避免FLAC3D和UDEC等數(shù)值模型不受實際地表形變約束的缺陷.該方法依托小基線集技術(shù)(small baseline subset,SBAS)獲取的研究區(qū)地面沉降信息作為觀測數(shù)據(jù),實現(xiàn)多源Okada模型在近距離多煤層參數(shù)反演中的應用.本文的目的是通過研究近距離多煤層開采引發(fā)地表變形的規(guī)律,以實測地表形變作為約束,反演得到近距離多煤層的幾何形狀、位置以及產(chǎn)狀等參數(shù)信息,對于我國普遍分布的煤采區(qū),尤其是大量已關(guān)閉或廢棄煤礦的地面沉陷范圍圈定以及多煤層采空位置、產(chǎn)狀等信息獲取具有重要借鑒意義,進而可為煤采區(qū),特別是已關(guān)閉或廢棄工礦區(qū)地面沉降監(jiān)測和綜合治理提供依據(jù).

1 近距離多煤層參數(shù)反演方法

近距離多煤層共同開采是目前全球范圍內(nèi)各礦場在實際開采作業(yè)過程中所面臨的一種復雜情況.研究其開采規(guī)律,并根據(jù)實際監(jiān)測的地表形變快速準確反演地下復雜的多煤層幾何形狀、位置、產(chǎn)狀等參數(shù)對于煤礦采空區(qū)的綜合治理極其重要.本文基于Okada模型[11],結(jié)合近距離多煤層分布的實際特征,依托疊加原理,提出了一種基于多源模型的近距離多煤層參數(shù)反演方法,其示意圖如圖1所示.

圖1 近距離多煤層的多源模型示意圖Fig.1 The multi-source model for the close multi-coalseam(a)—煤層分布等效圖； (b)—位錯模型示意圖； (c)—多源位錯模型分布圖.

圖1a是近距離多煤層開采過程中多個開采層和中間層的等效表示,圖1b是Okada有限矩形源模型,共有10個參數(shù),分別為開采面幾何中心P在地表的投影坐標(X,Y)、開采面沿著走向的長度L、沿傾向的寬度W、煤層開采中心深度depth、走向角strike(開采掘進方向與正北方向的夾角)、傾向角dip(開采面與水平方向的夾角)、rake(斷層破裂時相對于斷層走向移動的角度)、破裂方向的位錯分量slip、張性分量opening;圖1c是依托多煤層的分布情況,將Okada有限矩形源模型應用于多煤層的綜合圖.

基于多源模型的近距離多煤層參數(shù)反演方法首先使用基于梯度的自適應四叉樹采樣對InSAR監(jiān)測數(shù)據(jù)進行降采樣,通過對較高位移梯度的區(qū)域進行更精細的劃分來減少InSAR圖像內(nèi)部和邊緣的數(shù)據(jù)間隙的影響,在保證有效采樣的基礎上降低數(shù)據(jù)點位密度,提升反演計算效率.同時通過計算較高位移梯度和數(shù)據(jù)邊緣信息,再結(jié)合InSAR形變圖的空間分布特征和研究區(qū)多煤層的分布特征,預估多源模型的數(shù)量和模型參數(shù)的初始值等.然后利用已有的先驗信息和似然函數(shù)構(gòu)建形變觀測數(shù)據(jù)d和通過疊加理論計算得到的多源模型形變模擬數(shù)據(jù)G(t)之間的數(shù)學關(guān)系,接下來采用非線性貝葉斯反演算法[12-13]進行多煤層參數(shù)反演,在反演過程中,使用抽樣計算模型參數(shù)的后驗概率分布來判定多煤層參數(shù)擬合結(jié)果.多源模型形變模擬數(shù)據(jù)G(t)和后驗概率分布的計算公式如下:

G(t)=b1G(t)1+b2G(t)2+…+bsG(t)s,

(1)

(2)

其中:G(t)s和bs分別表示多源模型中不同源模擬得到的形變和所占權(quán)重;s為多源模型的數(shù)量;d為M維觀測數(shù)據(jù);t為N維模型參數(shù);p(d)為歸一化常數(shù);p(t)為多源模型參數(shù)的先驗信息;p(t/d)為在已知先驗信息下當前參數(shù)的后驗概率;p(d/t)為多源模型的似然函數(shù)(likelihood function).

本文依據(jù)模型參數(shù)擬合值的分布特征,選擇多維高斯表達式L(t)來表示似然函數(shù),如式(3)所示,其中Cd表示數(shù)據(jù)協(xié)方差矩陣.

L(t)=

(3)

在抽樣計算多源模型后驗概率時,采樣方法的選取對于后驗概率的計算結(jié)果和效率至關(guān)重要.本文采用基于Metropolis-Hastings接受準則的馬爾科夫鏈蒙特卡羅(Markov chain Monte Carlo,MCMC)算法[14]在空間上進行采樣,該方法可以大大提高采樣接受率,彌補了傳統(tǒng)蒙特卡羅方法的抽樣分布不能隨著模擬的進行而動態(tài)改變的缺陷.采樣時,首先根據(jù)先驗信息和模擬退火(simulated annealing,SA)算法估計參數(shù)并作為初始狀態(tài)參數(shù),然后在搜尋空間上建立一個平穩(wěn)分布π(i)的各態(tài)遍歷的馬爾科夫過程,如式(4)所示:

π(i)=P(i1=qi),i=1,2,…,N.

(4)

π(i)是t=1時刻處于狀態(tài)qi的概率.之后通過狀態(tài)轉(zhuǎn)移矩陣A生成大量樣本,如式(5)所示:

A=[aij]N×N.

(5)

aij=P(it+1=qi|it=qi),i,j=1,2,…,N.

(6)

其中：aij是在t時刻處于狀態(tài)qi的條件下向t+1時刻轉(zhuǎn)移到狀態(tài)qj的概率.其中i表示狀態(tài)序列,q表示可能的狀態(tài).

最后計算所有樣本多源模型反演參數(shù)的后驗概率,作為多煤層最優(yōu)參數(shù)的選擇指標.在尋找最優(yōu)參數(shù)的過程中,選擇后驗概率最大的一組多源模型參數(shù)即為最佳擬合參數(shù).

2 研究區(qū)InSAR形變監(jiān)測

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

康平煤田位于遼寧省沈陽市康平縣,屬于地下開采的隱蔽型多煤層聯(lián)合開采煤田(圖2).其中大平煤礦轄區(qū)面積28.61 km2,年生產(chǎn)能力330萬t,共劃分為7個采區(qū),其中有5個采區(qū)處于三臺子水庫下;小康煤礦轄區(qū)面積為30.08 km2,被劃分為6個采區(qū),核定年開采量270萬t.本文選用2018年1月至2019年8月間覆蓋康平煤田的43景歐空局C波段Sentinel-1B降軌數(shù)據(jù)進行形變提取.SAR衛(wèi)星的航向角為-166.19°,入射角為37.13°.DEM數(shù)據(jù)選用德國TanDEM-X獲取的空間分辨率為3 rad·s的DEM,用于影像配準、平地效應去除、軌道精煉等.

圖2 研究區(qū)位置圖Fig.2 Location of the studied area

2.2 InSAR形變場監(jiān)測

本文采用SBAS-InSAR技術(shù)對研究區(qū)域進行沉降監(jiān)測.SBAS技術(shù)由Berardino等[15]提出,為了減弱時空失相干效應的影響,小基線集技術(shù)從獲取的SAR影像中選擇時空基線均小于某一閾值的干涉像對生成差分干涉圖,同時根據(jù)基線情況將這些干涉像對劃分成若干小基線集,再基于最小二乘法解算每個小基線集的形變時間序列,然后利用奇異值分解法(singular value decomposition,SVD)將多個小基線集聯(lián)合求解,得到整個時間段上雷達視線向(line of sight,LOS)的形變時間序列[16-17].流程見圖3.

圖3 SAR數(shù)據(jù)處理流程圖Fig.3 Flow chart of the SAR data processing

本文在獲取研究區(qū)形變結(jié)果時,基于研究區(qū)近距離多煤層開采的地表實際覆蓋情況,充分考慮時空基線和多普勒效應的影響,設置空間基線閾值為臨界基線的2%;時間基線閾值為72 d,生成了216組干涉像對(圖4).

圖4 時間-空間基線圖Fig.4 The spatio-temporal baseline

基于SBAS-InSAR技術(shù),獲取了研究區(qū)在雷達視線向的累積沉降量,結(jié)果如圖5所示.在監(jiān)測范圍內(nèi),大平礦采區(qū)最大沉降速率為-99.7 mm/a,最大累積沉降量為-136.6 mm;小康礦近距離多煤層開采引起的最大沉降速率為-124.4 mm/a,最大累積沉降量為-194.8 mm;除此之外,其他區(qū)域沉降現(xiàn)象不明顯,這表明該區(qū)域整體較為穩(wěn)定,地面沉降主要是由煤炭開采所致.

為進一步分析兩個沉降區(qū)不同位置的形變過程,分別選取了Q1和Q2兩個特征點,刻畫其沉降時間序列,結(jié)果如圖6所示.藍色實線表示大平礦采區(qū)特征點Q1的沉降時間序列,紅色實線表示小康礦采區(qū)特征點Q2的沉降時間序列.可以發(fā)現(xiàn)2018年1月～2019年8月間,大平礦沉降區(qū)的Q1點累積沉降量為-111.7 mm;小康礦沉降區(qū)Q2點的累積沉降量最大,達到-178.4 mm.持續(xù)的地下煤炭開采已造成康平煤田區(qū)域產(chǎn)生了巨大的地表形變梯度和影響范圍,使當?shù)鼐用褡≌仍O施被積水淹沒或產(chǎn)生結(jié)構(gòu)性損壞而被廢棄，如圖5所示.

圖5 康平煤田累積沉降圖Fig.5 Cumulative subsidence of the Kangping coalfield

圖6 時間序列累積沉降圖Fig.6 Diagram of the time series cumulative subsidences

3 小康礦多煤層參數(shù)反演與分析

3.1 多源模型反演結(jié)果

本文以時序SBAS-InSAR方法獲取的小康礦形變監(jiān)測數(shù)據(jù)為基礎,綜合考慮形變反演過程中數(shù)據(jù)點位密度對反演計算效率的影響,選用1.44×10-4作為自適應四叉樹方差閾值進行數(shù)據(jù)降采樣,然后通過較高位移梯度和數(shù)據(jù)邊緣信息,選取2和3作為多源模型的模型數(shù),分別使用非線性貝葉斯反演算法對近距離多煤層參數(shù)進行反演.Okada模型中slip表示斷層破裂方向的位錯分量,rake表示斷層破裂時相對于斷層走向移動的角度,opening表示斷層垂直方向運動的張性分量.通過分析小康礦形變區(qū)的地表沉降特征,發(fā)現(xiàn)小康礦開采過程中多煤層傾角很緩,近乎水平開采,且地表沉降主要集中在垂直方向,因此本文簡化了Okada模型,假設斷層傾角為0°,只分析多煤層開采導致垂直方向的張性破裂引發(fā)的地表形變.在反演過程中使用獨立的測量坐標系,參考基準為[123.395°E, 42.65°N],根據(jù)小康礦區(qū)域煤系地層巖性,泊松系數(shù)確定為0.25.依托以上信息,使用基于雙源模型和三源模型的近距離多煤層參數(shù)反演方法分別反演小康礦近距離多煤層參數(shù),獲得最佳擬合結(jié)果如表1所示.

表1 小康礦多源模型最佳擬合參數(shù)Table 1 Best fitted parameters of the multi-source model of the Xiaokang mine

對比研究區(qū)近距離多煤層實際開采參數(shù)與多源模型反演最佳擬合參數(shù),評價多源模型參數(shù)反演方法的可靠性.本文以多煤層開采深度作為評價指標進行分析,小康礦開采資料顯示,該形變區(qū)域?qū)猄2N5和S2N7綜采面的位置,多煤層開采深度分別對應為427～593 m和534～614 m之間.由表1可知,雙源模型反演深度為530.40和579.24 m;三源模型反演得到多煤層深度分別為570.93,596.55和574.42 m;兩種多源模型反演深度均與實際開采深度一致,表明了多源模型在該區(qū)域多煤層參數(shù)反演中的可靠性.

3.2 對比分析

基于雙源和三源Okada模型反演得到的多煤層參數(shù),分別對地表形變進行正演,模擬地表形變結(jié)果,如圖7和圖8所示，分析模擬形變和觀測形變之間的殘差(圖7c和圖8c),發(fā)現(xiàn)三源模型的殘差控制更好,符合多煤層開采引起的實際形變.為了進一步對比觀測形變、雙源模型和三源模型模擬形變之間的差異,本文通過比較剖線A1A2A3上的形變信息,評價不同反演模型的精度.圖9a中,紅色實點表示觀測形變,粉色虛線表示雙源模型模擬形變,藍色虛線表示三源模型模擬形變.可以看出兩種模型都反演得到兩個漏斗式的沉降中心,但三源模型擬合程度更好.圖9b顯示了雙源模型和三源模型的誤差分布,其中藍色實線表示的三源模型誤差在更靠近0的附近振蕩,相較于雙源模型,反演精度更高.

圖7 觀測形變、雙源模型模擬形變和殘差結(jié)果對比圖Fig.7 Comparison of the observed deformation,simulated deformation of dual-source model and residual results(a)—InSAR觀測圖; (b)—雙源模型模擬圖; (c)—雙源模型殘差圖.

圖8 觀測形變、三源模型模擬形變和殘差結(jié)果對比圖Fig.8 Comparison of the observed deformation, simulated deformation of three-source model and residual results(a)—InSAR觀測圖; (b)—三源模型模擬圖； (c)—三源模型殘差圖.

圖9 沿剖線A1A2A3的觀測形變和模擬形變對比分析圖Fig.9 Comparison of observed and modelled deformation along profile A1A2A3(a)—觀測值與模擬值; (b)—模型誤差.

為了定量評價不同反演模型的精度,本文根據(jù)公式(7),(8)分別計算兩種反演模型值與InSAR觀測值在剖線A2A3上的均方根誤差(root mean squared error,RMSE)和線性擬合優(yōu)度(R2).

(7)

(8)

計算雙源模型的RMSE為15.9 mm,R2為0.399 0;三源模型的RMSE為6.7 mm,R2為0.791 2.對比發(fā)現(xiàn)三源模型的均方根誤差更小,與InSAR觀測結(jié)果的相關(guān)性更高,證明三源Okada模型更符合該區(qū)域近距離多煤層共同采動的實際形態(tài).

4 結(jié) 論

1) 研究區(qū)內(nèi)主要分布有兩個沉降漏斗,其他區(qū)域無明顯沉降現(xiàn)象.經(jīng)實地調(diào)查驗證,發(fā)現(xiàn)沉降區(qū)對應于地下煤炭開采的位置,說明該區(qū)域地面沉降主要是由煤炭采空所致.在2018年1月～2019年8月間,小康礦采區(qū)最大累積沉降量為-194.8 mm,最大沉降速率達到-124.4 mm/a.大平礦采區(qū)最大累積沉降量為-136.4 mm,最大沉降速率為-99.7 mm/a.

2) 為了解決近距離多煤層參數(shù)的反演問題,本文提出了一種基于多源模型的近距離多煤層參數(shù)反演方法.該方法首先依據(jù)數(shù)據(jù)降采樣計算得到的較高位移梯度和數(shù)據(jù)邊緣信息估計有限矩形源數(shù)目,然后以疊加理論為原則,將多個矩形源疊加作為反演模型,最終通過非線性貝葉斯反演算法計算多源模型參數(shù)的后驗分布概率,選取后驗概率最大的一組多煤層參數(shù)作為最佳擬合參數(shù).

3) 將所提出的方法應用于近距離多煤層開采的小康礦進行驗證.通過定量計算,發(fā)現(xiàn)三源模型的RMSE相較于雙源模型減少了9.2 mm,R2達到了0.791 2.結(jié)果表明三源模型反演精度更高,更符合小康礦近距離多煤層共同采動的實際特征.