基于模糊數學的露天鐵礦采空區地球物理探測評價研究

賈三石 付建飛 王恩德 郭 凱 門業凱

（1.東北大學秦皇島分校資源與材料學院，河北秦皇島066004；2.東北大學資源與土木工程學院，遼寧沈陽110819）

采空區作為一種礦產資源開采后的遺留產物，大量存在于金屬和非金屬礦山［1］，特別是在鐵礦開采歷史悠久的鞍本地區表現得最為典型和突出。該地區既有礦業整合開發前偷采濫挖、采富棄貧和無序開采留下的采空區，又有日偽時期掠奪式開采留下的采空區，且它們數量龐大、規模不等和賦存空間不明。但從20 世紀末以來，隨著資源開發利用的科學技術水平不斷提高，主要是采掘運輸能力和選礦工藝水平的大幅提升，礦產資源整合露天開采成為一種趨勢，具有開采量大、礦石回收率高、生產安全、效益高、成本低等優點［2］。由于歷史原因和資源開采方式的轉變，露天礦山多由散亂小規模井下開采變為整合大規模露天開采，導致采場多處于不明采空區的發育區。而采空區頂板厚度會隨著持續的開采擾動不斷減小，直至突然塌陷形成深坑，輕則造成采掘設備墜毀，重則造成生產運輸車輛和人員損失。且隨著露天采礦活動的持續進行，采坑深度不斷加大，形成了高陡邊坡作用下的礦山開采應力場不穩多變，導致礦區采空區安全問題日益突出，不斷誘發坍塌、變形移位或含水空區突水等礦山次生地質災害，成為時刻威脅露天金屬礦山生產安全的主要難題［3-6］。

為解決上述遺留采空區誘發的生產安全技術難題，近年來，科研院所和生產單位采用了多種技術方法和手段對采空區進行探測，如人工調查法、地球物理探測、工程鉆探、放射性氡氣測量、三維激光掃描和水下聲吶測量等［7］，尤以地球物理探測方法最受青睞。該類方法具有種類多、成本低、方便快捷高效、無損非接觸的地質CT 探測成像等特點，屬于超前的“綠色探測”范疇，適用于對面積性采場區域內的采空區異常進行快速掃面篩查和定位，確保后期高效的鉆探驗證和三維激光或水下聲吶精準測量工作有效開展。但地球物理探測成果影像顯示異常具有間接性、多解性和圈定指標的不確定性，這使得探測成果的解譯和應用需要有豐富經驗的地質與地球物理專家級的人才聯合參與。現階段，在露天礦山快速持續不斷的平臺作業生產中，每天采集的大量多方法物探影像數據難以得到具有豐富經驗的專家級人才及時處理，致使采礦作業平臺采空區異常得不到及時、準確圈定，造成鉆探驗證和采空區安全處理工作嚴重滯后，在制約礦山正常生產的同時威脅了礦山安全。針對該類應用地球科學問題僅依靠專家經驗解決的局限性，隨著計算機科學技術和數學地質建模技術的發展，大量學者提出了多種基于專家經驗和統計分析驅動的數據分析處理判別模型，實現了大部分地球科學數據成果的定性和半定量化分析研究與應用普及［8-9］。特別是模糊數學理論的出現，使得地球物理異常數據的處理和定量解譯成為可能［10］。采空區地球物理異常作為礦致地球物理異常之上的二次人造異常，更具多樣性、變異性和復雜性，造成的不確定性具有模糊性，尤其是采空區地球物理異常的最終解譯成果只能反映異常圈定區內有無采空區存在，而這類“有、“無”之間就是一種模糊關系。

為解決采空區地球物理異常圈定的不確定性和專家經驗圈定的滯后性，便于及時開展鉆探驗證和安全處理工作，基于多年的BIF 型鐵礦山采空區探測成果，特別是采空區類型和地球物理特征研究成果的總結，應用模糊數學理論和模型評價方法，基于GIS 平臺構建采空區空間地球物理異常影像快速處理與定量評價模型，以期實現金屬礦山不明采空區地球物理異常圈定的模型化和智能化。

1 露天鐵礦采空區類型和地球物理識別特征

