孟 麟， 王 智

(長江大學(xué)電子信息學(xué)院， 荊州 434023)

在實(shí)際煤炭的探勘和開采過程中，陷落柱是不可避免會(huì)遇到的一種災(zāi)害性的地質(zhì)異常體。它是由于地下水不斷對上方可溶性巖石進(jìn)行溶蝕和上方覆蓋巖層重力的共同作用下產(chǎn)生的溶巖塌陷體，給煤田的開采工作帶來了極大的安全隱患。為了預(yù)防安全隱患，準(zhǔn)確勘探出陷落柱位置是最有效的方法，許多學(xué)者利用陷落柱和周圍圍巖之間的密度、磁性和電性差異，在地表采用地震勘探[1]、瞬變電磁法[2]、直流電法[3]和高密度電法[4]等方法提前進(jìn)行陷落柱的勘探。受地表?xiàng)l件限制和陷落柱空間形態(tài)特征的影響，目前準(zhǔn)確高效地探測陷落柱位置信息還存在一定的困難。特別是對于直流電法勘探來說，地表?xiàng)l件會(huì)對勘探的深度產(chǎn)生比較大的影響，且直流電法勘探的精度會(huì)隨著地質(zhì)體埋深的增加而降低。為了解決勘探深度的問題，井-地聯(lián)合電阻率觀測裝置成為了較好的選擇。相比與井-井(或地-地)觀測裝置，井-地聯(lián)合觀測裝置對于地下深部的勘探有著異常幅值和勘探精度都更高的優(yōu)點(diǎn)，在勘探過程中有著很好的應(yīng)用前景。近年來，井-地聯(lián)合觀測裝置已在勘察領(lǐng)域得到廣泛應(yīng)用。湯井田等[5]利用井-地電阻率法得到歧離率來確定高阻率油藏邊界；趙廣茂等[6]用井-地觀測裝置來勘探地下陷落柱的存在和含水性；王智等[7]利用井-地觀測方法和改進(jìn)的異常電位法，確定了適用于起伏地表的自然邊界條件。

雖然井-地聯(lián)合電阻率觀測裝置已經(jīng)被廣泛應(yīng)用在實(shí)際工程勘探中，理論研究方法也存在許多，但是對于在起伏地表的情況下不同形態(tài)和位置的典型異常體的異常形態(tài)特征的系統(tǒng)研究不足。使用井-地聯(lián)合電阻率觀測裝置針對起伏地表情況下，含水陷落柱處于不同形態(tài)和位置所產(chǎn)生異常情況進(jìn)行系統(tǒng)的分析，以期為實(shí)際勘探含水陷落柱提供理論依據(jù)。

1 2.5D有限單元法數(shù)值模擬

在三維無限半空間Ω中，建立直角坐標(biāo)系。假設(shè)Z軸垂直向下，那么點(diǎn)電源O的坐標(biāo)為(xO,yO,zO),表示為O，則總電位滿足微分方程為[8]

(1)

(2)

式(2)中：u(x,y,z)簡寫為u；ΓS和Γ分別為三維無線半空間Ω的地表邊界和無窮遠(yuǎn)邊界；n為地表邊界ΓS的外法向向量；r表示該坐標(biāo)點(diǎn)到點(diǎn)電源之間的距離。

假設(shè)電性的分布沿Y軸方向不變，對式(1)、式(2)進(jìn)行傅氏變換，即將空間域的三維邊界問題變成波數(shù)域的二維邊界問題[8]：

(3)

式(3)中：U(x,ki,z)表示波數(shù)域的電位，簡寫為U，其中ki為波數(shù)(i=1,2,…,n)；K1和K0分別為一階和零階貝塞爾函數(shù)。式(3)對應(yīng)的變分問題[8]為

δF(U)=0

(4)

采用不規(guī)則三角形單元e對求解區(qū)域進(jìn)行剖分后，F(xiàn)e(U′)表示每個(gè)三角形單元上波數(shù)域電位U′的函數(shù)，并將其進(jìn)行總體合成后得

(5)

KU=P

(6)

