張 碩，趙文龍

(中鐵二院工程集團有限責任公司，成都 610031)

0 引言

我國工業(yè)和國民經(jīng)濟的持續(xù)高速發(fā)展，都離不開能源產(chǎn)業(yè)和礦產(chǎn)資源強有力的支撐。而煤炭資源的開采已是國家獲取能源礦產(chǎn)的主要途徑，隨著開采速度的加快，其開采深度和開采規(guī)模也不斷擴大，礦井突水、淹井事故已經(jīng)成為煤礦的重大災(zāi)害之一，給各大井田的安全生產(chǎn)帶來威脅。據(jù)統(tǒng)計國內(nèi)80%的礦井突水事故是由于構(gòu)造引起的[1-2]。陷落柱一旦成為地下水的通道，將會對礦井的安全生產(chǎn)造成巨大威脅，特別是陷落柱可以導通多層含水層，使整體的水力聯(lián)系增強，可將煤層底板下承壓的奧灰水導入礦井，造成重大災(zāi)害[3-4]。

近年來，直流電法技術(shù)蓬勃發(fā)展，在探測巷道頂?shù)装甯羲畬雍穸?，含水、導水?gòu)造，斷裂破碎帶以及工作面內(nèi)部小構(gòu)造等都取得了顯著的效果，但在井下探測采煤工作面底板圍巖中隱伏的含水，導水構(gòu)造方面仍是一個難題[5]。由于直流電的穿透特性及高阻煤層和巷道空間對電流的屏蔽，礦井直流電透視技術(shù)對于煤層底板內(nèi)的隱伏導水構(gòu)造具有良好的響應(yīng)規(guī)律。目前，已有幾種典型的建場方法，二維和三維條件的正演計算也頗有成果，資料處理的方法主要為曲線對比法和地電成像法，相應(yīng)的反演成像研究正在開展[6-7]，但相對還很欠缺。隨著計算機硬件的提升，直流電阻率法的三維反演技術(shù)不斷發(fā)展，在起伏地形條件，約束性反演等方面取得了較大進步[8-10]。為了進一步豐富直流電透視資料解釋方法，筆者以直流電法三維正反演程序與直流電透視方法相結(jié)合，采用巷道間AM法，直接測量電位，對隱伏含水陷落柱的探測進行了詳細研究。

1 正演邊界條件及數(shù)值模擬

1.1 正演邊界條件

在穩(wěn)定電流場中，由歐姆定律的微分形式出發(fā)，根據(jù)電荷守恒定理，高斯定理等可以得出點電源三維電場分布方程：

▽·(σ(x,y,z)▽u(x,y,z))=f

(1)

即：

(2)

因此正演要解決的問題，即是求解式(2)在滿足一定邊界條件下的解，筆者采用混合邊界條件，Dey等提出了混合邊界條件[14]：

(3)

其中：n為邊界?！尥夥ň€方向上的坐標變量；r為電流源到邊界上的距離；θ為n與r的夾角。該方法可以保持電位在邊界上的物理特性，可以減少邊界附近的網(wǎng)格數(shù)量，減少計算量。

1.2 模型設(shè)置

為了研究礦井直流電透視法對煤層底板陷落柱的響應(yīng)特征，采用等比原則，設(shè)置如圖1所示地電模型，采用有限差分法、采用A-M觀測方式進行正演計算：如圖1(a)所示為數(shù)值模型的三維結(jié)構(gòu)示意圖，從上到下分別為頂板圍巖、煤層、底板圍巖，圖中藍色方形體為擬設(shè)定的模擬導水陷落柱。圖1(b)為煤層底板的水平切片，圖1(c)為圖1(a)沿紅色線框切下的斷面圖。整個計算區(qū)域為190 m×130 m×124 m的長方體區(qū)域(網(wǎng)格為38×26×25)，整個區(qū)域存在兩種網(wǎng)格剖分形式，在煤層內(nèi)網(wǎng)格剖分為5 m×5 m×4 m，煤層厚度為4 m；其他區(qū)域網(wǎng)格剖分大小為5 m×5 m×5 m。煤層電阻率值設(shè)為1 000 Ω·m，頂?shù)装咫娮杪手翟O(shè)為100 Ω·m ，巷道的電阻率設(shè)為3 000 Ω·m ，回采工作面兩側(cè)巷道間距設(shè)為100 m。陷落柱的大小為20 m×20 m×30 m，電阻率值為10 Ω·m ，異常體模型的頂面位于煤層底板下深0 m處，在x方向上位于30個電極的中心，在y方向上位于兩巷道的中間位置。