1.1 露天鐵礦采空區類型

采空區作為一種地表以下一定深度礦產資源開采后形成的空間，目前主要存在如下類型劃分思路［1］。

（1）基于資源開采時間的采空區分類。根據采空區的形成時間，采空區可以分為3類：第1類為歷史時期形成的老采空區，第2類為目前正在開采形成的現采空區，第3類為未來計劃開采區域形成的采空區。

（2）基于開采礦產資源類型的采空區分類。根據開采礦產資源類型的不同，采空區可以分為金屬礦山采空區和非金屬礦山采空區兩大類，其中金屬礦山采空區依據采礦方法又可分為空場法采空區、崩落法采空區和充填法采空區。

（3）基于采空區規模大小分類。可以簡單分為獨立礦山采空區和整體礦山采空區兩類。其中，獨立礦山采空區按體積可以分為0.1～1 萬m3、1～3 萬m3和3～10萬m33類，整體礦山采空區按體積規模可以分為5類，分別為＜50萬m3、50～100萬m3、100～500萬m3、500～1 000 萬m3和＞1 000 萬m3。此外，還可以劃分為小規模的巷道式采空區和大規模的采場式采空區。

（4）基于采空區埋藏深度分類。根據采空區埋藏深度，可以將采空區分為兩類：第1 類為淺部采空區，埋藏深度在60 m 以內，是礦山開采和工程施工必須及時探明和安全處理的采空區；第2類為深埋采空區，埋藏深度大于60 m，需要進行安全評價、監測與處理的采空區。

（5）基于采空區賦存狀態分類。根據采空區的賦存狀態，既可以分為已知采空區和未知采空區，也可以分為已塌陷采空區和未塌陷采空區。未知采空區，即為不明采空區或隱伏采空區，是地球物理探測和礦山安全防治的重點對象。

（6）基于采空區充填介質類型的物性分類。不論何種采空區，一旦形成，在未被人為充填和完全塌陷之前，總保持一定的形態。對于其空間內的充填介質，除了常見的圍巖和礦石碎塊外，對其物性影響最大的介質即為四處流動的空氣和水。依據對采空區物性差異的影響程度，一般可以將不同充填介質類型的采空區劃分為兩大類：第1類為空氣充填型采空區，其內充填介質以空氣為主；第2類為水充填型采空區，其內充填介質以水為主。該分類方式也是不明采空區可以進行地球物理探測的前提條件和理論基礎。

1.2 采空區地球物理識別特征

不同類型的采空區具有不同的地球物理特征，而不同地球物理特征的采空區來自于不同物性差異對采空區類型及賦存圍巖類型的可區分程度，這種區分程度也是選擇不同地球物理探測方法的理論基礎。基于物性差異的認識，本研究將采空區分為兩大類，即空氣充填型采空區和水充填型采空區。

（1）空氣充填型采空區。空氣充填型采空區以一種極低密度絕緣體空氣為充填介質。該類采空區與賦存空間的巖礦石相比，具有極低密度、高電阻率和低電導率的地球物理識別特征，可以采用高密度電阻率法、瞬變電磁法和地震映像法進行組合探測。相關研究表明，空氣充填型采空區異常在高密度電阻率法和瞬變電磁法數據影像上均表現為高阻值異常圈閉空間，且阻值從異常邊界到異常核心逐漸增大，呈正向遞增梯度變化趨勢［11-12］；在對地震映像法探測數據進行解譯時，可通過抽取最佳偏移距道數據來快速定位采空區異常位置，具體表現為地震波形圖同相軸的突然增減和錯斷，且主要以同相軸的錯斷增多為采空區異常存在特征［13-14］。

（2）水充填型采空區。水充填型采空區以一種來源于礦山水循環系統、富含多種礦物質元素的水為充填介質，對外顯示為低密度的良導體。該類采空區與賦存空間的巖礦石相比，具有低密度、低電阻率和高電導率的地球物理識別特征，也可以選用分辨率高的高密度電阻率法、瞬變電磁法和地震映像法進行組合探測。相關研究表明，水充填型采空區異常在高密度電阻率法和瞬變電磁法數據影像上均表現為低阻值異常圈閉空間，且阻值從異常邊界到異常核心逐漸減小，呈負向遞減梯度變化趨勢［7，11］；對于地震映像法探測數據，水充填型采空區與空氣充填型采空區的地球物理異常識別特征一致［13-14］。

