煤礦采空區覆巖“三帶”智能識別方法

李建文，趙 文，吳振坤，徐小兵，王慶濤，段隆臣

(1.中煤科工生態環境科技有限公司，北京 100013；2.中國地質大學(武漢) 自動化學院，湖北 武漢 430074；3.中國地質大學(武漢) 工程學院，湖北 武漢 430074)

煤礦采空區治理是踐行礦山綠色發展理念的重要舉措，而對采空區的覆巖裂隙發育情況進行勘探是開展采空區治理的基礎條件。煤層開采后由于巖層移動、變形和破壞，在覆巖中形成采動裂隙，導致礦區的地質結構遭到破壞，極易引起地表塌陷、裂縫和下沉等地質災害，嚴重影響城市居民生活與建設發展。因此，探明采空區地下裂隙分布及覆巖破壞情況，對重新利用采空區場地和滿足城市建設開發需求有重要意義。

“三帶”型是一種典型的采空區覆巖發育破壞類型。根據覆巖中的裂隙發育情況，自下而上可分為垮落帶、斷裂帶和彎曲帶，簡稱“三帶”。其中，斷裂帶和垮落帶合稱導水裂隙帶，在《煤礦采空區巖土工程勘察規范(2017 年版)》(簡稱勘察規范)中也稱垮落斷裂帶[1]。

目前確定“三帶”分布與高度的研究主要包括實測“三帶”裂隙發育程度和預測導水裂隙帶高度兩方面。其中，常見的實測手段包括工程物探[2-4]、工程鉆探[5]、煤田測井[6]、鉆孔成像[7-8]、井下注水觀測[9-10]等。實測法是最可靠的“三帶”識別方式，但總體成本較高，并需要專門的技術人員對數據進行分析才能得到采空區覆巖裂隙發育情況。預測導水裂隙帶高度的研究方法主要包括經驗公式法[11-13]、相似材料模擬法[14]、應力計算法[15-18](也指數值模擬法)，以及機器學習法[19-20]等。其中，勘察規范中的經驗公式應用較為廣泛，可以得到斷裂帶和垮落帶的發育高度范圍，但不能得到垮落帶和斷裂帶的具體發育高度和“三帶”的具體分布深度，需要結合實測方法才能確定[21]。通過預測導水裂隙帶高度確定“三帶”分布與高度的主要優點是應用成本低，對實際工程有一定指導意義，但大多基于室內試驗獲得，因而不能完全反映實際情況[22-23]。

隨著采空區治理自動化、智能化的發展，將物聯網技術、數據融合和分析技術進行有機結合，從而預測斷裂帶高度，為鉆進過程中自動獲取“三帶”高度提供了新的思路[24]。目前鉆探現場可以安裝各種類型的傳感器，自動記錄鉆進時間、鉆進深度及沖洗液漏失量等鉆進參數，但缺少自動分析數據的功能。技術人員完鉆后對數據進行統一分析時具有滯后性，且分析結果受技術人員主觀判斷影響顯著。此外，施工現場對采空區的初步判別需要技術人員持續關注鉆進參數的變化，而當出現采空區誤判或漏判時會則影響下一步的實時決策，例如針對可能出現的孔內事故采取預防措施和確定完鉆孔深等。

綜上所述，本文結合實測“三帶”裂隙發育程度和預測導水裂隙帶高度兩類研究方法的特點，基于工程鉆探中實時獲取的鉆進數據，采用貝葉斯變化點檢測(Bayesian Online Changepoint Detection，BOCD)[25]算法在鉆進過程中自動分析數據，并融合勘察規范中的經驗公式，實現在鉆進過程中對覆巖“三帶”的智能識別。該方法可充分發揮鉆進數據的時效性，提高“三帶”劃分效率和準確性，并同時提高煤礦采空區治理的智能化水平。

1 “三帶”的鉆進參數響應特征

圖1 為采空區覆巖“三帶”的剖面。在工程鉆探中，沖洗液漏失量和鉆速的變化特征與“三帶”有很大相關性[7,26-28]。

圖1 采空區覆巖“三帶”剖面[7]Fig.1 Schematic cross section of overlying strata three zones of goaf[7]