解該方程組，便可求得節(jié)點(diǎn)m在對應(yīng)波數(shù)k中的電位U(m,ki)。然后按式(7)進(jìn)行傅里葉反變換便可得到空間中的電位u(m)：

(7)

式(7)中：波數(shù)(ki)與傅里葉反變換系數(shù)(gi)選自徐世浙[8]采用最優(yōu)化方法得到的5個(gè)波數(shù)和對應(yīng)的傅里葉反變化系數(shù)。

2 觀測方式與觀測裝置

井中和地表電極的聯(lián)合觀測方法可以分為：井-地聯(lián)合觀測方式和地-井聯(lián)合觀測方式。井-地聯(lián)合觀測方式是指將供電電極置于井中，測量電極置于地表的觀測方式；地-井聯(lián)合觀測方式是指將供電電極置于地表，測量電極置于井中的觀測方式。兩種觀測方式均可以通過分別移動(dòng)供電電極和測量電極和相對位置來研究不同水平位置和不同深度地下介質(zhì)的視電阻率分布情況。再根據(jù)在實(shí)際勘探工作中所使用的不同觀測裝置，可以將其細(xì)分為井-地二極裝置、地-井二極裝置、井-地三極裝置、地-井三極裝置、井-地四極裝置、地-井四極裝置等。因篇幅有限，研究選擇井-地二極裝置為觀測裝置。

井-地二極測深裝置，指供電電極A置于井中，另一供電電極B置于無窮遠(yuǎn)處；觀測電極M置于地表，另一觀測電極N置于無窮遠(yuǎn)處。再通過上下移動(dòng)供電電極A和水平移動(dòng)觀測電極M，便可以探測不同深度和水平位置的地下介質(zhì)分布。以M的橫坐標(biāo)為觀測記錄點(diǎn)的橫坐標(biāo)，A的縱坐標(biāo)為觀測記錄點(diǎn)的縱坐標(biāo)畫出測深斷面圖，觀測示意圖如圖1所示。

圖1 井-地二極觀測裝置示意[9]Fig.1 Sketch map of pole-pole on borehole-ground joint resistivity surveying system[9]

3 數(shù)值算例及分析

3.1 算法驗(yàn)證

在完成有限單元法中網(wǎng)格剖分的部分時(shí)，使用Gmsh軟件[10]來建立典型的地質(zhì)模型并完成局部加密的非結(jié)構(gòu)化網(wǎng)格剖分。在存儲(chǔ)大型稀疏矩陣時(shí)使用CSR(compressed sparse row)格式[11]，減少程序?qū)?nèi)存的占用。對于求解正演過程中的線性方程組，選用Eigen庫中的BiCGSTAB(biconjugate gradient stabilized method)求解器。在計(jì)算多個(gè)點(diǎn)電源時(shí)，使用MPI(message passing interface)進(jìn)行并行計(jì)算，提高計(jì)算效率。為了驗(yàn)證算法的正確性，首先選擇具有解析解的簡單模型進(jìn)行計(jì)算。

二層介質(zhì)的地電斷面模型如圖2所示。圖2中，ρ1、ρ2分別表示第一層與第二層介質(zhì)的電阻率，第一層斷面的高度為15 m。表2給出了單個(gè)點(diǎn)電源(O)時(shí)，地表不同測點(diǎn)(V)的電位。

由表1數(shù)據(jù)表明，計(jì)算解的最大相對誤差發(fā)生在距離點(diǎn)源最近的測量點(diǎn)上達(dá)到5%，其余測點(diǎn)的相對誤差均在1%左右，證明了編譯算法的正確性。

3.2 算例

使用井-地二級裝置進(jìn)行觀測。將鉆井井口與地質(zhì)體重心連線在地表上的投影，記為X測線。井中供電電極A′和地表觀測電極M′的點(diǎn)距均設(shè)置為1 m。測量電極在地表X測線0～50 m移動(dòng)。鉆井井口在X測線14 m處垂直向下，記為Z測線，供電電極在Z測線范圍5～20 m內(nèi)移動(dòng)。