2 模糊數學理論評價方法

模糊數學評價方法作為一種綜合評價方法，主要對受到多個不確定性因素影響的目標進行整體性評價［15-16］，這就需要首先確定評價對象的影響因素。該類因素對評價結果的影響程度不同，而需要確定該類因素對上層評價目標的權重，并經過模糊變換得到該類因素對評價總目標的隸屬度，最終完成對研究目標的綜合評價［17］。假設2 個有限集W=，且W 與K 的模糊關系可表示為K=WY。其中，Y 是W 到K 的一個模糊變換；W 為影響因素權重集；K 為模糊綜合評價成果集，表示最終綜合評價結果；Y 為模糊判別矩陣，描述因素權重集W和評價成果集K之間的模糊關系。

本研究考慮到物探數據的異常不確定性和間接性，為最大限度減少人為干預對評價結果的影響，故采用全程模糊綜合評價模型。首先，依據前期研究成果選出采空區異常評價指標因素集，通過數據重構提取每個評價指標因素，形成評價指標數據影像集；然后，對于獲取的每個評價指標數據進行自然間斷分級和模糊賦權分類，形成評價指標概率分布影像集；最后，通過模糊綜合評價模型判別采空區綜合地球物理異常，形成采空區異常分布概率影像圖。

3 露天鐵礦采空區異常評價指標

根據大量的BIF 型鐵礦采空區地球物理探測研究和鉆探驗證成果經驗總結［7，11-12，18］，優化組合出集高密度電阻率法、瞬變電磁法和地震映像法于一體的采空區異常探測方法，并利用方法中采空區與周圍巖礦石物性差異最大的一組物性參數組合來評價采空區的存在與否，如電法中的電阻率和電磁法中的電導率等。此外，依據探測數據中的隱含數據信息和不明采空區異常的空間賦存規律，如空氣充填型采空區的高阻異常圈閉性和地表的高阻異常不可能為采空區異常等，建立了如圖1所示的采空區地球物理異常評價指標體系。對于具體評價指標的數據特征和數據重構及其綜合評價，本研究以空氣充填型采空區為例進行分析。

3.1 高密度電阻率法數據特征和數據重構

利用高密度電阻率法探測采空區，主要是利用采空區與周圍巖礦石的電阻率差異，具體表現在電阻率數據二維影像上顯示為相對高阻值塊區和相對低阻值塊區（圖2（a）和圖2（b））。對于空氣充填型采空區，高密度電阻率法數據特征主要表現為高阻值異常圈閉區，從異常區的邊界到核心，阻值逐漸增高，表現為正向梯度變化。考慮到空氣充填型采空區具有正向梯度變化特征，可以采用離散傅里葉變換垂向求導獲得正向梯度值圈閉區域［19］（圖2（c）），即評價指標體系中的采空區異常梯度矢量特征G2T。對于電阻率高低的評價指標G1R，本研究按照傳統地球科學異常評價方法，在無遺漏采空區異常標準指導下，以每幅影像的電阻率平均值為背景值，高于背景值的區塊為正異常→高阻值圈閉區→空氣充填型采空區異常，低于背景值的區塊為負異常→低阻值圈閉區→水充填型采空區異常。此外，探測區域的地表接地條件、不明采空區的賦存位置和探測方法數據的空間分布等均可對采空區異常圈定的可靠性產生極大影響，如地表的正向梯度異常和高阻值圈閉難以成為采空區異常。為此，本研究提出了采空區異常埋深條件G3H作為評價指標，評價電阻率法探測數據點的空間分布可信度。

上述3 個評價指標G1R、G2T和G3H，由于數據格式和數據取值范圍存在很大差異，需要進行標準化的數據重構處理，主要是對各個評標指標涉及的參數值區間進行分類賦值。為使各個類所涵蓋的數值區間差異最大化，并體現出固有的自然分組圈閉特征，本研究采用了自然間斷點10 級分類法，形成了高密度電阻率法的3 個采空區異常評價指標影像分級圖（圖2（d）、圖2（e）和圖2（f））。