彎曲帶為斷裂帶頂部到地表之間的地層，基本呈整體移動，受采動影響不大。在垂直剖面上，地層上部很少出現離層裂隙；相比而言，地層下部更常見離層裂隙，但這些離層裂隙僅局部充水，且不與斷裂帶連通。因此，在彎曲帶中鉆進時，沖洗液漏失量和鉆速總體趨于平穩，局部存在輕微波動。

斷裂帶位于彎曲帶下方、垮落帶上方。其上部有裂隙，可以間接導水和積水，但因垂直裂隙不發育，一般不與下部裂隙溝通；其下部垂直裂隙逐漸發育增強，離層裂隙和垂向裂隙連通，導水性明顯增加，并能向下滲流至采空區。因此，鉆進至斷裂帶時，沖洗液漏失量會突然大幅增加，相較于正常地層可能增加數倍。對于垂直裂隙發育的地層，沖洗液甚至會全部漏失。與此同時，由于裂隙發育，鉆速會大幅提升，并存在較大波動。

垮落帶位于采空區覆巖最下部，通常為較大空洞或充填有垮落的巖塊，且巖塊之間存在很多空隙，連通性較強。因此，大多數“三帶”發育的采空區鉆孔在鉆進至垮落帶時，沖洗液漏失量很大，甚至全部漏失，而鉆速相較于正常地層也會顯著提高，且可能存在明顯波動。

當鉆孔進入采空區底板后，由于地層受采動影響不大、相對完整，鉆速相較于垮落帶會明顯降低，且逐漸趨于穩定。

2 “三帶”識別方法

通過“三帶”鉆進響應特征分析可知，從彎曲帶進入斷裂帶時沖洗液漏失量和鉆速均會明顯增加，而進入采空區底板時鉆速會明顯降低。因此，將這兩個特征的提取看作對變化點的識別，并采用BOCD 算法來識別彎曲帶下限和采空區底板。

然而，僅通過變化點檢測無法獲得實際的“三帶”界限，主要原因為：(1)識別變化點與目標變化點為多對一，需要從多個變化點中篩選出目標變化點；(2)通過鉆進參數僅能確定彎曲帶下限深度、采空區底板深度和垮落帶下限深度，無法確定斷裂帶下限深度。因此，本文在BOCD 算法識別變化點的基礎上，結合勘察規范中有關斷裂帶和垮落帶高度的經驗公式，以及相關鉆前資料，以此自動識別準確的“三帶”界限，并計算得到“三帶”分布的高度值。

2.1 BOCD 算法原理

BOCD 算法基于貝葉斯框架，可以實時識別數據隨時間的變化情況，目前已被廣泛應用于金融分析、生物識別等領域。BOCD 算法的核心思想是：假設時間序列x1:t(x1:t為序列x1,x2,···,xt的集合，t∈N*為序列中數據點所對應的時間或序號)中任意兩個相鄰變化點之間的數據相互獨立，且分布類型相同，則新觀測數據xt+1是否為變化點的判斷依據是該數據的預測條件概率值P(xt+1|x1:t)的大小。若新觀測的數據與上一變化點之后的數據非常相似，則其預測條件概率值將會很高，即新的觀測數據不是變化點。反之，其預測條件概率值會很低，即新的觀測點為變化點。

定義行程長度rt為從上一變化點到t時刻走過了r個單位時間，表示t時刻行程長度為r時的數據集合，則P(xt+1|x1:t)的表達式為：

將風險函數f看作離散指數分布，即f=1/λ，其中λ為尺度參數，與數據的采樣頻率有關。P(rt,x1:t)可以通過上一時刻的概率迭代計算得到。

2.2 “三帶”智能識別方案

假設待識別鉆孔的設計孔深為d，垮落帶、斷裂帶分布在設計孔深往上Hu之內，候選變化點的有效識別深度區間為[d-Hu,d]，所有“三帶”界限深度均位于該區間內。Hu可按如下公式計算：

