張淼淼， 石顯新

(1.煤炭科學(xué)研究總院， 北京 100013； 2.中國煤炭科工集團(tuán)西安研究院有限公司， 西安 710077)

直流電阻率法是以巖層中巖石和礦石的電性差異為物質(zhì)基礎(chǔ)，由電極向地下供入電流，通過觀測與研究人工建立的地中電流場的分布規(guī)律來進(jìn)行地質(zhì)勘探的一種電法勘探方法[1]。直流電阻率勘探無論是在金屬、非金屬礦產(chǎn)普查和地質(zhì)構(gòu)造研究方面，還是在水文地質(zhì)調(diào)查、工程地質(zhì)調(diào)查等方面，均發(fā)揮著重要的作用[2]。劉偉等[3]對低阻富水構(gòu)造的多極距聯(lián)合剖面曲線進(jìn)行了二維正演模擬，對模型在二維中的傾斜狀態(tài)進(jìn)行了闡述。王志鵬等[4]利用高密度電法對斷層進(jìn)行探測，該方法對斷層識別具有一定的有效性和準(zhǔn)確性，但是對逆斷層的探測效果明顯好于正斷層。孟麟等[5]采用井地電阻率法對起伏地形下陷落柱模型進(jìn)行正演模擬，結(jié)果表明地形會改變陷落柱異常響應(yīng)的形態(tài)和位置。通過上述研究可以發(fā)現(xiàn)，對于電阻率法探測斷層，主要是基于在二維空間下的單一方位變化進(jìn)行研究，而在其他地質(zhì)異常體的研究中也主要是對地形變化對探測結(jié)果的影響。所以在三維空間下異常體不同方位對探測結(jié)果的影響需要進(jìn)一步研究。

在運(yùn)用直流電阻率勘探中，準(zhǔn)確確定異常體實(shí)際位置與校正其位置偏差都非常重要，為了提高勘探效果，就需要對地質(zhì)異常體的影響因素進(jìn)行研究，其中對地質(zhì)異常體的產(chǎn)狀變化也很有必要，對于二維空間下異常體的空間變化只是單一方向的變化，而在實(shí)際探測中是具有不同空間方位的變化的。因此，現(xiàn)采用COMSOL有限元軟件進(jìn)行三維空間下不同方位的異常體進(jìn)行直流電阻率研究，首先對點(diǎn)電源進(jìn)行直流電法正演模擬誤差驗(yàn)證，然后用板狀體模擬斷層，采用單極-偶極與偶極-偶極裝置研究低阻板狀體在空間方位上變化的異常響應(yīng)特征，最后結(jié)合工程實(shí)例進(jìn)一步分析驗(yàn)證，以期為直流電阻率法后期處理中準(zhǔn)確識別斷層等提供理論依據(jù)。

1 穩(wěn)定電流場的基本理論

在穩(wěn)定電流場的任意點(diǎn)上，電流密度矢量j與電場強(qiáng)度矢量E，則介質(zhì)電阻率ρ與電流密度，電位U和電場強(qiáng)度之間有著如下關(guān)系[2]：

(1)

(2)

(3)

將式(1)和式(2)代入式(3)，則可以得

(4)

在無源均勻介質(zhì)空間中，ρ為常數(shù)，則式(4)為

(5)

式(5)即拉普拉斯方程，它表示穩(wěn)定電流場中不包含電流源在內(nèi)的空間任何一點(diǎn)的電位只是空間位置的函數(shù)，在直角坐標(biāo)系(x，y，z)中拉普拉斯方程表達(dá)式為

(6)

在球(r，θ，φ)坐標(biāo)系中，有

(7)

在有源的均勻介質(zhì)空間中，ρ為常數(shù)，有一點(diǎn)電源A(x0，y0，z0)，其電流強(qiáng)度為I，通過引入狄拉克函數(shù)δ(A)，得到滿足電位與電流之間關(guān)系的微分方程，即泊松方程形式：

(8)

式(8)中：δ(A)是以點(diǎn)A(x0，y0，z0)為中心的狄拉克函數(shù)。

求出滿足給定邊界條件的泊松方程的解，即穩(wěn)定電流場的邊值問題，如圖1所示，在求解電場問題時，邊值問題有三種[6-8]。

在地面邊界Γs上，電位的法向?qū)?shù)為

(9)

在無窮遠(yuǎn)邊界Γ∞上，電位是點(diǎn)源電位：

(10)