3.2 瞬變電磁法數據特征和數據重構

瞬變電磁法通過測量地下介質由于一次發射脈沖電磁場產生的二次感應電磁場隨時間變化的衰減特征，來獲取地下不均勻體的導電性能和位置（圖3（a）），其反演結果經過處理后如圖3（b）所示，但由于瞬變電磁法存在關斷時間效應，其淺部為探測盲區。空氣充填型采空區在瞬變電磁法成果圖中也表現為高阻值異常圈閉特征，從圈閉異常的邊緣到核心，電阻率逐漸增高，同樣可以采用離散傅里葉變換垂向求導獲得正向梯度值圈閉區域（圖3（c）），即評價指標體系中的采空區異常梯度矢量特征S2T。對于電阻率評價指標S1R，可采用探測成果電阻率的平均值為背景值，高于背景值的區塊即可判定為空氣充填型采空區異常。由于瞬變電磁法存在探測盲區，且隨著探測信號的衰減造成深部成果的可信度較低，故本研究提出采空區異常埋深條件S3H來綜合評價探測數據點的空間異常可信度。

為使瞬變電磁法的探測成果可與高密度電阻率法的探測成果進行模糊綜合評價，需對其3個評價指標（S1R、S2T和S3H）進行標準化的數據重構處理。數據重構處理同樣采用自然間斷點10 級分類法，形成了瞬變電磁法的3 個評價指標影像分級圖（圖3（d）、圖3（e）和圖3（f））。

3.3 地震映像法數據特征和數據重構

基于最佳偏移距的地震映像法主要用于快速確定采空區異常存在與否和空間位置，其埋深和規模可通過與高密度電阻率法和瞬變電磁法綜合評價獲得。采空區在地震映像法中主要表現為同相軸的錯斷和增多（圖4（a）），個別顯示為同相軸的錯斷和減少。為此，本研究將空氣充填型采空區在地震映像法中的數據異常特征分為3 個類別進行位置區間賦值重構，無采空區異常位置區段賦值為1，采空區異常顯示一般區段賦值為3，采空區異常顯示強區段賦值為5，得到地震映像法的采空區異常評價指標D1W（圖4（b））。

4 基于模糊數學模型的露天鐵礦采空區探測地球物理異常評價

依據高密度電阻率法、瞬變電磁法和地震映像法的探測成果，提取了7個采空區地球物理異常評價指標，分別為G1R、G2T、G3H、S1R、S2T、S3H和D1W，且當它們每個指標在一個指定區段內都達到權重概率預定值以上，方可綜合判定為采空區異常。同時，需在保證無遺漏采空區的條件下縮小探測數據異常圈定范圍，以降低鉆探驗證成本。因此，有必要對7 個不確定性評價指標合理提取關鍵參數并進行賦權，確保評價模型的可靠性和拓展性，從而可適用于復雜多變的礦山地球物理探測評價，達到采空區地球物理異常評價的定量化和智能化。本研究采用已驗證區域的空氣充填型采空區的震—電—磁探測數據來開展研究（圖5），同時，為了確保一致性和系統性，文中用到的震—電—磁地球物理探測成果數據為同一組數據。

4.1 評價指標模糊分類

7 個采空區異常評價指標可以分為4 類：第1 類為探測數據值指標，為G1R和S1R；第2類為探測數據提取指標，為G2T和S2T；第3 類為有效探測區間指標，為G3H和S3H；第4類為采空區探測異常位置指標，為D1W。對于上述每類評價指標的模糊分類，可使已進行數據重構的評價指標歸一化和概率化，即將采空區異常評價指標分類數據轉換到0～1 的概率區間內。該區間內，“0”分配給確定不屬于指定采空區異常集合的分類值，“1”分配給確定屬于指定采空區異常集合的分類值，0～1 區間內的概率值分配給可能為采空區異常的分類值。

4.1.1 第一類采空區異常評價指標模糊分類