圖2 二層介質(zhì)地電斷面Fig.2 Two-layers geo-electric cross section

陷落柱是煤層下伏可溶性巖石被地下水溶蝕后引起上覆巖層冒落而形成的古巖溶塌陷體。含水陷落柱在電阻率特性上顯示低阻特性，為了系統(tǒng)地對比各種裝置和不同位置模型的觀測視電阻率異常，圍巖電阻率ρ1為100 Ω·m，含水陷落柱電阻率ρ2為5 Ω·m。模擬的含水陷落柱的模型類似為“子彈形”，頂部埋深為8 m，底部深度16 m。通過改變模型中心距離點(diǎn)源所在的鉆井的垂直距離來分析其對觀測視電阻率異常的影響，設(shè)定的距離分別為11、21 m。若含水陷落柱模型頂點(diǎn)的位置確定時(shí)，改變模型的傾斜角度，分別為左傾30°、直立90°和右傾30°，對比傾斜角度的不同對觀測視電阻率異常的影響。地質(zhì)模型的網(wǎng)格剖分采用局部加密的方式，平坦地形下直立含水陷落柱模型及網(wǎng)格剖分如圖3所示。

3.2.1 平坦地形

由圖4可知，當(dāng)在平坦地形的情況下含水陷落柱模型的視電阻率斷面圖中異常只存在獨(dú)立的橢圓形的低阻異常，存在這樣的低阻異常主要是因?yàn)楹萋渲Ｐ拖鄬鷰r電阻率為低阻，在電場中吸引電流。在X測線方向，不論距離鉆井的距離和傾斜角度如何變化，低阻異常中心始終位于含水陷落柱模型的右側(cè)；在Z測線方向，隨著測點(diǎn)深度的增加，在含水陷落柱模型的上部分的深度區(qū)間中，在水平方向上的等值線相對于模型下部分的深度區(qū)間中的等值線更加密集，在垂直方向等值線的變化較為平緩。

表1 算法所得的電位計(jì)算解與解析解

圖3 平坦地形下直立含水陷落柱模型示意Fig.3 The upright water-bearing collapse column model in flat terrain

圖4 含水陷落柱模型視電阻率斷面(平坦地形)Fig.4 The resistivity section of water-bearing collapse column model(flat terrain)

對比圖4(a)、圖4(c)、圖4(e)和圖4(b)、圖4(d)、圖4(f)可以看出，低阻異常的最大值和距離鉆井的距離有關(guān)，在含水陷落柱模型距離鉆井11 m的位置時(shí)，低阻異常的最大值為13%左右；在含水陷落柱模型距離鉆井21 m位置時(shí)，低阻異常的最大值減小到10%左右。同時(shí)模型在距離鉆井11 m的位置時(shí)，低阻異常中心的深度與含水陷落柱模型頂端的深度一致，水平方向上差距2～4 m，當(dāng)距離鉆井的距離為21 m時(shí)，低阻異常的中心發(fā)生上移，但在水平方向上的差距基本不變化。

在含水陷落柱模型距離鉆井的位置發(fā)生改變時(shí)，改變模型傾斜角度所產(chǎn)生影響的趨勢基本不發(fā)生變化。當(dāng)模型左傾30°時(shí)，低阻異常被向右上方擠壓；當(dāng)模型右傾30°時(shí)，低阻異常被向左下方拉伸。但無論如何改變傾斜角度，在含水陷落柱模型底部向下方向的視電阻率等值線是最為稀疏的。

圖5 純山谷地形視電阻率斷面Fig.5 The resistivity section of valley terrain

3.2.2 山谷地形