Hli與采空區頂板傾角α、采空區頂板巖層的碎脹系數k等有關[1,7]，該公式給出了Hli的最大和最小值，為了擴大候選變化點的有效深度區間，盡量覆蓋更多的候選點，這里取Hli的最大值；與Hli的取值類似，這里取最大值，取最小值；M為已知值，Hn為設計值。

由“三帶”鉆進中的沖洗液漏失量和鉆速變化特征可知，區間[d-Hu,d]內的第一個變化點所對應的深度即為彎曲帶下限深度d1；最后一個變化點所對應的深度即為采空區底板深度d4，則垮落帶下限深度d3為d4-M。

上述界限深度確定后，接下來需要確定斷裂帶下限深度d2。由于斷裂帶下限深度附近可能存在多個候選變化點，應用垮落斷裂帶和垮落帶高度的經驗公式[1]，結合彎曲帶下限深度d1和垮落帶下限深度d3，由下式計算斷裂帶下限深度d2的識別深度范圍[d2]：

式(5)中，括號內的值為估算的Hm變化比例范圍，正負號分別代表范圍的上下限；尋找距離[d2]區間中心最近的變化點，即為斷裂帶下限深度d2；若[d2]區間內無變化點，則取該區間的中心點為斷裂帶下限深度d2。

根據上述“三帶”界限識別結果，可得斷裂帶實際高度為d2-d1、垮落帶實際高度為d3-d2。至此，“三帶”的界限深度和高度已經全部識別完成，完整的識別流程如圖2 所示。

圖2 “三帶”識別流程Fig.2 Three zones identification flow chart

按照圖2 流程，彎曲帶下限深度、采空區底板深度可以實時確定。在通過采空區底板確定垮落帶下限后，結合勘察規范中經驗公式的結果可以確定斷裂帶下限深度。當采空區底板深度確定后，即可確定完鉆深度。在實際工程中，該流程中涉及到的鉆前地質資料等參數的值(如M、d、k、α、Hn等)可以在鉆進前通過工控機或物聯網設備輸入給定，鉆進時無需其他輔助即可自動識別“三帶”結果。需要指出的是，上述“三帶”識別方法是基于回轉鉆進的鉆進參數特征實現的，如果鉆進方式發生明顯改變，則其鉆進過程中的“三帶”識別還需要進一步開展相關適配性研究。

上述BOCD 算法和“三帶”識別過程可以通過Python 代碼實現，并將其作為Web 應用部署在服務器上，從而通過無線網絡分別與數據感知端和用戶應用端相連，實現“三帶”識別的物聯應用。“三帶”識別的Web 應用如圖3 所示，數據感知端包含不同的鉆井監測站點，每個站點由各類傳感器(如鉆速傳感器、流量傳感器等)組成，用于獲取基礎數據。Web 服務器是主要功能的實現載體，如數據存儲功能、“三帶”識別功能等，其通過無線網絡接收數據感知端的數據和將識別結果等數據傳輸到應用終端。應用終端直接面向用戶，包括移動終端、大屏等，用于顯示識別結果或輸入數據。

圖3 “三帶”識別的Web 應用Fig.3 Illustrating the application of three zones identification for goaf in Web applications

3 工程應用

3.1 項目概況

為了驗證所提方法的應用效果，選擇山東某礦區一口勘察井的監測數據進行測試。礦區采用立井分水平開拓(主井、副井、風井)，中央并列式通風，條帶短壁式綜合機械化一次采全高采煤工藝，頂板自由垮落法管理[7]。測試鉆孔的采空區頂板為砂質泥巖、泥巖，單軸飽和抗壓強度平均值分別約為40.77 和31.66 MPa，屬于較硬巖。煤層平均傾角α=6°，頂板巖層的碎脹系數為k=1.35，平均采高為M=2.8 m。鉆孔為垂直鉆孔，設計孔深d=570.00 m。

該項目所在勘察區使用回轉鉆進工藝，并采用PDC 鉆頭作為碎巖工具。測試鉆孔只統計純鉆進時間內的沖洗液漏失量和鉆速，其中沖洗液漏失量通過流量計測量進漿量和出漿量得到，深度和鉆速通過無線射頻位移傳感器得到。從鉆進孔深為266.00 m 時開始，每鉆進約0.5 m 記錄一次數據，鉆進至孔深565.87 m 時停止記錄，共657 組數據。

