起伏地形下偶極-偶極激電測深二維反演軟件開發及應用

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

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

0 前言

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

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

1 二維正反演算法

1.1 正演模擬

1.1.1 波數域點源電位的變分問題[4-6]

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

(1)

此時波數域點源二維直流電位U所對應的變分問題為式(2)

(2)

式(1)和式(2)中，Ω為研究區域；Γs和Γ∞為研究區域的地面邊界和無窮遠邊界；σ為介質電導率；κ為傅立葉變換波數；r為電源點至Γ∞積分點的距離；K0、K1為零階、一階第二類修正貝塞爾函數；n是Ω的外法向。

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

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

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

(3)

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

κU=s

(4)

解線性方程組(4)，求得各節點波數域中的總電位U。

1.1.3 視電阻率和視極化率的計算

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

(5)

由各節點的總電位V及對應的裝置系數，可得偶極-偶極測量裝置的視電阻率表達式[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)為偶極裝置的裝置系數；V(A,M)、V(A,N)、V(B,M)和V(B,N)分別為A點供電M、N處電位和B點供電M、N處電位，供電電流為常數，一般假設為1 A。計算出各節點的電位后，即可實現偶極-偶極裝置視電阻率的計算，采用等效電阻率法求取視極化率[9]。

1.2 視電阻率和視極化率數據的反演 [1]

采用帶先驗信息的最小二乘反演方法，目標函數為

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

(7)

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

(8)

圖1 正演計算區域網格剖分Fig．1 Subdividing of region for forward modeling

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

(9)

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

或寫為矩陣形式

ηa=Aη

(10)

式中ηa是視極化率響應矢量；η為本征極化率參數矢量；A是偏導數矩陣。最小二乘反演方法的目標函數為

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

(11)

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

2 軟件實現及模型檢驗

軟件運行環境：Windows XP/Win7。軟件開發工具：MS VC6.0,Compaq Fortran6.6。

2.1 軟件的功能組成及集成實現

二維反演軟件系統主要包括數據輸入、數據顯示、反演計算和反演結果的輸出等功能(圖2)。

軟件集成環境基于VC6.0，數據輸入、數據顯示和結果輸出模塊采用VC6.0編寫，可在源代碼級別上基于類對象方式由集成程序編譯使用。二維反演計算模塊由FORTRAN語言編寫，使用Compaq Fortran編譯形成動態庫，提供接口由程序主進程調用。

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

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

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

2.2 模型檢驗

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

2.2.1 模型一

裝置參數：80根電極，點距為5 m，最小隔離系數為“1”，最大隔離系數為20。

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

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

2.2.2 模型二

裝置參數：60根電極，點距為5 m，最小隔離系數為“1”，最大隔離系數為20。

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

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

3 應用實例

作者采用本文編制的偶極—偶極激電測深二維反演軟件，對內蒙某工區的500線相位激電資料進行二維反演。

圖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

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

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

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

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

圖11 實測視電阻率擬斷面Fig．11 Pseudosection of observed apparent resistivity

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

圖13 實測視相位擬斷面Fig．13 Pseudosection of observed apparent phase

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

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

4 結論

1)在有限元點源二維正演模擬中，采用第二類齊次邊界條件結合相應的剖分技術。將二維正演模擬引入二維反演中，為偶極—偶極激電測深二維反演軟件的實現奠定基礎。

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

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

參考文獻：

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

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

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

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

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

[6] 阮百堯. 三角單元部分電導率分塊連續變化點源二維電場有限元數值模擬[J].廣西科學，2001，8(1)：1-5.

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

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

[9] 李金銘.地電場與電法勘探[M].北京：地質出版社，2005.