起伏地形下偶極-偶極激電測(cè)深二維反演軟件開發(fā)及應(yīng)用

顧觀文 ， 吳文鸝 ， 林品榮 ， 梁 萌

(中國地質(zhì)科學(xué)院 地球物理地球化學(xué)勘查研究所， 廊坊 065000)

0 前言

基于陣列方式的激電測(cè)深觀測(cè)，由于工作效率高、觀測(cè)數(shù)據(jù)密度大，在資源勘查、工程和環(huán)境勘察等方面應(yīng)用廣泛。其觀測(cè)數(shù)據(jù)的處理解釋，目前國內(nèi)主要依賴于國外商用高密度電阻率法數(shù)據(jù)處理解釋軟件，國外有代表性的高密度電法二維成像軟件為Res2dinv和Earthimage 2D,Res2dinv為美國Geotomo公司開發(fā)，目前已授權(quán)給瑞典ABEM公司和德國DMT公司，軟件可對(duì)視電阻率和視激電數(shù)據(jù)進(jìn)行二維反演成像，其突出特點(diǎn)是能靈活支持多種觀測(cè)裝置。Earthimage 2D為美國AGI公司與其高密度電法儀配套的視電阻率/視激電數(shù)據(jù)二維反演成像軟件。

作者在前人工作基礎(chǔ)上[1-3]，采用第二類齊邊界條件結(jié)合相應(yīng)網(wǎng)格剖分技術(shù)實(shí)現(xiàn)偶極-偶極激電測(cè)深二維正演計(jì)算，并將其引入偶極—偶極激電測(cè)深二維反演中，同時(shí)在可視化開發(fā)環(huán)境下將二維反演與可視化技術(shù)進(jìn)行有機(jī)集成，形成集數(shù)據(jù)輸入、數(shù)據(jù)可視、反演計(jì)算和結(jié)果輸出為一體的二維反演軟件，通過理論模型合成數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)的反演表明，研制形成的二維反演軟件具備一定的實(shí)用性，并能與物化探所研制的大功率多功能電法儀器配套使用。

1 二維正反演算法

1.1 正演模擬

1.1.1 波數(shù)域點(diǎn)源電位的變分問題[4-6]

三維點(diǎn)源電場(chǎng)的總電位通過沿走向方向的傅立葉變換至波數(shù)域后，可將三維點(diǎn)源直流電場(chǎng)的邊值問題轉(zhuǎn)化為2.5維直流電場(chǎng)的邊值問題(式(1))來進(jìn)行處理。

(1)

此時(shí)波數(shù)域點(diǎn)源二維直流電位U所對(duì)應(yīng)的變分問題為式(2)

(2)

式(1)和式(2)中，Ω為研究區(qū)域；Γs和?！逓檠芯繀^(qū)域的地面邊界和無窮遠(yuǎn)邊界；σ為介質(zhì)電導(dǎo)率；κ為傅立葉變換波數(shù)；r為電源點(diǎn)至Γ∞積分點(diǎn)的距離；K0、K1為零階、一階第二類修正貝塞爾函數(shù)；n是Ω的外法向。

1.1.2 齊次邊界條件的引入及模型剖分[4,7]

為了適應(yīng)簡(jiǎn)化后的邊界條件，在模型網(wǎng)格剖分時(shí)，將計(jì)算區(qū)域劃分為觀測(cè)區(qū)和擴(kuò)展區(qū)，觀測(cè)區(qū)為地面電極(包括供電和測(cè)量電極)布置區(qū)，擴(kuò)展區(qū)是在觀測(cè)區(qū)的基礎(chǔ)上分別沿測(cè)線方向和深度方向擴(kuò)展形成的區(qū)域(圖1)，以近似無窮遠(yuǎn)邊界。觀測(cè)區(qū)沿測(cè)線和深度方向按測(cè)點(diǎn)間距等間隔剖分，在已剖分觀測(cè)區(qū)的基礎(chǔ)上，沿觀測(cè)區(qū)測(cè)線方向左右兩邊及深度方向擴(kuò)展網(wǎng)格數(shù)通常為13，擴(kuò)展區(qū)網(wǎng)格步長按1.3倍遞增[8]。

邊界條件簡(jiǎn)化后，式(2)的變分問題簡(jiǎn)化為：

(3)