式(10)中：c為常數(shù)；r為點(diǎn)源到邊界距離。對式(10)求偏導(dǎo)，消去常數(shù)c，得

(11)

Γs為地面邊界；?！逓闊o窮遠(yuǎn)邊界；r為任意點(diǎn)A到無窮遠(yuǎn)邊界?！?的距離圖1 點(diǎn)電源邊界示意圖Fig.1 Sketch of point electric source boundary

2 誤差分析

COMSOL軟件是一款以有限元法為基礎(chǔ)，通過求解偏微分方程的多物理場仿真軟件，用來分析電磁學(xué)、結(jié)構(gòu)力學(xué)、聲學(xué)、流體流動和化工等眾多領(lǐng)域的實(shí)際工程問題，COMSOL核心產(chǎn)品可以輕松實(shí)現(xiàn)建模流程的各個環(huán)節(jié)，為科研人員從常規(guī)有限元編程中解脫出來，COMSOL還擁有強(qiáng)大的求解器，可自動檢測要求解的物理場，以及問題大小，并針對具體問題選擇求解器，可以直接或迭代求解。在COMSOL中，用戶可以自由控制模型的各個方面，其主要功能提供了基本的幾何建模工具，也可以從AutoCAD制圖軟件導(dǎo)入模型。此外，還可以直接在圖形用戶界面使用方程式和表達(dá)式中輸入用戶自定義參數(shù)來進(jìn)行仿真。COMSOL軟件是基于有限元分析來模擬的，所以主要的分析過程與有限元大致相同，主要有三個步驟：前處理過程、求解過程和后處理過程。COMSOL分析的流程圖如圖2所示。

COMSOL的前處理階段主要是建立幾何模型、給定材料性質(zhì)、設(shè)置邊界條件、網(wǎng)格剖分。在COMSOL模擬仿真中，其模型建立的好壞將直接影響收斂的速度以及精度，所以在建立幾何模型時，應(yīng)該對實(shí)際模型進(jìn)行簡化處理，這樣會減少求解的時間。在設(shè)置邊界條件中，在距離點(diǎn)電源無窮遠(yuǎn)處可以假設(shè)電位為零，這樣就可以滿足第一類邊界條件。但由于在COMSOL中建立太遠(yuǎn)邊界，將會使得計算范圍增大，從而導(dǎo)致計算時間太長，所以通過在模型周圍建立無限元域即可防止模型過大。在地表邊界處，空氣的電阻率無限大而導(dǎo)致電流不能穿過，所以在模型中設(shè)置地表面電流密度的法向分量為零，滿足第二類邊界條件[9]。在COMSOL軟件中無需設(shè)置內(nèi)部連續(xù)條件，軟件自動滿足。

在有限元建模中網(wǎng)格的選擇將會決定模擬的精度，在COMSOL中通過自適應(yīng)網(wǎng)格優(yōu)化解決這個問題，它首先在初始網(wǎng)格上求解，然后通過迭代將元素插入估計誤差高的區(qū)域，然后重新求解模型，這樣可以最大限度地減少模型中的總體誤差。在COMSOL 軟件中，可以根據(jù)需要來創(chuàng)建自由網(wǎng)格、掃掠網(wǎng)格和映射網(wǎng)格等，利用網(wǎng)格剖分工具，可以生成三角形、四邊形、四面體和六面體等網(wǎng)格單元，而四面體網(wǎng)格的適應(yīng)性較強(qiáng)，且容易實(shí)現(xiàn)局部網(wǎng)格加密，所以一般采用四面體網(wǎng)格剖分[10]。COMSOL有限元軟件中的全局自適應(yīng)網(wǎng)格細(xì)化方法通過使用誤差估計，確定建模域中局部誤差最大的點(diǎn)，在兼顧整個模型的局部誤差的同時，同時在局部誤差更顯著的區(qū)域使用較小的單元。這種方法的好處在于，所有的網(wǎng)格細(xì)化工作均自動完成；其缺點(diǎn)是無法對網(wǎng)格進(jìn)行控制。

圖2 COMSOL分析流程圖Fig.2 Chart of COMSOL analysis

首先為了驗(yàn)證 COMSOL 軟件對于地電模型正演的可靠性，通過建立三維均勻半空間的簡單地電模型對其進(jìn)行數(shù)值模擬。建立一個半徑為10 m的半球體模型模擬均勻半空間，設(shè)置均勻半空間的背景電阻率為ρ=100 Ω·m，供電電極坐標(biāo)(0，0，0)處，電流強(qiáng)度I=1 A，求取數(shù)值解并將與理論解進(jìn)行對比。在三維均勻半空間中，點(diǎn)電源電位U解析公式為

