石曉宇， 張燕海， 楊可明， 姚樹一， 王劍

(1.中國礦業(yè)大學(xué)(北京) 地球科學(xué)與測繪工程學(xué)院， 北京 100083；2.淮北礦業(yè)(集團(tuán))股份有限公司 通防地測部， 安徽 淮北 235000)

0 引言

煤炭是目前我國主要的供給能源[1]。地下煤炭資源的大面積開采使得采區(qū)周圍巖層原始應(yīng)力平衡被打破[2]，引起巖層移動和變形，不僅直接造成礦區(qū)地表沉陷，還會帶來一系列潛在的災(zāi)害問題。對礦區(qū)地表沉陷規(guī)律的研究是采礦工程中一個重大課題，其中如何更高效、更系統(tǒng)、更安全地對礦區(qū)地表進(jìn)行沉陷監(jiān)測與預(yù)計備受關(guān)注。

地下采煤引起的巖層移動和變形是一個復(fù)雜的時空過程，傳統(tǒng)的地表沉陷預(yù)計方法(如概率積分法、剖面函數(shù)法等)僅能獲得地表沉陷最終穩(wěn)定后的空間變化，不能較好地反映開采沉陷隨時間變化的動態(tài)過程[3]。眾多學(xué)者綜合考慮因巖體本身復(fù)雜性而難以全面研究采煤工作面上覆巖層移動及地表沉陷動態(tài)過程的情況，以局部觀測點(diǎn)為切入點(diǎn)，建立沉陷區(qū)單點(diǎn)沉陷量時間序列模型，從而研究礦區(qū)地表動態(tài)沉陷規(guī)律[4-6]。GM(1,1)與灰色Verhulst模型為灰色系統(tǒng)理論中2個基礎(chǔ)的時間序列模型，可在小樣本、貧信息、不確定性等劣勢情況下，鑒別系統(tǒng)因素間發(fā)展趨勢的相異程度，尋找系統(tǒng)變化規(guī)律，并對事物發(fā)展趨勢作出模糊性的短期、中長期預(yù)測[7-8]。GM(1,1)是用一階微分方程對1個變量建立的模型[9]，主要解決具有指數(shù)變化規(guī)律的生成序列，只能描述單調(diào)變化過程；灰色Verhulst模型是由德國生物學(xué)家Verhulst在研究物種繁殖趨勢時提出的一階非線性動態(tài)模型[10]，其基本思想是生物繁殖數(shù)量受到周圍環(huán)境的影響，繁殖曲線會以“S”形呈現(xiàn)，最終達(dá)到一個平衡穩(wěn)定的狀態(tài)。將這2種灰色模型(Grey Model，GM)應(yīng)用于礦區(qū)地表沉陷預(yù)計，建立沉陷量時間序列模型，可較好地揭示地表單點(diǎn)連續(xù)動態(tài)變化的過程，實(shí)現(xiàn)地表沉陷動態(tài)預(yù)計。然而，現(xiàn)有文獻(xiàn)只針對2種GM在礦區(qū)沉陷預(yù)計方面的應(yīng)用進(jìn)行研究，并沒有結(jié)合模型自身特點(diǎn)，分析和比較二者在礦區(qū)地表沉陷預(yù)計中的適用性。此外，目前礦區(qū)地表單點(diǎn)動態(tài)沉陷預(yù)計方法主要基于傳統(tǒng)的水準(zhǔn)測量數(shù)據(jù)[11]，監(jiān)測方法單一，成本高，且觀測點(diǎn)易破壞，不能保證地表形變信息的實(shí)時性。合成孔徑雷達(dá)差分干涉測量(Differential Interferometry Synthetic Aperture Radar，D-InSAR)技術(shù)通過雷達(dá)干涉圖的差分獲得地表形變信息，以高精度、全天候、連續(xù)空間覆蓋、實(shí)時性等優(yōu)勢成為新興的形變監(jiān)測手段，被廣泛應(yīng)用于礦區(qū)地表沉陷監(jiān)測中[12-15]。將D-InSAR技術(shù)與GM相結(jié)合，避免監(jiān)測與預(yù)計獨(dú)立進(jìn)行，可較好地形成系統(tǒng)化及模塊化流程，實(shí)現(xiàn)礦區(qū)地表沉陷監(jiān)測與動態(tài)預(yù)計一體化。