用有限單元法求變分問題(3)[5]，歸結(jié)到求波數(shù)域中各節(jié)點(diǎn)上電位的線性方程組

κU=s

(4)

解線性方程組(4)，求得各節(jié)點(diǎn)波數(shù)域中的總電位U。

1.1.3 視電阻率和視極化率的計(jì)算

通過對(duì)波數(shù)域中電位U,由式(5)做傅立葉余弦逆變換求得主剖面上的三維空間中的總電位V[5]。

(5)

由各節(jié)點(diǎn)的總電位V及對(duì)應(yīng)的裝置系數(shù)，可得偶極-偶極測(cè)量裝置的視電阻率表達(dá)式[9]為：

ρs(A,B,M,N)=K(A,B,M,N)[V(A,M)-V(A,N)-V(B,M)+V(B,N)]

(6)

式(6)中：K(A,B,M,N)為偶極裝置的裝置系數(shù)；V(A,M)、V(A,N)、V(B,M)和V(B,N)分別為A點(diǎn)供電M、N處電位和B點(diǎn)供電M、N處電位，供電電流為常數(shù)，一般假設(shè)為1 A。計(jì)算出各節(jié)點(diǎn)的電位后，即可實(shí)現(xiàn)偶極-偶極裝置視電阻率的計(jì)算，采用等效電阻率法求取視極化率[9]。

1.2 視電阻率和視極化率數(shù)據(jù)的反演 [1]

采用帶先驗(yàn)信息的最小二乘反演方法，目標(biāo)函數(shù)為

φ=‖Wd(Δd-AΔm)‖2+‖Wm(m-mb+Δm)‖2

(7)

式(7)中 ：m為預(yù)測(cè)模型參數(shù)矢量；mb為基本模型參數(shù)矢量；Δd為實(shí)測(cè)視電阻率與正演視電阻率的數(shù)據(jù)差矢量；A為偏導(dǎo)數(shù)矩陣；m、mb、Δd和A為取對(duì)數(shù)后的值。Wd=diag(1/σ1，1/σ2，…，1/σN)為數(shù)據(jù)的擬方差矩陣；σi為第i個(gè)數(shù)據(jù)的均方誤差；Wm是模型加權(quán)矩陣，即模型先驗(yàn)信息。式(7)對(duì)Δm求導(dǎo)并令其等于零，可得到線性方程組(8)

(8)

圖1 正演計(jì)算區(qū)域網(wǎng)格剖分Fig．1 Subdividing of region for forward modeling

解方程組(8)得到的模型修改量Δm,將Δm與m相加，得新的預(yù)測(cè)模型參數(shù)矢量，繼續(xù)下一次迭代，直至平均均方誤差滿足要求。平均均方誤差rms定義為

(9)

按照Seigel(1959)理論，如果地下空間由M塊不同電阻率ρj和本征極化率ηj的巖礦組成(j=1,2,…,M)，視極化率響應(yīng)可表示為

或?qū)憺榫仃囆问?/p>

ηa=Aη

(10)

式中ηa是視極化率響應(yīng)矢量；η為本征極化率參數(shù)矢量；A是偏導(dǎo)數(shù)矩陣。最小二乘反演方法的目標(biāo)函數(shù)為

φ=‖Wd(ηa-Aη)‖2+‖Wη(η-ηb)‖2

(11)

極化率模型增量求解與電阻率模型增量的求解過程相同，偏導(dǎo)數(shù)陣A已在視電阻率反演過程中求得，因此視極化率的反演計(jì)算量減小。

2 軟件實(shí)現(xiàn)及模型檢驗(yàn)

軟件運(yùn)行環(huán)境：Windows XP/Win7。軟件開發(fā)工具：MS VC6.0,Compaq Fortran6.6。

2.1 軟件的功能組成及集成實(shí)現(xiàn)

二維反演軟件系統(tǒng)主要包括數(shù)據(jù)輸入、數(shù)據(jù)顯示、反演計(jì)算和反演結(jié)果的輸出等功能(圖2)。

軟件集成環(huán)境基于VC6.0，數(shù)據(jù)輸入、數(shù)據(jù)顯示和結(jié)果輸出模塊采用VC6.0編寫，可在源代碼級(jí)別上基于類對(duì)象方式由集成程序編譯使用。二維反演計(jì)算模塊由FORTRAN語言編寫，使用Compaq Fortran編譯形成動(dòng)態(tài)庫，提供接口由程序主進(jìn)程調(diào)用。