G1R代表高密度電阻率法的直接電阻率測量值，對于空氣充填型采空區，在無遺漏采空區異常標準指導下，探測的電阻率值若大于背景值或探測影像的電阻率平均值，即可認定為采空區異常，可以采用模糊較大值計算模型進行分析：

式中，x為輸入計算值；f1為定義過渡區形狀和特征的展開參數，本研究取3；f2為定義的起算中點，用于確定交叉點（起算中點分配概率值為0.5，大于起算中點的值為采空區異常的概率較高，小于起算中點的值為采空區異常的概率較低）。

該模型以評價指標G1R中探測影像電阻率平均值所代表的分類值（圖2（e））為起算中點f2，其起算概率為0.5，將柵格數據G1R（圖2（e））圖像值作為輸入值進行模糊計算，得到采空區異常評價指標G1R指示采空區異常存在的概率分布圖（圖6（a）），其概率值越大，屬于采空區的可能性便越大。同理，可以采用模糊較大值計算模型對瞬變電磁法采空區異常評價指標S1R（圖3（e））進行采空區異常概率分布計算，其起算中點f2為探測影像電阻率均值所代表的分類值，可以得到采空區異常評價指標S1R指示采空區異常存在的概率分布圖（圖6（b））。

4.1.2 第2類采空區異常評價指標模糊分類

第2 類采空區異常評價指標是高密度電阻率法的G2T和瞬變電磁法的S2T。G2T為電阻率的空間變化梯度值，正梯度值為高阻值異常圈閉區顯示，代表空氣充填型采空區異常。該評價指標可以采用模糊高斯模型進行計算：

式中，x為輸入計算值；f1為定義過渡區形狀和特征的展開參數，本研究取3；f2為定義的起算中點，用于確定交叉點（起算中點分配概率值為1，大于起算中點的值或小于起算中點的值對應為采空區異常的概率均較低）。

模糊高斯模型以指標G2T平均值所代表的分類值為起算中點f2，起算概率為1，小于或大于起算中點的概率均較低，該模型可以有效消除地表接觸不良引起的正梯度值超大“采空區異常”，也可以規避將正梯度值過低的采空區異常邊界過渡區確定為采空區異常。本研究將柵格數據G2T（圖2（d））圖像值作為輸入計算值進行模糊計算，得到了采空區異常評價指標G2T指示采空區異常存在的概率分布圖（圖6（c））。

對于瞬變電磁法的采空區異常評價指標S2T，由于存在關斷時間導致的淺部探測盲區，可有效消除因地表接觸不良造成的“高阻采空區異常圈閉”，其正梯度值（圖3（d））變化越大，顯示為采空區異常的概率便越大，也可采用模糊較大值計算模型進行計算，其起算中點f2為正梯度值平均值所代表的分類值，計算得到的采空區異常評價指標S2T指示采空區異常存在的概率分布如圖6（d）所示。

4.1.3 第3類采空區異常評價指標模糊分類

第3 類采空區異常評價指標為有效探測區間指標，為G3H和S3H（圖2（f）和圖3（f））。該類評價指標可以有效篩除淺部接觸不良導致的異常數據、可信度低的深邊部探測數據和關斷時間導致的探測盲區數據，以降低采空區異常誤判概率。該類指標可采用模糊高斯模型（式（2））進行計算，并依據前期探測和驗證成果，將有效探測數據深度區間的中心深度值所代表的分類值作為起算中點f2，賦權概率為1，而淺部或深部異常數據為采空區的概率均較低，得到G3H和S3H的采空區異常概率分布圖（圖6（e）和6（f））。

4.1.4 第4類采空區異常評價指標模糊分類

地震映像法的采空區異常位置指標D1W，可依據其異常分類特征采用模糊較大值模型（式（1））計算，f1取值為2，其起算中點f2為3 級異常分類中的1 級和2 級采空區異常的中間值2（圖4（b）），起算概率為0.5，圖像計算得到的采空區異常評價指標D1W的采空區異常概率分布圖如圖6（g）所示。

4.2 評價模型構建