本文以淮北礦業(yè)(集團(tuán))股份有限公司袁店二礦7221工作面為試驗(yàn)區(qū)域，融合D-InSAR技術(shù)和GM(1,1)與灰色Verhulst模型，通過建立沉陷量與時間的GM進(jìn)行地表沉陷動態(tài)預(yù)計，實(shí)現(xiàn)監(jiān)測與預(yù)計一體化，并針對2種GM的特點(diǎn)，結(jié)合預(yù)計結(jié)果分析和比較其在礦區(qū)地表沉陷預(yù)計中的適用性。

1 礦區(qū)地表沉陷D-InSAR監(jiān)測

1.1 試驗(yàn)區(qū)域與相關(guān)數(shù)據(jù)

袁店二礦7221工作面上方地表主要覆蓋物為村莊和農(nóng)田，四鄰無其他采掘活動。該工作面為俯斜開采，工作面寬約110 m，走向長約544 m，煤層傾角為10～20°，基巖厚度為132～162 m，采高為4.5 m，于2017-12-06開始回采，2018-06-03停采。為了控制地表采動損害在建筑物等安全保護(hù)等級范圍之內(nèi)，分別在地表公路或建筑體上設(shè)置水準(zhǔn)測量與變形監(jiān)測點(diǎn)，以監(jiān)測地表及重要建筑物變形情況。但該區(qū)域內(nèi)人工溝渠縱橫分布，監(jiān)測點(diǎn)稀少、空間配置不合理且常被人為破壞，無法準(zhǔn)確對重要建筑物進(jìn)行損害評估。因此，有必要在工作面開采過程中對該區(qū)域進(jìn)行地表沉陷D-InSAR監(jiān)測。

為了盡可能減少D-InSAR技術(shù)相位失相干對地表沉陷監(jiān)測結(jié)果精度的影響，選取冬季和初春時段(該時段植被蕭疏且地物波譜特征無明顯變化)的9景哨兵衛(wèi)星(Sentinel-1A)雷達(dá)影像作為D-InSAR監(jiān)測的數(shù)據(jù)源。影像時間為2017-11-04—2018-02-08，相鄰時間間隔(12 d)的影像可以形成8個干涉對，見表1，其中VV為同向垂直極化方式，IW為干涉測量寬幅模式。外部數(shù)字高程模型選用分辨率為90 m的SRTM(Shuttle Radar Topography Mission，航天飛機(jī)雷達(dá)地形測繪使命)數(shù)據(jù)。

表1 Sentinel-1A雷達(dá)影像干涉對

1.2 D-InSAR監(jiān)測結(jié)果及分析

根據(jù)表1數(shù)據(jù)，利用D-InSAR技術(shù)二軌法獲取8組7221工作面時間序列沉陷分布。為方便分析，以2017-11-04雷達(dá)影像為參考對象，將各組沉陷量依次累加，即可得到各組時間序列沉陷分布。由于2017-11-04—28地形無明顯變化，所以僅給出后6組時間序列沉陷分布，如圖1所示?？煽闯鲈诨夭沙跗?2017-12-06—10)，沉陷量為-10～0 mm。隨著工作面推進(jìn)，進(jìn)入沉陷發(fā)展階段(圖1(b)—圖1(f))，在2017-12-22開始形成下沉盆地，與7221工作面推進(jìn)方向一致。下沉盆地自西南向東北方向擴(kuò)張，沉陷量和平面范圍隨之增大。至2018-02-08，最大沉陷量為-85 mm，位于7221工作面切眼附近，離注漿孔(注3、注4、注5，位于工作面東北部)較近的村莊未發(fā)生明顯沉陷，說明注漿充填開采工藝有很好的局域減沉效果。

(a) 2017-11-04—2017-12-10 (b) 2017-11-04—2017-12-22 (c) 2017-11-04—2018-01-03

(d) 2017-11-04—2018-01-15 (e) 2017-11-04—2018-01-27 (f) 2017-11-04—2018-02-08

圖1 7221工作面時間序列沉陷分布

Fig.1 Subsidence distribution of 7221 working face in time series