集成后形成的軟件界面如圖3所示。界面中左邊部分顯示信息包括電極個(gè)數(shù)、點(diǎn)距、數(shù)據(jù)點(diǎn)數(shù)、最大最小隔離系數(shù)和反演迭代誤差等信息，右邊部分上部為實(shí)測(cè)數(shù)據(jù)擬斷面，中部為最后一次迭代反演的正演數(shù)據(jù)擬斷面，下部為反演模型斷面，通過界面上迭代誤差大小、實(shí)測(cè)和正演擬斷面的相似程度可衡量反演結(jié)果的有效程度。

圖2 軟件功能組成Fig．2 Composition of software function

圖3 二維反演軟件界面Fig．3 Interface of 2D inversion software

2.2 模型檢驗(yàn)

對(duì)施加5% 隨機(jī)噪聲的兩個(gè)地電模型(模型來源于“2005年全國電法及電磁法勘探正反演軟件推優(yōu)會(huì)”)的合成數(shù)據(jù)進(jìn)行二維反演試算，以檢驗(yàn)二維反演程序的有效性。

2.2.1 模型一

裝置參數(shù)：80根電極，點(diǎn)距為5 m，最小隔離系數(shù)為“1”，最大隔離系數(shù)為20。

模型描述：在均勻圍巖中存在一個(gè)低阻高極化異常體。圍巖電阻率為100 Ω·m，極化率為“1”。異常體電阻率為10 Ω·m，極化率為10。異常體的寬度和向下延伸長度分別為20 m和10 m，上頂距地表深度為10 m。異常體從測(cè)線坐標(biāo)190 m開始，延續(xù)到坐標(biāo)210 m。模型示意圖見圖4。

反演結(jié)果：圖5和圖6為迭代5次的反演結(jié)果，RMS(均方誤差)為6%，從反演斷面圖可以看出低阻高極化異常體得到了很好地反映，反演模型斷面和真實(shí)理論模型一致。

2.2.2 模型二

裝置參數(shù)：60根電極，點(diǎn)距為5 m，最小隔離系數(shù)為“1”，最大隔離系數(shù)為20。

模型描述：在均勻圍巖中分別存在一個(gè)低阻高極化體和一個(gè)高阻高極化體。圍巖電阻率為10 Ω·m,極化率為10。異常體1電阻率為2 Ω·m，極化率為50，異常體2電阻率為 50 Ω·m，極化率為50。異常體1和異常體2的寬度和向下延伸長度均為15 m和10 m,上頂距地表的深度為10 m。橫向上第一個(gè)異常體從100 m開始延續(xù)到115 m，第二個(gè)異常體從210m開始延續(xù)到225 m。模型示意圖見圖7。

反演結(jié)果：圖8和圖9為迭代5次的反演結(jié)果，RMS(均方誤差)為7%，從圖中可以看出，兩個(gè)異常體在橫向和縱向上都得到了很好的歸位，反演斷面和真實(shí)模型斷面(圖7)一致。

3 應(yīng)用實(shí)例

作者采用本文編制的偶極—偶極激電測(cè)深二維反演軟件，對(duì)內(nèi)蒙某工區(qū)的500線相位激電資料進(jìn)行二維反演。

圖4 模型一示意圖Fig．4 Sketch of model 1

圖5 反演電阻率斷面圖Fig．5 Resistivity section of inversion model

圖6 反演極化率斷面圖Fig．6 IP section of inversion model

圖7 模型二示意圖Fig．7 Sketch of model 2

圖8 反演電阻率斷面圖Fig．8 Resistivity section of inversion model

圖9 反演極化率斷面圖Fig．9 IP section of inversion model

在工作區(qū)的500線上，采用軸向偶極—偶極裝置，開展了頻率域激電的測(cè)深測(cè)量，觀測(cè)參量為視電阻率ρs和絕對(duì)相位φs。工作頻率4 Hz、2 Hz、1 Hz、2 S、4 S。野外工作裝置如圖10所示。供電極距AB=80 m、接收極距MN=80 m，點(diǎn)距為40 m，隔離系數(shù)n=1～12。