由圖5可見，純山谷地形的視電阻率斷面圖會(huì)產(chǎn)生高阻、低阻異常并存的情況，均呈橢圓形分布，與前人研究得出的地形與視電阻率呈左右對稱的規(guī)律相符，這是由于在山谷地形中存在空氣，阻礙了電流的通過。低阻異常位于山谷地形的右下方，當(dāng)山谷地形臨近鉆井時(shí)，低阻異常中心距離山谷地形更近，位于X=28～30 m，Z=5 m處，異常最大值為22%；若山谷地形遠(yuǎn)離鉆井時(shí)，低阻異常中心向下偏移，位于X=38～40 m，Z=6 m處，異常最大值為14%。高阻異常的中心位置在山谷地形的左下方，臨近鉆井時(shí)，其中心位于X=14～15 m，Z=5～6 m處，異常最大值為14%；遠(yuǎn)離鉆井時(shí)，其中心略微向右偏移，位于X=25～26 m，Z=5 m處，異常最大值為11%。在圖5中等值線最密集的地方都位于山谷地形底端的中心處。

根據(jù)圖6(a)、圖6(c)、圖6(e)和圖5(a)對比，即在臨近鉆井的純山谷地形下方存在含水陷落柱模型。視電阻率斷面圖中依舊是高阻和低阻異常并存的情況，改變傾斜角度對高阻和低阻異常的最大值并無影響。對比純山谷地形，此時(shí)低阻異常的最大值為36%，高阻異常的最大值改變?yōu)?0%。當(dāng)含水陷落柱模型左傾30°時(shí)，模型位于山谷右部低阻異常下方，視電阻率異常基本無變化；當(dāng)含水陷落柱模型直立時(shí)，由于模型的低阻特性吸引電流，是視電阻率等值線向高阻異常方向擠壓。模型的上部分面積小，相對應(yīng)地對電流的吸引力小，等值線變化小，下部分面積相對大，對等值線影響增大；當(dāng)含水陷落柱模型右傾30°時(shí)，斷面圖視電阻率異常不再顯示為高阻異常和低阻異常左右對稱，在模型的影響下阻斷了高阻異常對下方的影響，產(chǎn)生了一個(gè)新的低阻異常，異常中心與模型的最低點(diǎn)平行，最大值為6%。

根據(jù)圖6(b)、圖6(d)、圖6(f)和圖5(a)對比，低阻異常的最大值變?yōu)?0%，高阻異常的最大值變?yōu)?2%，斷面圖的視電阻率等值線趨勢基本無變化，是因?yàn)樯焦鹊匦螌ο路胶萋渲Ｐ彤a(chǎn)生的異常產(chǎn)生了遮擋效果。

根據(jù)圖7(a)、圖7(c)、圖7(e)和圖5(b)對比，即在遠(yuǎn)離鉆井的純山谷地形下方存在一個(gè)含水陷落柱模型。與純山谷地形產(chǎn)生的一個(gè)低阻異常和一個(gè)高阻異常不同，此時(shí)視電阻率斷面上顯示兩個(gè)低阻異常和一個(gè)高阻異常，從左到右依次為低阻-高阻-低阻，分別是由含水陷落柱模型阻斷高阻異常產(chǎn)生和山谷地形產(chǎn)生。純山谷地形產(chǎn)生的高阻和低阻異常中心位置不變動(dòng)，低阻異常最大值改變?yōu)?2%，高阻異常的最大值改變?yōu)?%，基本和圍巖視電阻率相等。由含水陷落柱模型阻斷山谷高阻異常而產(chǎn)生的低阻異常的中心位于模型的左側(cè)，在Z測線方向和模型的中心平行，在X測線方向和模型中心距離5 m。由于含水陷落柱在有低阻特性，對電流有吸引作用使山谷地形中部垂直X測線的視電阻率等值線向模型方向凸起，而偏移角度基本不影響此趨勢。當(dāng)存在偏移角度時(shí)對由模型產(chǎn)生的低阻異常的形態(tài)有所影響，當(dāng)右傾30°時(shí)，對低阻異常由向左下方的拉伸趨勢；當(dāng)左傾30°時(shí)，對低阻異常向右方拉伸趨勢。

根據(jù)圖7(b)、圖7(d)、圖7(f)和圖6(b)、圖6(d)、圖6(f)對比可得，當(dāng)遠(yuǎn)離鉆井即遠(yuǎn)離測點(diǎn)電源的時(shí)候，含水陷落柱模型對純山谷地形產(chǎn)生的影響減弱，但基本的趨勢仍然存在。