本研究采用數據重構方法和模糊數學計算模型提取了7 個采空區異常評價指標G1R、G2T、G3H、S1R、S2T、S3H和D1W，綜合評價采空區異常存在與否。鑒于上述每個采空區異常評價指標均來自不確定性判別數據集或依據專家經驗總結得出，且評價指標模糊分類中的概率賦值也選用了最低可能性，以免遺漏采空區異常，具體表現在有效探測數據區間內，空區存在的區域必有采空區異常，而有采空區異常的區域未必全是空區。因此，考慮到指標體系的不確定性，本研究空區異常評價采用模糊綜合評價模型（圖6（h））［20］，即模糊疊加模型，具體可為模糊And 函數評價模型：

式中，arg1,arg2,arg3,…,argn 為不同的評價指標取值，本研究n ∈[ ]1,7，Fuzzy And 疊加類型將返回單元位置所屬集合的最小值。

模糊And疊加模型會創建一個輸出柵格，為像元位置所屬的每個集合的每個像元值給定了最小模糊值。本研究中，采空區異常評價模型是7個指標條件中的每個條件都已相對于其屬于合適集合的隸屬度被模糊化，模糊And模型會識別多個條件內像元屬于一個合適集合的最小可能性。對于采空區地球物理異常評價，主要依據無遺漏采空區異常標準和模糊分類的概率賦值原則，只需要選擇符合所有指標中概率值大于或等于0.5 的具有采空區異常的空間位置集合，疊加生成采空區異常指標綜合評價概率分布圖，圖中概率值大于0.5 的區塊圈閉即為采空區異常，需要進一步進行鉆探驗證。

4.3 評價結果

基于上述綜合評價結果，采用ArcGIS 軟件中的空間分析模塊進行弓長嶺鐵礦區多采空區異常評價指標圖像矩陣運算，計算模型為Fuzzy And，結果見圖7。圖7中測線起點以高密度電阻率法G8測線起點為0 點，與其平行的瞬變電磁法和地震映像法測線則以相對于電法起點疊加距離作為起點計算空間坐標。在圖7 中，根據0.5 及其以上概率值作為異常圈定標準，圈定出了一個采空區異常區（即空氣充填型采空區）。在對異常區進行驗證時，可將采空區異常區上部邊界與驗證得出的采空區頂板埋深進行對比分析。圖7 中采空區異常區圈定的頂板埋深為35.7 m，經過鉆探驗證埋深為33.8 m，兩者相差1.9 m，相比于傳統經驗法的圈定誤差（±5.0 m 以內），精度提高了3.0 m以上。

通過對空氣充填型采空區的地球物理探測數據進行指標化提取，并采用本研究提出的綜合評價模型，在ArcGIS 平臺上定量圈定了采空區異常邊界，通過鉆探成果相比，可以在無遺漏采空區的前提下，大大提高采空區圈定效率和準確度。在上述評價方法指導下，對弓長嶺鐵礦區水充填型采空區進行了綜合評價，圈定的采空區異常頂板埋深18.3 m，鉆探驗證結果為20.1 m，兩者相差在2 m 以內。可見，本研究選擇的采空區異常圈定指標具有合理性，即對于不確定性的地球物理異常可以采用模糊數學方法進行定量研究，實現地球物理探測數據評價的指標化、定量化和智能化。此外，本研究評價方法可以通過ArcGIS 平臺對采集到的大量探測數據進行流程式的智能化快速綜合評價，并可依據評價結果進行快速鉆探驗證，確保礦山生產安全。

5 結 論

（1）劃分了多種采空區類型，依據采空區充填介質和物性差異的不同，將采空區分為空氣充填型采空區和水充填型采空區兩大類，并分別建立了相應的地球物理識別特征和震—電—磁探測方法組合。

（2）根據高密度電阻率法、瞬變電磁法和地震映像法的采空區異常數據識別特征，構建并提取了7個采空區異常評價指標G1R、G2T、G3H、S1R、S2T、S3H和D1W。

（3）通過對7個采空區異常指標進行模糊分類和模糊數學模型綜合評價能夠定量圈定出采空區異常邊界，可使探測精度保持在±2.0 m以內，為采空區異常快速定量化和智能化評價提供了技術手段。