圖1 數(shù)值模型參數(shù)設(shè)置說明圖Fig.1 Numerical model parameter setting instructions(a)裝置模型示意圖;(b)煤層底板切片;(c)YOZ斷面

模型中電極分別布設(shè)于兩側(cè)巷道的底板上，供電電極A布設(shè)30個，分別編號為A1、A2…、A30，點距為5 m，測量電極同樣布設(shè)30個，分別編號為M1、M2…、M30，測量點距為5 m。當A1點供電時，M1到M30號測點分別進行觀測，同樣當A2點供電時，M1到M30號測點分別進行觀測，這樣直到A30點供電完成，一次正演共測得900個數(shù)據(jù)。

1.3 數(shù)值模擬結(jié)果分析

選取A1、A10、A15供電時測量的數(shù)據(jù)繪制成電阻率曲線圖(圖2)。當A1電極供電時，從M1號電極到M30號電極所測得視電阻率值逐漸增大，最大值達到105 Ω·m 以上；當A10電極供電時，從M1號電極到M30號電極所測得視電阻率值逐漸增大，但整體值較A1號電極供電時要??；當A15號電極供電時，測得視電阻率曲線關(guān)于M15號測點對稱。整體計算視電阻率值均大于底板圍巖電阻率，即底板存在低阻異常體時二極裝置測得視電阻率值出現(xiàn)反異?，F(xiàn)象[15-16]。

圖2 A1、A5、A15供電視電阻率曲線圖Fig.2 Apparent resistivity curve of A1、A5、A15 Power Supply(a)觀測模型;(b)視電阻率曲線

2 反演成像理論

筆者采用擬高斯-牛頓法進行反演，運用共軛梯度算法求取模型修改量，提高計算速度[8]。同時，為了減少內(nèi)存的消耗，程序采用不直接求取雅克比矩陣J的計算方法[9]。

2.1 反演方程

設(shè)參數(shù)模型為m=(ρ1,ρ2,...,ρk)T，其中k為網(wǎng)格數(shù)量，觀測數(shù)據(jù)dobs=(ρs1,ρs2,...,ρsn)T，其中n為觀測數(shù)據(jù)量，模型正演計算數(shù)據(jù)為dcal。首先構(gòu)造反演目標函數(shù)：

(4)

引入拉格朗日常數(shù)λ和光滑度矩陣C，得到：

λ(CΔm)T(CΔm)

(5)

如果函數(shù)φ(m)的值趨于極小，則函數(shù)φ(m)的導數(shù)等于零：

(6)

因為dcal(mk+1)=dcal(mk)+J·(Δmk)，所以式(6)可寫為：

(JTJ+λCTC)Δm=JTΔd

(7)

當完成非線性問題的線性化，求解該方程組是三維電阻率反演成像中的重要的問題。解上述反演方程的困難在于三維問題的雅克比偏導數(shù)矩陣維數(shù)很大，無法直接求取，因此采用不直接求取J的方法。

2.2 雅克比矩陣J的處理

采用不直接求取J的方法，利用一次擬正演求取雅克比偏導數(shù)矩陣J與一個向量的乘積，從而大大降低了計算機的內(nèi)存占用，提高了計算速度[11-14]。

