基于多種數據融合方法的采空區富水性研究

周子鵬，王施智，沈福斌

(1．國土資源部煤炭資源勘查與綜合利用重點實驗室，陜西 西安 710021；2.陜西省煤田物探測繪有限公司，陜西 西安 710005)

0 引言

2020年中國發生各類煤礦死亡事故122起，死亡人數225人，其中重特大事故3起，包括湖南省衡陽市耒陽市導子煤業有限公司源江山煤礦“11·29”重大透水事故，造成13人死亡。盡管防治水工作已成為重中之重，但單一勘探方法不可避免地具有片面性及多解性，因此將各種方法聯合或許會取得較為滿意的結果。張長明等[1]將井下瞬變電磁技術和礦井音頻電透視技術結合使用探查煤層底板富水情況，發揮了各自的方法優勢，提高了勘探精度；牟義等[2]利用槽波地震、無線電波透視、瞬變電磁法、音頻電透視探查了晉煤礦區采煤工作面的災害情況，包括斷層、陷落柱、巖層富水情況，解釋成果準確性顯著提升；江微娜等[3]采用地面瞬變電磁與井下音頻電透視聯合探測了晉城老礦區的采空區情況，取得了一定的成果；高波、范濤等[45]聯合了地面瞬變電磁及高密度電法，較好的劃分出了巖層富水區，提高了劃分的精確度。楊文采院士[6]曾指出，只有聯合反演才是地球物理解釋的出路。基于以上研究，使用井下瞬變電磁法、直流電法以及音頻電透視法進行煤礦采空區的探查，后期根據其原理之間的聯系和差異，將不同方法的參數進行轉化，利用數據融合方法最終劃分了勘探區內采空區的積水范圍。

1 數據融合

作為一種數據處理方法，數據融合的中心思想是從多源、多類別的數據之中最大程度地提取有用信息[7]。

1.1 數據標準化

作為融合數據的第1步，須將各數據作標準化變化進行規整，使得數據更加穩定。

假設獲得電阻率數據為ρ=(ρ1,ρ2,ρ3,…,ρn)，則標準化公式為

(1)

1.2 數據歸一化

歸一化的目的是為了讓數據統一歸屬于(0，1)范圍內，同時也將數據進行了尺度、量綱統一，便于不同單位的數據能夠進行比較和運算。

線性歸一化公式為

(2)

式中，ρ歸為電阻率歸一化結果；ρi為視電阻率，i=1，2，3,…,n；ρmin為視電阻率最小值；ρmax為視電阻率最大值。

1.3 相關性分析

按照上述方法將不同采集數據進行變換后，將數據與原始建模的、包含采空區數據的模型進行對比，分析其相關性，獲得關聯度。當數據較為冗余龐大且關聯性較高時，可以采用主成分分析法獲取集中的有效信息變量[8]。

簡單相關系數

基性含長結構帶，其中可見黑云母化或金云母化，部分見于角閃石角巖內，其中SiO2為帶入組分，而Ai、Ti、Fe、Ca、Mn、Mg、P、K、Na、Cr、V為帶出組分。其厚度范圍較廣，可達800～1 000 m。

(3)

式中，Cov(x,y)為x與y的協方差；Var[x]為x的方差；Var[y]為y的方差。

1.4 主成分分析(PCA)

主成分分析的內容主要是將原本有關聯的幾組或多組信息經過變化重新組合成新的無關聯的綜合指標，這些新的指標成為主成分[89]。

在整個線性組合中，Var[f1] 越大，表示f1包含的信息越多，故選取的f1應該為maxVar[f]，則稱f1為第一主成分。若第一主成分不足以代表原來n個指標的全部信息，即考慮選取f2，為了簡潔、有效地反映原始信息，要求Cov(f1，f2)=0，即重合的部分不需要再次使用，稱f2為第二主成分，依此類推可以構造出第三，第四，……，第n個主成分。

主要的分析步驟如下。

計算相關系數矩陣：主要目的是確定相關系數大小。

(4)

式中，rij為原始數據變量之間的相關系數，即(3)式的各元素計算結果；R為實對稱矩陣。

求取特征值和特征向量：根據特征方程|λE-R|=0，利用Jacobi法求出特征值，再根據特征值λi求出對應的特征向量ei(i=1，2，3,…,n)。

2 工程布置

本次選取某處具有采空區的煤礦113202工作面進行實踐。分別采用瞬變電磁法、直流電法、音頻電透視法進行探測，工作面及測點相對位置如圖1所示。

圖1 工作面及測點相對位置示意Fig.1 Schematic diagram of the relative positions of the working face and measuring points

儀器及工程布置如下：

瞬變電磁探測使用YCS512本安型瞬變電磁儀，中心回線裝置，發射線圈為1.5 m的正方形線圈，接收線圈為0.8 m的正方形線圈，發射和接收線圈分別為10匝和20匝，發射電流3 A，點距10 m。共設計物理點808個。

直流電法探測采用YD32(A)高分辨電法儀，根據工作區地質資料，3-2號煤層與2-2號煤層間距約42 m，確定供電極距OA采用算術極距，極間距為10 m，最小OA選為15 m，最大OA選為75 m，無窮遠極OB不小于300 m。接收極距MN/2為5 m。設計物理點52個。

