直流電透視三維反演技術(shù)在探測(cè)含水陷落柱中的應(yīng)用研究

張 碩，趙文龍

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

0 引言

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

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

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

1.1 正演邊界條件

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

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

(1)

即：

(2)

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

(3)

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

1.2 模型設(shè)置

為了研究礦井直流電透視法對(duì)煤層底板陷落柱的響應(yīng)特征，采用等比原則，設(shè)置如圖1所示地電模型，采用有限差分法、采用A-M觀測(cè)方式進(jìn)行正演計(jì)算：如圖1(a)所示為數(shù)值模型的三維結(jié)構(gòu)示意圖，從上到下分別為頂板圍巖、煤層、底板圍巖，圖中藍(lán)色方形體為擬設(shè)定的模擬導(dǎo)水陷落柱。圖1(b)為煤層底板的水平切片，圖1(c)為圖1(a)沿紅色線框切下的斷面圖。整個(gè)計(jì)算區(qū)域?yàn)?90 m×130 m×124 m的長方體區(qū)域(網(wǎng)格為38×26×25)，整個(gè)區(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個(gè)電極的中心，在y方向上位于兩巷道的中間位置。

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

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

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

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

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

2 反演成像理論

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

2.1 反演方程

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

(4)

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

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

(5)

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

(6)

因?yàn)閐cal(mk+1)=dcal(mk)+J·(Δmk)，所以式(6)可寫為：

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

(7)

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

2.2 雅克比矩陣J的處理

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

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

(8)

所以

(9)

在正演過程中由于

(10)

因此

(11)

(12)

(13)

所以

(14)

其中

(15)

所以兩邊求導(dǎo)得

(16)

由于A為對(duì)稱矩陣，所以

(17)

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

(18)

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

2.3 共軛梯度算法

3 數(shù)值模型實(shí)驗(yàn)

3.1 實(shí)驗(yàn)?zāi)Ｐ?/h3>

建立三維空間的含水陷落柱模型，將觀測(cè)區(qū)域剖分為14 m×14 m×12 m的網(wǎng)格，設(shè)陷落柱模型規(guī)格為2 m×2 m×2 m，電阻率值為10 Ω·m ，底板圍巖電阻率為100 Ω·m ，覆蓋煤層電阻率為2 000 Ω·m 。陷落柱的中心位置坐標(biāo)為(0 m，0 m，-1 m)，兩平行巷道間距離為9 m，左側(cè)為一連通巷道，空間位置如圖3所示，圖3(b)和圖3(c)分別為導(dǎo)水陷落柱在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é)果進(jìn)行切片，如圖4所示。藍(lán)色低阻區(qū)域電阻率小于60 Ω·m ，呈橢球狀，反演結(jié)果可以很好地反映含水體的位置。

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

4 物理實(shí)驗(yàn)及結(jié)果分析

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)論與建議

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

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

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

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

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