為更好地對D-InSAR監(jiān)測結(jié)果進(jìn)行定量研究，提取下沉盆地中沿7221工作面傾向l1和走向l2(l1，l2位置如圖2所示)的沉陷剖面線，如圖3所示。l1長807 m，方向?yàn)樽晕飨驏|，設(shè)其最西端為零點(diǎn)；l2長1 089 m，方向?yàn)樽阅舷虮?，設(shè)其最南端為零點(diǎn)。從圖3(a)可看出，2017-11-04—2018-02-08，l1最大沉陷量為-81 mm，2017-11-28—2018-01-03沉陷剖面線變陡，沉陷量顯著增大，后期進(jìn)入沉陷衰減階段，開始微小沉陷。從圖3(b)可看出，2017-11-04—2018-02-08，l2最大沉陷量為-84 mm，2017-11-28—2018-01-03沉陷速度較快，動態(tài)最大沉陷量不斷增加，2018-01-03—02-08沉陷速度開始衰減，但平面影響范圍不斷增大。對比圖3(a)和圖3(b)可看出，開采對傾向沉陷影響范圍小于走向，且傾向沉陷階段發(fā)育早于走向，下沉盆地近似出現(xiàn)平底狀態(tài)。

2 GM(1,1)與灰色Verhulst模型動態(tài)預(yù)計

2.1 模型建立與預(yù)計

基于D-InSAR監(jiān)測獲取的2017-11-16—2018-02-08 7221工作面地表沉陷量，從各期沉陷分布(以表1中干涉對組合為參考， 2017-11-16—28為第1期觀測數(shù)據(jù)，2017-11-16—12-10為第2期觀測數(shù)據(jù)，以此類推至2018-02-08，共7期數(shù)據(jù))中提取43個水準(zhǔn)點(diǎn)的像元信息，并繪制沉陷量曲線，如圖4所示。選取H2，H3，E2，G5，F(xiàn)4為觀測點(diǎn)，對7期數(shù)據(jù)建立GM(1,1)與灰色Verhulst模型方程，獲得前7期沉陷量擬合值,并預(yù)計第8，9期(2018-02-20，2018-03-04)沉陷量。1—7期各觀測點(diǎn)沉陷實(shí)際值與2種GM擬合值對比如圖5所示。

圖2 l1，l2位置

(a) 傾向l1

(b) 走向l2

2.2 預(yù)計結(jié)果精度檢驗(yàn)

GM精度檢驗(yàn)一般分為事中檢驗(yàn)和事后檢驗(yàn)[16]。事中檢驗(yàn)包含殘差檢驗(yàn)、關(guān)聯(lián)度檢驗(yàn)和后驗(yàn)差檢驗(yàn)。殘差檢驗(yàn)屬于算術(shù)檢驗(yàn)，定義為模型值與實(shí)際值之間的差，是一種較為客觀的逐點(diǎn)檢驗(yàn)法；關(guān)聯(lián)度檢驗(yàn)屬于幾何檢驗(yàn)，表征2個事物之間的關(guān)聯(lián)程度，關(guān)聯(lián)值越大(越接近1)則相關(guān)性越好，模型等級越高；后驗(yàn)差檢驗(yàn)屬于統(tǒng)計檢驗(yàn)，是按照均方差比值C和小誤差概率P來評估灰色模型優(yōu)劣，C越小越好，相反P越大越好。事后檢驗(yàn)即預(yù)測檢驗(yàn)，通常采用殘差檢驗(yàn)和相對誤差檢驗(yàn)來了解預(yù)測誤差。本文對于GM預(yù)計結(jié)果的精度檢驗(yàn)分為2步：采用殘差檢驗(yàn)、關(guān)聯(lián)度檢驗(yàn)和后驗(yàn)差檢驗(yàn)對GM(1,1)和灰色Verhulst模型進(jìn)行擬合精度檢驗(yàn)，結(jié)果見表2；對預(yù)計的第8,9期沉陷值進(jìn)行殘差和相對誤差檢驗(yàn)，結(jié)果見表3和表4。

(a) 水準(zhǔn)點(diǎn)位置

(b) 各水準(zhǔn)點(diǎn)沉陷量曲線

(a) H2

(b) H3

(c) E2

(d) G5

(e) F4

表2 GM(1,1)和灰色Verhulst模型擬合精度檢驗(yàn)Table 2 Fitting precision checking of GM(1,1) and grey Verhulst model

表3 第8期沉陷量預(yù)計精度檢驗(yàn)Table 3 Prediction precision checking of subsidence value in phase 8

2.3 預(yù)計結(jié)果分析

