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

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

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

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

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

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

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

近距離多煤層共同開采是目前全球范圍內(nèi)各礦場在實際開采作業(yè)過程中所面臨的一種復(fù)雜情況.研究其開采規(guī)律,并根據(jù)實際監(jiān)測的地表形變快速準(zhǔn)確反演地下復(fù)雜的多煤層幾何形狀、位置、產(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在地表的投影坐標(biāo)(X,Y)、開采面沿著走向的長度L、沿傾向的寬度W、煤層開采中心深度depth、走向角strike(開采掘進(jìn)方向與正北方向的夾角)、傾向角dip(開采面與水平方向的夾角)、rake(斷層破裂時相對于斷層走向移動的角度)、破裂方向的位錯分量slip、張性分量opening;圖1c是依托多煤層的分布情況,將Okada有限矩形源模型應(yīng)用于多煤層的綜合圖.

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

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

L(t)=

(3)

在抽樣計算多源模型后驗概率時,采樣方法的選取對于后驗概率的計算結(jié)果和效率至關(guān)重要.本文采用基于Metropolis-Hastings接受準(zhǔn)則的馬爾科夫鏈蒙特卡羅(Markov chain Monte Carlo,MCMC)算法[14]在空間上進(jìn)行采樣,該方法可以大大提高采樣接受率,彌補(bǔ)了傳統(tǒng)蒙特卡羅方法的抽樣分布不能隨著模擬的進(jìn)行而動態(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ù)的選擇指標(biāo).在尋找最優(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ù)進(jìn)行形變提取.SAR衛(wèi)星的航向角為-166.19°,入射角為37.13°.DEM數(shù)據(jù)選用德國TanDEM-X獲取的空間分辨率為3 rad·s的DEM,用于影像配準(zhǔn)、平地效應(yīng)去除、軌道精煉等.

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

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

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

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

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

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

基于SBAS-InSAR技術(shù),獲取了研究區(qū)在雷達(dá)視線向的累積沉降量,結(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)定,地面沉降主要是由煤炭開采所致.

為進(jìn)一步分析兩個沉降區(qū)不同位置的形變過程,分別選取了Q1和Q2兩個特征點,刻畫其沉降時間序列,結(jié)果如圖6所示.藍(lán)色實線表示大平礦采區(qū)特征點Q1的沉降時間序列,紅色實線表示小康礦采區(qū)特征點Q2的沉降時間序列.可以發(fā)現(xiàn)2018年1月～2019年8月間,大平礦沉降區(qū)的Q1點累積沉降量為-111.7 mm;小康礦沉降區(qū)Q2點的累積沉降量最大,達(dá)到-178.4 mm.持續(xù)的地下煤炭開采已造成康平煤田區(qū)域產(chǎn)生了巨大的地表形變梯度和影響范圍,使當(dāng)?shù)鼐用褡≌仍O(shè)施被積水淹沒或產(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ù)為基礎(chǔ),綜合考慮形變反演過程中數(shù)據(jù)點位密度對反演計算效率的影響,選用1.44×10-4作為自適應(yīng)四叉樹方差閾值進(jìn)行數(shù)據(jù)降采樣,然后通過較高位移梯度和數(shù)據(jù)邊緣信息,選取2和3作為多源模型的模型數(shù),分別使用非線性貝葉斯反演算法對近距離多煤層參數(shù)進(jìn)行反演.Okada模型中slip表示斷層破裂方向的位錯分量,rake表示斷層破裂時相對于斷層走向移動的角度,opening表示斷層垂直方向運動的張性分量.通過分析小康礦形變區(qū)的地表沉降特征,發(fā)現(xiàn)小康礦開采過程中多煤層傾角很緩,近乎水平開采,且地表沉降主要集中在垂直方向,因此本文簡化了Okada模型,假設(shè)斷層傾角為0°,只分析多煤層開采導(dǎo)致垂直方向的張性破裂引發(fā)的地表形變.在反演過程中使用獨立的測量坐標(biāo)系,參考基準(zhǔn)為[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ù)反演方法的可靠性.本文以多煤層開采深度作為評價指標(biāo)進(jìn)行分析,小康礦開采資料顯示,該形變區(qū)域?qū)?yīng)S2N5和S2N7綜采面的位置,多煤層開采深度分別對應(yīng)為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ù),分別對地表形變進(jìn)行正演,模擬地表形變結(jié)果,如圖7和圖8所示，分析模擬形變和觀測形變之間的殘差(圖7c和圖8c),發(fā)現(xiàn)三源模型的殘差控制更好,符合多煤層開采引起的實際形變.為了進(jìn)一步對比觀測形變、雙源模型和三源模型模擬形變之間的差異,本文通過比較剖線A1A2A3上的形變信息,評價不同反演模型的精度.圖9a中,紅色實點表示觀測形變,粉色虛線表示雙源模型模擬形變,藍(lán)色虛線表示三源模型模擬形變.可以看出兩種模型都反演得到兩個漏斗式的沉降中心,但三源模型擬合程度更好.圖9b顯示了雙源模型和三源模型的誤差分布,其中藍(lán)色實線表示的三源模型誤差在更靠近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ū)對應(yīng)于地下煤炭開采的位置,說明該區(qū)域地面沉降主要是由煤炭采空所致.在2018年1月～2019年8月間,小康礦采區(qū)最大累積沉降量為-194.8 mm,最大沉降速率達(dá)到-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) 將所提出的方法應(yīng)用于近距離多煤層開采的小康礦進(jìn)行驗證.通過定量計算,發(fā)現(xiàn)三源模型的RMSE相較于雙源模型減少了9.2 mm,R2達(dá)到了0.791 2.結(jié)果表明三源模型反演精度更高,更符合小康礦近距離多煤層共同采動的實際特征.