圖4 為現場獲取的沖洗液漏失量和鉆速數據所繪制的曲線，兩條曲線均可以劃分為兩段。在第一段曲線中，沖洗液漏失量整體數值在110 L/m 左右，部分數據的值較高，主要是由于監測誤差和這一段地層中存在的破碎帶造成；鉆速曲線整體數值在11 mm/min 左右，部分數據的值較低，考慮主要是由于監測誤差造成。在第二段曲線中，沖洗液漏失量和鉆速值都明顯高于第一段曲線，這主要與鉆頭進入垮落帶和斷裂帶有關。

圖4 原始鉆進監測數據Fig.4 Raw drilling monitoring data

3.2 BOCD 變化點識別

利用鉆進數據識別“三帶”時，同一采樣頻率的λ大致相同，且采樣頻率越大，λ也越大。實際應用中可通過多次測試確定某一采樣頻率對應的λ值，本文中的λ取值為50。給定上述條件，即可以識別鉆進數據的變化點。

在識別前，為了提高識別效果，選擇中值濾波的方法降低異常數據對變化點識別結果的干擾。中值濾波窗口大小設置為5。圖5 為BOCD 算法的運行結果，每條垂直虛線都與一個行程的開始對應，即對應一個變化點，其中圖5a 和圖5b 分別為沖洗液漏失量和鉆速濾波后的曲線及對應的變化點識別結果。與原始數據相比，濾波后沖洗液漏失量和鉆速的第一段曲線更加平穩，基本沒有異常值，第二段曲線同樣波動減小，并且曲線的變化特點更加明顯，表明濾波可以有效地去除異常值。如圖5 所示，沖洗液漏失量曲線識別出來的變化點對應深度為494.22、511.36、515.36、543.50 m，鉆速曲線識別出來的變化點對應深度為269.72、505.68、510.42、515.71、544.26、562.87 m，所有這些深度值都是用來確定“三帶”界限的候選變化點。

圖5 濾波后鉆進監測數據變化點識別結果Fig.5 Identification results of Bayesian Online Changepoint Detection from the drilling monitoring data after filtering

3.3 “三帶”識別結果與分析

試驗鉆孔所在地區已探明“三帶”的鉆孔中，已知實際垮落斷裂帶高度的最大值為87.30 m，設計鉆孔要求鉆入采空區底板深度Hn不少于10 m，采高M為2.8 m。

由給定的鉆前地質資料可知，采空區頂板屬于較硬巖，則根據勘察規范附錄L 可以確定試驗鉆孔垮落斷裂帶高度Hli、垮落帶高度Hm的經驗計算公式如下：

取 max(Hli)為61.14 m。本實例中，區域內只有厚度為2.8m 的單一煤層，因此根據規范計算獲得的Hli和相同。取min為43.34m，由式(4)可得Hu為135.95 m，則候選變化點的有效深度區間為444.04～570.00 m。這一范圍內的第一個變化點所對應深度為494.22 m，則彎曲帶下限深度d1為494.22 m。這一范圍內的最后一個變化點對應深度為562.87 m，則采空區底板深度d4為562.87 m。進而確定垮落帶下限深度d3為560.07 m。

由式(5)可得d2的識別深度范圍為[550.17，552.65]，該界限內沒有變化點，則斷裂帶下限深度d2取識別區間中心點所在位置，即551.41 m。

至此，已得到“三帶”的下限深度自上而下分別為494.22、551.41、560.07 m。為評價本文識別方法的準確性，選擇與本文依托相同項目和測試鉆孔的文獻[7]中的“三帶”結果作為對比。該文獻綜合利用工程地質鉆探、孔內電視與煤田測井3 種方法，得到測試鉆孔的“三帶”界限深度自上而下分別為493.55、551.10、559.55 m。圖6 為智能識別得到的“三帶”結果(灰色垂直虛線)與文獻[7]所得結果(黑色垂直虛線)的對比情況，可以看出兩者較吻合。“三帶”界限深度的智能識別結果與實際結果對比見表1，誤差均小于1 m。“三帶”高度的智能識別結果和經驗公式計算結果與其實際高度的對比見表2。其中，智能識別結果的高度誤差均小于3%，而經驗公式計算得到的斷裂帶和垮落高度誤差分別為9.23%和4.85%。結果表明，“三帶”智能識別方法顯著優于經驗計算方式，且具有較高的識別精度。