(12)

式(12)中：ρ為背景電阻率；R為觀測距離。

表1給出了解析解和數(shù)值解的對比，從表1可知，在距離點(diǎn)源1 m處，相對誤差為6.293 1%；在距離點(diǎn)源5 m處，相對誤差為0.404 7%；而在距離點(diǎn)源10～50 m處，相對誤差都小于0.100 0%。可以看出，在距離點(diǎn)電源較近的范圍(小于3 m)，誤差相對較大，越遠(yuǎn)離點(diǎn)電源，其誤差越小。所以，可以看出COMSOL 軟件對于直流電法三維地電模型的數(shù)值模擬結(jié)果是可靠的。

表1 三維點(diǎn)電源模型電位誤差對比Table 1 Comparison of potential errors in three-dimensional point power supply models

3 低阻板狀體直流電阻率模擬

3.1 板狀體不同旋轉(zhuǎn)方向數(shù)值模擬

3.1.1 模型1

設(shè)水平均勻大地電阻率背景值ρ1=100 Ω·m，在水平地面下方有一電阻率為ρ2=10 Ω·m的板狀低阻異常體，模擬含水構(gòu)造異常體，截面積10 m×10 m，厚度2 m，以板狀異常體中心為圓心，沿Y軸順時針旋轉(zhuǎn)，旋轉(zhuǎn)角度為0°、30°、60°、90°。根據(jù)模型采用單極-偶極與偶極-偶極裝置進(jìn)行觀測計算。首先建立的模型一，原點(diǎn)位于地面正中心。其次先采用單極-偶極裝置在地面A點(diǎn)(-70，0，0)供電，電流強(qiáng)度為1 A，M、N為測量電極，極距為5 m；然后采用偶極-偶極裝置，A(-70，0，0)和B(-65，0，0)兩點(diǎn)供電，電流強(qiáng)度分別為1 A和-1 A，極距5 m，兩種方式測線距離為-60～60 m與-55～60 m。模型如圖3所示。

圖4給出了當(dāng)埋深h=10 m時采用單極-偶極沿Y軸旋轉(zhuǎn)的正演結(jié)果視電阻率曲線。從圖4中的4條曲線可以看出，除曲線前端位置出現(xiàn)不同程度的波動以外，整體呈現(xiàn)出不對稱形態(tài)。當(dāng)長方體未進(jìn)行旋轉(zhuǎn)時即Y軸旋轉(zhuǎn)0°曲線中可以看出，在極值異常兩側(cè)呈現(xiàn)左高右低的雙峰形態(tài)，極值異常最??；隨著順時針旋轉(zhuǎn)角度的增大，在旋轉(zhuǎn)30°、60°、90°時極值異常兩側(cè)都呈現(xiàn)出左低右高的雙峰形態(tài)，極值異常逐漸增大。正演結(jié)果顯示，當(dāng)板狀異常體沿Y軸順時針旋轉(zhuǎn)時，雙峰的高低與旋轉(zhuǎn)方向呈現(xiàn)對應(yīng)關(guān)系；且極值異常大小隨著旋轉(zhuǎn)角度的增大，異常幅值也在增大。

圖3 模型1示意圖Fig.3 Schematic diagram of model 1

圖4 模型1單極-偶極裝置視電阻率曲線Fig.4 Apparent resistivity curve of pole-dipole device in model 1

圖5給出了當(dāng)埋深h=10 m時采用偶極-偶極裝置沿Z軸旋轉(zhuǎn)的正演結(jié)果視電阻率曲線。從圖5中的4條曲線可以看出，除繞Z軸旋轉(zhuǎn)30°曲線前端位置出現(xiàn)波動以外，其余3條曲線較光滑，而整體呈現(xiàn)出不對稱曲線形態(tài)。而正演結(jié)果與單極-偶極模擬結(jié)果很相似，當(dāng)板狀異常體沿Z軸順時針旋轉(zhuǎn)時，雙峰的高低與旋轉(zhuǎn)方向呈現(xiàn)對應(yīng)關(guān)系；極值異常大小隨著旋轉(zhuǎn)角度的增大，異常幅值也隨之增大。

3.1.2 模型2