以一個3×3的矩陣進行運算，假設(shè)觀測數(shù)據(jù)個數(shù)為3，模型網(wǎng)格數(shù)量為3，以此為例，雅克比偏導數(shù)矩陣表達式為：

(8)

所以

(9)

在正演過程中由于

(10)

因此

(11)

(12)

(13)

所以

(14)

其中

(15)

所以兩邊求導得

(16)

由于A為對稱矩陣，所以

(17)

最終JTΔd可寫成如下形式：

(18)

這樣JTΔd的計算量就相當于以as1、as2、…、asK為場源的一次正演。按照這種方法，同樣可以利用一次擬正演計算出J·x的值，大大降低了反演過程中計算機內(nèi)存開銷，縮短計算時間。

2.3 共軛梯度算法

3 數(shù)值模型實驗

3.1 實驗?zāi)Ｐ?/h3>

建立三維空間的含水陷落柱模型，將觀測區(qū)域剖分為14 m×14 m×12 m的網(wǎng)格，設(shè)陷落柱模型規(guī)格為2 m×2 m×2 m，電阻率值為10 Ω·m ，底板圍巖電阻率為100 Ω·m ，覆蓋煤層電阻率為2 000 Ω·m 。陷落柱的中心位置坐標為(0 m，0 m，-1 m)，兩平行巷道間距離為9 m，左側(cè)為一連通巷道，空間位置如圖3所示，圖3(b)和圖3(c)分別為導水陷落柱在xoy平面和zoy平面上的投影。

圖3 含水陷落柱及電極布設(shè)模型Fig.3 The numerical model of collapsed column and electrode location(a)三維地電模型;(b)沿xoy切片；(c)沿zoy切片

3.2 反演結(jié)果分析

反演中取模型光滑因子為0.05，經(jīng)過五次迭代得到結(jié)果進行切片，如圖4所示。藍色低阻區(qū)域電阻率小于60 Ω·m ，呈橢球狀，反演結(jié)果可以很好地反映含水體的位置。

圖4 反演電阻率結(jié)果Fig.4 Inversion resistivity results

4 物理實驗及結(jié)果分析

對上述的數(shù)值模擬結(jié)果進行物理實驗的驗證，采用擬高斯-牛頓法對物理實驗數(shù)據(jù)反演。由于高阻煤層和巷道空間對電流的屏蔽作用，在巷道底板布設(shè)電極的電場分布特征與地面半空間中布設(shè)電極的電場分布特征極為一致，所以在物理模擬中采用了地面半空間模型，采用一定濃度的含食鹽泥漿模擬實際的含水陷落柱，如圖5所示,其中觀測電極為A1、A2、…、A30號，相鄰供電點距離為0.1 m，測量電極命名為M1、M2、…、M30號，相鄰測點距離為0.1 m，陷落柱大小為40 cm×40 cm×40 cm，電阻率值約5 Ω·m 。

圖5 物理模擬實驗設(shè)計與現(xiàn)場圖Fig.5 Design and site of physical simulation experiment(a)物理模擬實驗示意圖;(b)實驗現(xiàn)場

4.1 實測數(shù)據(jù)分析

在擬定異常區(qū)域開挖之前進行一次觀測，得出不存在導水陷落柱時的觀測值，然后根據(jù)實驗設(shè)計，按照設(shè)定的尺寸開挖、充水，再次進行測量，得到兩組視電阻率曲線的對比圖。

實驗數(shù)據(jù)提取M點測得的電位值，采用二極(AM)裝置系數(shù)計算電阻率值。為了較為直觀地對實驗數(shù)據(jù)進行分析，選取A15號供電電極供電時由30個測量電極測得的數(shù)據(jù)繪制電阻率曲線圖，如圖6所示，從M1到M30點觀測時，整體視電阻率表現(xiàn)出先增大后減小的趨勢，在M14號電極觀測的視電阻率值最高。物理實驗結(jié)果與數(shù)值模擬結(jié)果具有一致性。