表1 識別“三帶”界限與實際“三帶”界限[7]對比Table 1 Comparison of the intelligent identification results and actual results[7] boundary depths in the three zones of goaf

表2 “三帶”高度結果對比Table 2 Comparison of heights in the three zones of goaf

圖6 “三帶”識別結果與實際人工劃分結果[7]對比Fig.6 Comparison of the intelligent identification results and actual results[7] in the three zones of goal

應用結果表明，在變化點識別的基礎上，結合勘察規范中垮落斷裂帶和垮落帶的經驗公式，可將斷裂帶下限的范圍由彎曲帶下限和垮落帶下限之間縮到更小的范圍，使得識別結果更有針對性。因此，融合了BOCD算法的采空區“三帶”智能識別方法比傳統方法更接近實際情況。

4 結論

a.提出的基于沖洗液漏失量和鉆速等鉆進數據的變化點檢測方法，可以實現在鉆進過程中的“三帶”智能識別，無需等完鉆后再進行數據分析，充分發揮了鉆進數據的時效性，有效避免了技術人員的主觀判斷和經驗對結果產生的影響。

b.與實際多種勘察手段綜合劃分的結果對比表明，智能識別的“三帶”深度誤差均小于1 m，高度誤差均小于3%，具有較高的識別精度。將實測方法和預測方法相結合，將沖洗液漏失量和鉆速這兩個與“三帶”有密切關系的數據作為輸入參數，提高了識別的精度。提出的“三帶”智能識別方法可以代替多種手段綜合劃分的方式，有助于提高“三帶”劃分效率，降低采空區勘察成本。

c.提出的“三帶”智能識別方法適用于開采方式、地質條件、煤層厚度、底板位置等相關信息比較明確的采空區勘探。現有監測參數僅考慮了鉆速和沖洗液漏失量，且其響應特征主要針對以沖洗液為循環介質的回轉鉆進，未來可進一步研究其他鉆進參數及其在不同鉆進方式中與“三帶”之間的響應關系，以提高智能識別方法的適用性。

符號注釋：

d為設計孔深，m；d0為地表深度，m；d1為彎曲帶下限深度，m；d2為斷裂帶下限深度，m；d3為垮落帶下限深度，m；d4為采空區底板深度，m；f為風險函數；Hn為待識別鉆孔規定鉆進采空區底板的深度，m；Hli為按勘察規范經驗公式計算得到的待識別鉆孔垮落斷裂帶高度，m；為區域內已探明的實際垮落斷裂帶高度，m；為按勘察規范經驗公式計算得到的已完成鉆孔垮落斷裂帶高度，m；Hm為按勘察規范中經驗公式計算得到的待識別鉆孔垮落帶高度，m；Hu為待識別鉆孔彎曲帶下限與設計孔深的最大距離，m；k為采空區頂板巖層的碎脹系數；max(Hli)、min(Hli)、ave(Hli)分別為按經驗公式計算得到的待識別鉆孔垮落斷裂帶的最大高度、最小高度和平均高度，m；為實際垮落斷裂帶高度的最大值，m；min()為按勘察規范中經驗公式計算得到的已完成鉆孔垮落斷裂帶高度的最小值，m；M為采高，m；N*為正整數集；P(rt|rt-1)為變化點先驗概率；rt-1、rt為行程長度，分別表示從上一變化點到t-1、t時刻走過的單位時間；t為觀測數據點所對應的時間或序號，從1 開始；x1:t、x1:t-1分別為從第一個時間單位到第t、t-1 個時間單位的觀測數據集；xt+1、xt分別為t+1、t時刻的觀測數據；為t時刻對應行程長度為r內的數據集合；α為采空區頂板傾角，(°)；λ為離散指數分布的尺度參數。