在Z測線方向上，5～20 m區(qū)域內(nèi)為觀測區(qū)域，此區(qū)域使用Surfer軟件根據(jù)數(shù)據(jù)自動(dòng)填充圖6 含水陷落柱模型視電阻率斷面(山谷地形臨近鉆井)Fig.6 The resistivity section of water-bearing collapse column model(valley terrain near drilling)

在Z測線方向上，5～20 m區(qū)域內(nèi)為觀測區(qū)域，此區(qū)域使用Surfer軟件根據(jù)數(shù)據(jù)自動(dòng)填充圖7 含水陷落柱模型視電阻率斷面(山谷地形遠(yuǎn)離鉆井)Fig.7 The resistivity section of water-bearing collapse column model(valley terrain away from drilling)

3.2.3 山脊地形

由圖8可知，在純山脊地形的視電阻率斷面圖中產(chǎn)生了高阻、低阻異常并存的情況，地形右邊存在高阻異常，左邊存在低阻異常與純山谷地形產(chǎn)生的視電阻率異常正好相反，這也符合前人的研究結(jié)果。當(dāng)山脊地形臨近鉆井時(shí)，低阻異常中心位于X=15 m，Z=5 m處，異常值最大為10%，高阻異常中心位于X=26～27 m，Z=5 m處，異常最大值為17%；當(dāng)山脊地形遠(yuǎn)離鉆井時(shí)，低阻異常中心位于X=26～27 m，Z=5 m處，異常的最大值為10%左右，高阻異常中心位于X=37～38 m，Z=10 m處，異常值最大為10%。在改變山脊地形與鉆井的距離的過程中，對低阻異常造成的影響不大，但是使高阻異常發(fā)生了明顯的改變，使其向下拉伸的趨勢嚴(yán)重且產(chǎn)生的高阻異常最大值也下降近1/2。

根據(jù)圖8(a)、圖9可得，在臨近鉆井的純山脊地形下方存在含水陷落柱模型時(shí)無論模型離鉆井的距離和傾斜角度如何變化，在斷面中都顯示兩個(gè)低阻異常和一個(gè)高阻異常，從左到右依次為低阻-高阻-低阻。

由圖9(a)、圖9(c)、圖9(e)可以看出，在含水陷落柱模型距離鉆進(jìn)11 m時(shí)，左部低阻異常最大值為13%，同時(shí)低阻異常中心向下拉伸，在Z測線方向基本與異常體模型的中心平行，含水陷落柱模型傾斜角度為右傾30°時(shí)，低阻異常有向左拉伸的趨勢。中部高阻異常最大值下降為7%，此時(shí)含水陷落柱模型將高阻異常向下的趨勢阻斷。右部低阻異常最大值為7%，低阻異常的中心隨含水陷落柱模型的角度傾斜而改變，左傾30°時(shí)低阻異常中心向下偏移；直立狀態(tài)時(shí)低阻異常中心和含水陷落柱模型的中心平行；右傾30°時(shí)低阻異常中心向上偏移。由圖9(b)、圖9(d)、圖9(f)可以看出，當(dāng)含水陷落柱模型距離鉆井21 m時(shí)，左部低阻異常的最大值和形態(tài)與純山脊地形對比都無明顯區(qū)別。中部高阻異常最大值下降到12%，右部低阻異常的最大值為7%。右部低阻異常中心隨含水陷落柱模型傾斜角度的變化與圖9(a)、圖9(c)、圖9(e)顯示的變化趨勢相同。

根據(jù)圖8(b)、圖10可得，在遠(yuǎn)離鉆井的純山脊地形下方存在含水陷落柱模型的斷面圖與遠(yuǎn)離鉆井的純山脊地形的斷面圖相比，都存在一個(gè)左部的低阻異常和一個(gè)右部的高阻異常，不同之處在于存在右部的高阻異常由于受到異常體模型的影響其異常中心向下偏移嚴(yán)重。