(1) 從圖5(a)、圖5(b)可看出，實(shí)際沉陷量曲線后期呈平底飽和狀態(tài)，此時GM(1,1)和灰色Verhulst模型前期(1,2,3期)擬合結(jié)果差別不大，但后期(4—7期)灰色Verhulst模型較GM(1,1)收斂快，擬合結(jié)果與實(shí)際沉陷量吻合較好。圖5(c)中實(shí)際沉陷量曲線后期包含非單調(diào)擺動過程，使用GM(1,1)擬合誤差較大。圖5(d)、圖5(e)中實(shí)際沉陷量曲線呈近似單峰型，GM(1,1)擬合曲線更接近實(shí)際值，且后期GM(1,1)和灰色Verhulst模型擬合結(jié)果相差不大。

表4 第9期沉陷量預(yù)計精度檢驗(yàn)Table 4 Prediction precision checking of subsidence value in phase 9

(2) 從表2可看出，觀測點(diǎn)H2,H3使用灰色Verhulst模型擬合結(jié)果明顯優(yōu)于GM(1,1)，平均殘差分別減小了約26.5%,36.7%。根據(jù)后驗(yàn)差精度等級(表5)，觀測點(diǎn)H2的GM(1,1)擬合結(jié)果只能達(dá)到合格，而其灰色Verhulst模型擬合結(jié)果精度等級達(dá)到良好；觀測點(diǎn)G5采用GM(1,1)擬合時平均殘差較灰色Verhulst模型減小了1.249 mm，關(guān)聯(lián)度提高了約3.2%；觀測點(diǎn)E2，F(xiàn)4采用2種GM擬合的效果相近，精度等級均達(dá)到良好，且平均殘差之差均小于1 mm。

表5 后驗(yàn)差精度等級Table 5 Precision level of posterior error

(3) 從表3、表4可看出，采用GM(1,1)和灰色Verhulst模型預(yù)計第8，9期5個觀測點(diǎn)的沉陷量時精度與表2中擬合精度較一致，證明可根據(jù)各點(diǎn)擬合精度(殘差、關(guān)聯(lián)度和后驗(yàn)差)選擇沉陷預(yù)計適用模型。GM(1,1)對觀測點(diǎn)H3的預(yù)計誤差最大，最大殘差絕對值為56 mm，顯然該點(diǎn)沉陷量曲線更適宜選用灰色Verhulst模型來描述。灰色Verhulst模型對第8期觀測點(diǎn)G5的預(yù)計精度最低，最大相對誤差為21%，而使用GM(1,1)預(yù)計時殘差減少了6 mm，且相對誤差減小了5%，預(yù)計效果更佳。

可見對于單點(diǎn)沉陷曲線后期呈平底飽和狀態(tài)的觀測點(diǎn)(如H2,H3)，地表進(jìn)入沉陷衰退階段時選用灰色Verhulst模型進(jìn)行沉陷動態(tài)預(yù)計，精度明顯優(yōu)于GM(1,1)，最大相對誤差減少了55%；地表進(jìn)入沉陷發(fā)展階段后，對于單點(diǎn)沉陷量曲線呈近似單峰型的觀測點(diǎn)(如G5，F(xiàn)4)，選用GM(1,1)預(yù)計沉陷量更合適，最大相對誤差可減少35%。

3 結(jié)論

(1) 井下采煤地表移動是一個先加速后減速的過程，即隨著工作面推進(jìn)和開采時間推移，地表點(diǎn)會經(jīng)歷開始移動、劇烈移動、緩慢移動及停止移動4個階段，其移動軌跡是一個隨機(jī)變化的過程。

(2) 礦山開采沉陷開始至活躍前期，地表單點(diǎn)沉陷量曲線呈近似單峰型，采用GM(1,1)進(jìn)行短期預(yù)計效果更佳；當(dāng)?shù)V區(qū)地表單點(diǎn)沉陷量曲線呈平底飽和狀態(tài)，沉陷進(jìn)入衰退階段，此時采用灰色Verhulst模型預(yù)計的結(jié)果與實(shí)際值誤差較小，適用于中長期預(yù)計。

(3) GM(1,1)和灰色Verhulst模型較傳統(tǒng)地表沉陷預(yù)計方法而言，是假設(shè)地表沉陷歷程符合已知時間序列所描述的動態(tài)過程，通過先驗(yàn)信息來研究或預(yù)測未知領(lǐng)域的時間函數(shù)模型，而開采引起的地表沉陷是一個復(fù)雜的時空變化過程，僅依靠時間序列模型并不能實(shí)現(xiàn)準(zhǔn)確預(yù)計，后續(xù)會將時間序列模型與某種常規(guī)靜態(tài)方法(如概率積分法)相結(jié)合，建立更精確的地表任意點(diǎn)動態(tài)沉陷預(yù)計模型。