以板狀異常體中心為圓心，沿Z軸順時針旋轉(zhuǎn)，旋轉(zhuǎn)角度為0°、30°、60°、90°，其余設(shè)置與模型1一致。根據(jù)模型采用單極-偶極與偶極-偶極裝置進(jìn)行觀測計算，模型2如圖6所示。

圖7給出了當(dāng)埋深h=10 m時采用單極-偶極沿Y軸旋轉(zhuǎn)的正演結(jié)果視電阻率曲線。從圖7中的4條曲線可以看出，4條曲線前端位置都出現(xiàn)了一定波動，整體呈現(xiàn)出平滑的不對稱曲線。當(dāng)長方體沿Z軸旋轉(zhuǎn)0°曲線中可以看出，在極值異常兩側(cè)呈現(xiàn)左高右低的雙峰形態(tài)；隨著順時針旋轉(zhuǎn)角度的增大，在旋轉(zhuǎn)30°、60°、90°時極值異常兩側(cè)呈現(xiàn)出左低右高的雙峰形態(tài)，極值逐漸減小。正演結(jié)果顯示，當(dāng)板狀異常體沿Y軸順時針旋轉(zhuǎn)時，雙峰的高低與旋轉(zhuǎn)方向呈現(xiàn)對應(yīng)關(guān)系，隨著順時針旋轉(zhuǎn)角度增大，極值左邊的單峰逐漸減?。磺覙O值異常大小隨著旋轉(zhuǎn)角度的增大，異常幅值也在增大。

圖5 模型1偶極-偶極裝置視電阻率曲線Fig.5 Apparent resistivity curve of dipole-dipole device in model 1

圖6 模型2示意圖Fig.6 Schematic diagram of model 2

圖8給出了當(dāng)埋深h=10 m時采用偶極-偶極裝置沿Z軸旋轉(zhuǎn)的正演結(jié)果視電阻率曲線。從圖8中的4條曲線可以看出，有3條曲線前端位置都出現(xiàn)不同程度的波動，而整體呈現(xiàn)出不對稱曲線形態(tài)。當(dāng)板狀異常體沿Z軸順時針旋轉(zhuǎn)時，雙峰的高低與旋轉(zhuǎn)方向呈現(xiàn)對應(yīng)關(guān)系，隨著順時針旋轉(zhuǎn)角度增大，極值左邊的單峰逐漸減??；且極值異常大小隨著旋轉(zhuǎn)角度的增大，異常幅值也在增大。

通過低阻板狀體在沿Y軸與Z軸旋轉(zhuǎn)的模擬結(jié)果對比可以發(fā)現(xiàn)，兩種模型雙峰的高低與旋轉(zhuǎn)方向都呈現(xiàn)相對應(yīng)關(guān)系，但是在沿Z軸旋轉(zhuǎn)時，需旋轉(zhuǎn)角度較大(>60°)才有明顯變化；兩種旋轉(zhuǎn)都隨著旋轉(zhuǎn)角度的增大，異常幅值也在增大。

圖7 模型2單極-偶極裝置視電阻率曲線Fig.7 Apparent resistivity curve of pole-dipole device in model 2

圖8 模型2偶極-偶極裝置視電阻率曲線Fig.8 Apparent resistivity curve of dipole-dipole device in model 2

3.2 不同旋轉(zhuǎn)變化橫向比較分析

上述分析主要是在每個角度的異常幅值以及形態(tài)變化上進(jìn)行對比分析，但在橫向上差異并不明顯。因此，將視電阻率曲線圖轉(zhuǎn)換成視電阻率斷面圖，進(jìn)一步分析異常體角度變化在橫向上的差異。以單極-偶極裝置視電阻率斷面為例，如圖9所示。

由圖9(a)可知，當(dāng)異常體繞Y軸旋轉(zhuǎn)為0°時，可以準(zhǔn)確反映異常體位置，但異常范圍較?。浑S著旋轉(zhuǎn)角度增加，低阻異常所在位置相對于真實(shí)位置有所提前，且異常范圍開始變大。由圖9(b)可知，當(dāng)異常體繞Z軸旋轉(zhuǎn)，異常范圍變化很小，但隨著旋轉(zhuǎn)角度的增加，異常體的位置有所提前。

因此，在用直流電阻率勘探中，應(yīng)綜合考慮地質(zhì)構(gòu)造的產(chǎn)狀，從而對異常體進(jìn)行校正、歸位。

4 工程實(shí)例