由圖10(a)、圖10(c)、圖10(e)可以看出，在含水陷落柱模型距離鉆井11 m時(shí)，左部的低阻異常的最大值上升到19%，右部高阻異常的最大值下降到2%。當(dāng)含水陷落柱模型右傾30°和直立時(shí)，視電阻率等值線變化不明顯，但是含水陷落柱模型左傾30°時(shí)，低阻異常左部視電阻率等值線向左拉伸趨勢明顯。由圖10(b)、圖10(d)、圖10(f)可以看出，在含水陷落柱模型距離鉆井21 m時(shí)，左部的低阻異常的最大值上升到14%，右部高阻異常的最大值下降到3%。相比左傾30°和右傾30°，在直立狀態(tài)時(shí)高阻異常中心下降的趨勢最大。

通過對井-地二極裝置各種地形情況和不同位置的含水陷落柱模型的正演結(jié)果分析，可以得到如下結(jié)果。

(1)在平坦地形情況下，含水陷落柱模型會(huì)產(chǎn)生一個(gè)低阻異常。隨著模型的位置和傾斜角的改變，視電阻率斷面圖中的低阻異常的形態(tài)和最大值也跟著發(fā)生改變。

(2)在純山谷或純山脊地形的情況下，在視電阻率斷面圖中從左到右依次產(chǎn)生高阻-低阻或低阻-高阻異常。隨著地形距離鉆井位置改變，視電阻率斷面圖中顯示的趨勢不變，但是產(chǎn)生異常的形態(tài)和最大值發(fā)生改變。

在Z測線方向上，5～20 m區(qū)域內(nèi)為觀測區(qū)域，此區(qū)域使用Surfer軟件根據(jù)數(shù)據(jù)自動(dòng)填充圖8 純山脊地形視電阻率斷面Fig.8 The resistivity section of ridge terrain

在Z測線方向上，5～20 m區(qū)域內(nèi)為觀測區(qū)域，此區(qū)域使用Surfer軟件根據(jù)數(shù)據(jù)自動(dòng)填充圖9 含水陷落柱視電阻率斷面(山脊地形臨近鉆井)Fig.9 The resistivity section of water-bearing collapse column model(ridge terrain near drilling)

在Z測線方向上，5～20 m區(qū)域內(nèi)為觀測區(qū)域，此區(qū)域使用Surfer軟件根據(jù)數(shù)據(jù)自動(dòng)填充圖10 含水陷落柱視電阻率斷面(山脊地形遠(yuǎn)離鉆井)Fig.10 The resistivity section of water-bearing collapse column model(ridge terrain away from drilling)

(3)當(dāng)含水陷落柱模型處于山谷地形或山脊地形的產(chǎn)生的低阻異常區(qū)域時(shí)，模型在視電阻率斷面圖中顯示的異常不明顯。起伏地形對其產(chǎn)生了遮擋作用。

(4)當(dāng)含水陷落柱模型處于山谷地形或山脊地形的產(chǎn)生的高阻異常區(qū)域時(shí)，視電阻率斷面圖中產(chǎn)生低阻-高阻-低阻的異常情況，根據(jù)除地形產(chǎn)生的兩個(gè)異常外的另一個(gè)低阻異?？梢耘袛喈惓ｓw的形態(tài)和位置信息。

4 結(jié)論

針對起伏地形條件下，存在含水陷落柱對地下產(chǎn)生的異常無系統(tǒng)分析的問題，采用有限單元法對井-地電極聯(lián)合觀測裝置電阻率測深的二維地質(zhì)體視電阻率斷面進(jìn)行數(shù)值計(jì)算。共包含平坦地形、山谷地形和山脊地形三種地形條件，在兩種起伏地形中又討論了其距離鉆井不同位置的影響。每種地形條件下均設(shè)置了距離鉆井11 m或21 m，傾斜角度為左傾30°、直立90°和右傾30°共6種情況的含水陷落柱模型。詳細(xì)地分析了在這些不同情況下，含水陷落柱模型在視電阻率斷面中產(chǎn)生的影響，研究結(jié)果為實(shí)際含水陷落柱勘察過程提供了一定的理論依據(jù)。