圖6 A15號電極供電視電阻率曲線圖Fig.6 Apparent resistivity curve of A15 power supply(a)A15號電極供電觀測;(b)視電阻率曲線

結(jié)合數(shù)值模擬結(jié)果分析，測區(qū)AM之間存在低阻異常時，中間電極A供電時所測得電阻率值反而增大，出現(xiàn)了反異?，F(xiàn)象。該現(xiàn)象是由于低阻體在地表，對電流產(chǎn)生吸引，地表AM間整體電阻下降，導致測點M處測得電位升高。沈平等[15]在井間電阻率成像方法的研究中也提出了這種反異常的現(xiàn)象，并對此進行了解釋。湯井田等[16]在研究不同裝置下點源球體的近似解與精確解對比時，對A-M法出現(xiàn)正反異常的條件做了詳細的論述，認為極距AM應(yīng)小于異常體的埋深與半徑之比，否則會出現(xiàn)反異常的現(xiàn)象。

4.2 實測數(shù)據(jù)反演

我們用以上理論對上述物理實驗數(shù)據(jù)進行反演成像處理，反演初始模型電阻率值設(shè)為22 Ω·m ，模型網(wǎng)格為35×35×30，經(jīng)過5次迭代用時42 s，擬合誤差下降到5%。同一模型參數(shù)，采用直接求取雅克比偏導數(shù)矩陣J的方式進行計算迭代一次的時間為765 s。利用一次擬正演計算雅克比矩陣和一個向量乘積的方法，可以大大縮短計算時間。

物理模型實驗反演結(jié)果如圖7所示，實際模擬陷落柱在x方向和y方向坐標為(0 m，0 m)，中心深度為0.2 m。由反演結(jié)果可以看出，在探測區(qū)域的中心位置(實際模擬含水陷落柱所在位置)，出現(xiàn)了一個明顯的低阻異常區(qū)域。通過z方向上切片可以看出，低阻異常區(qū)的中心深度為0.6 m，模擬低阻異常區(qū)域?qū)嶋H挖掘的深度為0.2 m，考慮到在實驗進行中泥漿水不斷下滲，所以實際低阻異常的中心位置下移。該區(qū)域在水平位置垂向位置上都和模擬低阻體所在位置具有一致性，反演結(jié)果能較好地反映模擬低阻異常體所在三維空間中的位置。從反演結(jié)果中可以看出存在一處明顯的高阻異常區(qū)域，由實際實驗區(qū)域電阻率分布不均造成。

圖7 實測數(shù)據(jù)反演結(jié)果切片圖Fig.7 The map of 3D inversion used measured data

5 結(jié)論與建議

通過對陷落柱模型的數(shù)值模擬和物理模擬，以及實驗結(jié)果的三維反演成像研究得出以下結(jié)論：

1)通過對導水陷落柱模型的正演計算結(jié)果可以看出，采用常規(guī)視電阻率計算公式進行計算時，在低阻陷落柱區(qū)域計算結(jié)果為高阻，出現(xiàn)反異?，F(xiàn)象，測點距離陷落柱越遠，視電阻率越接近圍巖電阻率值。

2)在進行三維反演時，采用不直接求取雅克比矩陣J的方法，利用一次擬正演求取雅克比矩陣J與一個向量的乘積，從而大大降低了計算機的內(nèi)存占用，加快了計算速度。

3)對物理實驗數(shù)據(jù)進行三維反演，結(jié)果很好地反映出了導水陷落柱的實際位置，數(shù)值模擬的視電阻率計算結(jié)果和實測數(shù)據(jù)的反演結(jié)果，為實際電透視探測陷落柱的解釋提供了借鑒。

4)以后將繼續(xù)針對不同的觀測系統(tǒng)，不同空間位置的異常體進行研究，進一步總結(jié)探測規(guī)律。