采用二維反演軟件對(duì)4 s時(shí)的視電阻率和視相位進(jìn)行反演，反演初始模型為均勻半空間，均勻半空間的物性值分別為記錄點(diǎn)的平均視電阻率和平均視相位。圖11為迭代5次的反演結(jié)果，RMS(均方誤差)為13%，從圖11、圖12、圖13和圖14可以看出，實(shí)測(cè)和正演視電阻率擬斷面以及實(shí)測(cè)和正演視相位擬斷面形態(tài)相似程度高。

圖15為500線相位和電阻率反演斷面(由藍(lán)色漸變至紅色對(duì)應(yīng)低值漸變至高值)及中梯激電剖面曲線，從圖15中可以看出，礦致異常在相位和電阻率反演斷面有明顯反映，激電反演異常位置與中梯剖面曲線激電異常位置對(duì)應(yīng)一致。

圖10 相位激電軸向偶極-偶極裝置示意圖Fig．10 Sketch of dipole-dipole array IP

圖11 實(shí)測(cè)視電阻率擬斷面Fig．11 Pseudosection of observed apparent resistivity

圖12 第五次迭代正演視電阻率擬斷面Fig．12 Forward apparent resistivity pseudosection of the fifth iteration

圖13 實(shí)測(cè)視相位擬斷面Fig．13 Pseudosection of observed apparent phase

圖14 第五次迭代正演視相位擬斷面Fig．14 Forward apparent phase pseudosection of the fifth iteration

圖15 500測(cè)線上獲取的礦致異常Fig．15 Mineralization anomalies of the line 500(a)激電中梯測(cè)量剖面；(b)500 線激電深相位反演斷面；(c)500 線激電深電阻率反演斷面

4 結(jié)論

1)在有限元點(diǎn)源二維正演模擬中，采用第二類齊次邊界條件結(jié)合相應(yīng)的剖分技術(shù)。將二維正演模擬引入二維反演中，為偶極—偶極激電測(cè)深二維反演軟件的實(shí)現(xiàn)奠定基礎(chǔ)。

2)利用VC 6.0和Compaq Fortran6.6開發(fā)工具，將二維正反演算法與可視化編程技術(shù)有機(jī)集成，形成了集數(shù)據(jù)輸入、數(shù)據(jù)成圖、反演計(jì)算和結(jié)果多方式輸出為一體的偶極—偶極激電測(cè)深二維反演軟件。

3)通過理論模型合成數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)的反演表明，研制形成的二維反演軟件對(duì)時(shí)間域激發(fā)極化數(shù)據(jù)和相位激電數(shù)據(jù)的反演是有效的，并能與物化探所研制的大功率多功能電法儀器配套使用。

參考文獻(xiàn)：

[1] 阮百堯，村上裕，徐世浙. 激發(fā)極化數(shù)據(jù)的最小二乘二維反演方法[J]. 地球科學(xué)—中國地質(zhì)大學(xué)學(xué)報(bào)，1999，24(6)：620-623.

[2] 吳文鸝. 電法勘探工作站軟件系統(tǒng)簡(jiǎn)介[J]. 地質(zhì)與勘探，2003(增刊)：147-149.

[3] 顧觀文. 電法勘探工作站程序設(shè)計(jì)技術(shù)[J]. 地質(zhì)與勘探，2003(增刊)：152-154.

[4] 顧觀文，吳文鸝，高艷芳，等. 電阻率/激電測(cè)深二維人機(jī)交互正演模擬[J]. 物探化探計(jì)算技術(shù)，2007,29：89-92.

[5] 徐世浙. 地球物理中的有限單元法[M]. 北京：科學(xué)出版社，1994.

[6] 阮百堯. 三角單元部分電導(dǎo)率分塊連續(xù)變化點(diǎn)源二維電場(chǎng)有限元數(shù)值模擬[J].廣西科學(xué)，2001，8(1)：1-5.

[7] 黃俊革.三維電阻率/極化率有限元正演模擬與反演成像[D].長沙:中南大學(xué)，2003.

[8] 羅延鐘，萬樂，董浩斌，等.不平地形條件下高密度電阻率法的2.5維反演[J].地質(zhì)與勘探，2004(增刊)：172-175.

[9] 李金銘.地電場(chǎng)與電法勘探[M].北京：地質(zhì)出版社，2005.