為了進(jìn)一步說明斷層對直流電阻率法結(jié)果的影響，選取野外探測實(shí)例進(jìn)行說明。為保證廣西南寧至玉林城際鐵路某段建設(shè)安全，在該段內(nèi)進(jìn)行斷層等地質(zhì)災(zāi)害的勘測。

該段地表植被發(fā)育，地表為第四系覆蓋，主要巖性有白堊系、寒武系泥巖、砂巖，物性參數(shù)如表2所示。由表2可知，第四系、白堊系砂巖夾泥巖之間存在一定的電性差異，較完整巖體破碎、較破碎或含水巖體之間在電阻率值和等值線形態(tài)有明顯的區(qū)別，因此該工區(qū)具備開展直流電阻率法的地球物理勘探條件。

表2 物性參數(shù)統(tǒng)計表Table 2 Statistical table of physical parameters

在該區(qū)段布置一條#1高密度測線，選用溫納裝置進(jìn)行觀測，電極間距5 m，最大隔離系數(shù)為31，通過逐次擴(kuò)大電極距a值，使探測深度逐漸加深，觀測測點(diǎn)處垂直方向由淺到深的電阻率變化。資料處理采用瑞典高密度二維反演軟件RES2DINV，最后將反演數(shù)據(jù)利用Surfer及 AutoCAD軟件繪制成果圖。

圖9 單極-偶極視電阻率斷面圖Fig.9 Profile of apparent resistivity of pole-dipole device

圖10 視電阻率斷面圖及地質(zhì)斷面圖Fig.10 Apparent resistivity section and geological section

測線視電阻率斷面圖和地質(zhì)斷面圖如圖10所示，由淺及深電阻率可以分3層，淺地表電阻率較為不均勻，推測其主要為第四系、覆蓋層，電阻率較低，為10～30 Ω·m；測點(diǎn)距離0～400 m：中間層為白堊系砂巖夾泥巖，較破碎，電阻率ρ=10～20 Ω·m，下部為白堊系砂巖夾泥巖，較完整，電阻率ρ=20～60 Ω·m；測點(diǎn)距離400～800 m：中間層為寒武系砂巖夾泥巖，較破碎，電阻率ρ=10～100 Ω·m，下部為寒武系砂巖夾泥巖，較完整，電阻率ρ=100～500 Ω·m；測點(diǎn)距離300～400 m，電阻率等值線變陡呈低阻縱向延伸，通過后期鉆探資料得知，由斷裂引起，斷層破碎帶，與基巖傾角約60°，電阻率ρ=10～50 Ω·m。

從圖10(a)與圖10(b)對比中可知，在地質(zhì)斷面圖中斷層的真實(shí)位置大概在測線300～400 m，而圖10(a)中通過高密度電法所測得的斷層范圍顯示在測線400～550 m。因此，在進(jìn)行直流電阻率法測量斷層等異常時，傾角變化會導(dǎo)致其真實(shí)位置發(fā)生偏移，且異常范圍也會增加，在后期應(yīng)進(jìn)行校正處理才能精確定位。

5 結(jié)論

(1)采用COMSOL軟件對三維半空間模型正演模擬，通過驗(yàn)證數(shù)值解與解析解的誤差，可以知道利用COMSOL軟件進(jìn)行三維半空間模型正演計算是可行的。

(2)在模型模擬中，采用單級-偶極裝置和偶極-偶極裝置來模擬低阻板狀體繞Y軸與Z軸的變化，繞Y軸旋轉(zhuǎn)其極值兩側(cè)雙峰的高低與旋轉(zhuǎn)方向有很大關(guān)系，但是在沿Z軸旋轉(zhuǎn)時，需旋轉(zhuǎn)角度較大(>60°)才有明顯變化；兩種旋轉(zhuǎn)都隨著旋轉(zhuǎn)角度的增大，異常幅值也在增大。

(3)兩種旋轉(zhuǎn)方式在橫向上的變化，隨著繞Y軸旋轉(zhuǎn)角度增加，低阻異常所在位置相對于真實(shí)位置有所提前，且異常范圍增大；繞Z軸旋轉(zhuǎn)，異常范圍變化較小，位置有所提前。

(4)實(shí)際數(shù)據(jù)得到的反演結(jié)果與數(shù)值模擬的結(jié)論相吻合，當(dāng)有傾角斷層時，會導(dǎo)致直流電阻率法結(jié)果出現(xiàn)斷層真實(shí)位置偏移，異常范圍增加，所以應(yīng)對后期結(jié)果進(jìn)行校正才能準(zhǔn)確定位。