采用YT120(A)音頻電穿透儀，工作面寬320 m，為確保信噪比，發射、接收頻率選擇為15 Hz。發射點距為50 m，接收點距為10 m。設計物理點52個。

3 數據提取及變換

3.1 原始數據分析

瞬變電磁法、直流電法及音頻電透視法實測數據空間分布示意圖如圖2所示，由實測數據分布圖可以大致判斷突變點的位置、干擾位置及數據質量和差異，為下一步剖面解釋打好基礎工作。如圖2(a)所示，瞬變電磁數據質量較好，存在少數尖點，推測為偶然干擾，在大號點方向視電阻率值較大，且由小到大發生突變；從圖2(b)可以看出，直流電法受干擾較小，整體較為平緩，視電阻率在兩端遠距離處較大；圖2(c)顯示音頻電透視數據在近距離處變化較大，排除干擾的可能性后可根據其變化尋找差異，圈定異常區域。

圖2 3種方法實測數據空間分布示意Fig.2 Schematic diagram of the spatial distribution of measured data for the three methods

剖面圖可以反應數據橫向、縱向的整體情況，測區內瞬變電磁典型剖面如圖3所示，以此來判斷這一測線上數據的完整性和準確性，同時確定方法的有效性。

圖3 瞬變電磁典型視電阻率等值線斷面(探測方向與頂板呈90°)Fig.3 Transient electromagnetic typical apparent resistivity contour section (detection direction is 90° to the roof)

圖3中黑色虛線為頂板2-2煤層縱向位置。如圖，在0～80 m范圍內，視電阻率等值線呈低阻圈閉，根據礦井回采資料，該范圍內頂板2-2煤層沒有采空，因此分析該低阻異常為頂板砂巖含水的反映；在350～600 m和700～850 m之間，視電阻率值相對較低，等值線呈低阻圈閉，該范圍內頂板為2-2煤采空區，故分析此2處低阻異常為采空區積水的反應。確定瞬變電磁法圈定的異常范圍確定為“高敏感區”。

3.2 插值方法選擇

由于3種方法的探測方式不同，所得到的數據格式、數量也不同，因此在實際操作過程中，要選擇一種插值方法來對數據進行插值獲得擁有相同數據坐標的數據點。

在眾多的插值方法中，由于克里金插值算法能夠描述區域化變量的空間結構性變化，又能描述其隨機性變化，其變量之間都具有空間上的相關性，且隨機誤差也具備二階平穩性，因此成為地學上空間內插值的不二之選[12]。其計算公式為

(5)

式中，x1,…,xn為系列空間上的觀測點；z(x1),…,z(xn)為相應觀測值；z*(x0)為區域變化量在xi處的值；δi為變量，其選取符合無偏性和估計方差最小的原則。

3.3 數據主成分分析

首先對插值得到的數據進行標準化、歸一化，由于進行的是線性變換，因此在變換前后其基本特征不變，即數據失真的程度較小[9，11]。

在對數據進行標準化和歸一化之后，將3種數據共引入進行主成分分析，根據前述以瞬變電磁數據為主，減少數據的冗余度后[16]，討論音頻和直流2種方法與其相關性，得到擬視電阻率特征見表1，原始數據與主成分相關系數見表2，取特征值大于1的成分作為主成分，小于1的舍去，即最后獲得主成分為F1、F2。

表1 擬視電阻率數值特征一覽

由表2可以看出，主成分之間相關度不高，因此可以進行獨立分析。F1與瞬變電磁視電阻率值相關系數最高，F2與音頻電透視視電導率相關系數最高，根據前文主要使用F1進行地質分析。

表2 原始數據與主成分相關性系數

根據實際井下情況，F1反應的地質區域與采空區富水可能性高的區域聯系緊密，后面就F1的幅值特征來對區內的采空區特性進行解釋。

3.4 平面解釋

最終解釋成果圖如圖4所示。大于0.5值初步判斷為高阻異常，小于-0.5值判斷為低阻異常，以此為依據，以視電阻率-0.9共圈定富水異常區3處，Ⅰ號異常區范圍內頂板滴水較大，且煤層未采空，分析為頂板砂巖富水異常；Ⅱ號異常區范圍內雖然視電阻率值較低，但煤層未采空，故分析為頂板砂巖富水異常；Ⅲ號異常區靠近工作面輔運順槽一側，此處異常范圍較大，根據地層巖性分析后，判斷其為采空積水區，且積水量較大。

圖4 解釋成果Fig.4 Interpretation results

經疏放水工作驗證，本次勘探圈定的異常區較為準確，故融合算法也可以作為常用的反演算法。

4 結論

(1)瞬變電磁法對低阻含水體較為敏感，能夠探查大范圍的富、含水異常，可作為掃面性工作的主要方法。

(2)直流電法可作為瞬變電磁法的輔助方法，其勘探結果與瞬變電磁法結果的關聯性較大。

(3)主成分分析法的結果需要結合區域地質特征進行異常區的劃分。

(4)數據融合計算后所得解釋成果圖，較單一原始成果更為準確，且計算過程快速高效，可以應用于實際工作中。