肖 波，劉海飛，戴前偉

(1.廣東省電力設(shè)計(jì)研究院，廣州 510600；2.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙 410083)

0 前言

固定點(diǎn)源激電測(cè)深法(FPS)又稱(chēng)固定點(diǎn)源雙邊三極測(cè)深法(圖1)。由前蘇聯(lián)柯馬羅夫B.A[1]于60年代提出。工作中供電正極A固定不動(dòng)，供電負(fù)極B置于無(wú)窮遠(yuǎn)處，測(cè)量電極MN在A(yíng)極兩邊移動(dòng)。該方法測(cè)量極距可根據(jù)實(shí)際情況靈活變化，導(dǎo)體極化場(chǎng)位置隨著供電點(diǎn)位置變化而變化，且在工作中能采用多臺(tái)接收機(jī)同時(shí)測(cè)量，在具有較高工作效率的同時(shí)可更換確定極化體的空間位置。但該裝置在數(shù)據(jù)處理過(guò)程中，常規(guī)的相對(duì)強(qiáng)度-圓弧交匯定量解釋法對(duì)于簡(jiǎn)單地電條件下的模型解釋通常能取得較好效果，但在起伏地電條件下，其解釋效果通常會(huì)受到較為明顯的地形影響，從而未能精確確定極化體的中心深度及上界面位置，為勘探解釋工作帶來(lái)困難。

圖1 固定點(diǎn)源測(cè)深法裝置簡(jiǎn)圖Fig．1 The fixed point source IP sounding equipment

對(duì)此從有限元正演模擬和最小二乘反演方法的角度入手，對(duì)起伏地形條件下的地電模型進(jìn)行實(shí)例計(jì)算，并將最小二乘反演結(jié)果與圓弧交匯—相對(duì)強(qiáng)度計(jì)算結(jié)果相對(duì)比。由對(duì)比結(jié)果表明：最小二乘反演方法能較好地適應(yīng)起伏地形條件進(jìn)行反演計(jì)算，其反演效果優(yōu)于相對(duì)強(qiáng)度-圓弧交匯定量解釋法，能更有效地體現(xiàn)出實(shí)際異常體的真實(shí)信息。

1 正演模擬

正演模擬采用有限單元法進(jìn)行計(jì)算。將二維地電斷面區(qū)域離散成許多相互連接的網(wǎng)格單元，對(duì)各個(gè)網(wǎng)格單元節(jié)點(diǎn)電位進(jìn)行求解得到地電斷面的電位分布，通過(guò)將電位分布轉(zhuǎn)化為視電阻率從而得到斷面各個(gè)節(jié)點(diǎn)的電阻率值，根據(jù)“等效電阻率”公式求取出視極化率值。

在二維地電條件下，點(diǎn)源場(chǎng)中各網(wǎng)格節(jié)點(diǎn)的電位計(jì)算可歸結(jié)為對(duì)若干個(gè)給定波數(shù)λ求解電位的傅氏變換v(x,λ,z)所滿(mǎn)足的二維偏微分方程的邊值問(wèn)題：

(1)

(2)

采用有限單元法對(duì)變分問(wèn)題求解，得出變換電位V(x,λ,z)，作傅里葉逆變換：

將電位分布值轉(zhuǎn)化為視電阻率值，通過(guò)“等效電阻率”求取，從而得到視極化率值。

2 解釋方法

2.1 常規(guī)定量解釋方法

在電法勘探中，通過(guò)以?xún)蓚€(gè)以上的供電點(diǎn)到ρs/ηs曲線(xiàn)極值點(diǎn)P的距離為半徑，分別在二維地電斷面上做圓弧，圓弧的交點(diǎn)位置近似為異常體的空間中心位置，該作圖方法稱(chēng)為圓弧交匯法[3-4]，如圖2所示。

圖2 圓弧交匯法示意圖Fig．2 Diagram of arc intersection method

相對(duì)強(qiáng)度法是以圓弧交匯法為基礎(chǔ)，將二維地電斷面網(wǎng)格化并以點(diǎn)源為圓心，以點(diǎn)源到各網(wǎng)格節(jié)點(diǎn)的距離為半徑作圓弧，取圓弧與水平地表交點(diǎn)位置對(duì)應(yīng)的視電阻率/視極化率在地表的投影值為該網(wǎng)格節(jié)點(diǎn)的值，位于點(diǎn)源坐標(biāo)左、右邊的網(wǎng)格節(jié)點(diǎn)值分別由視電阻率/視極化率左、右支實(shí)測(cè)值確定。對(duì)于多個(gè)點(diǎn)源，先對(duì)各節(jié)點(diǎn)上 求和進(jìn)而計(jì)算出他們的平均值 找出斷面上最大的平均值 ，用各網(wǎng)格節(jié)點(diǎn)的平均值 去除以 ，從而得到了相對(duì)強(qiáng)度值 ；對(duì)于低阻體則相反，提取斷面上最小平均值 并用該值除以網(wǎng)格各節(jié)點(diǎn)的平均值 ，從而達(dá)到突出斷面異常體的中心位置的目的。相對(duì)強(qiáng)度法的作圖如圖3所示。

圖3 相對(duì)強(qiáng)度法的成圖示意圖Fig．3 Diagram of relative intensity method

2.2 最小二乘反演

反演采用基于圓滑約束的最小二乘反演方法，利用實(shí)測(cè)數(shù)據(jù)和正演模型構(gòu)造一目標(biāo)函數(shù)，并使其達(dá)到極小。圓滑約束最小二乘是基于以下方程[5～6]：

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

(3)

式中 Δm=m-m0：模型參數(shù)修改矢量(m為N維模型參數(shù)矢量，m0為初始模型參數(shù)矢量)；mb為基本模型參數(shù)矢量；Δd=d-d0：數(shù)據(jù)殘差矢量(d為M維實(shí)測(cè)視電阻率向量，d0為初始模型的觀(guān)測(cè)值)；J為偏導(dǎo)數(shù)矩陣;λ為拉格朗日乘數(shù)；C為光滑度矩陣。

對(duì)方程組(3)求解從而得到模型參數(shù)修正矢量Δm，并將其代入：

m(k)=m(k-1)+Δm

(4)

從而得到預(yù)測(cè)模型參數(shù)矢量m(k)。如此重復(fù)，直至實(shí)測(cè)數(shù)據(jù)與模型數(shù)據(jù)之間的平均均方誤差滿(mǎn)足要求為止，電阻率反演過(guò)程結(jié)束[7-9]。平均均方誤差的計(jì)算公式為：

(5)

3 模型計(jì)算

采用有限單元法進(jìn)行正演計(jì)算，分別用圓弧交匯—相對(duì)強(qiáng)度法、最小二乘反演進(jìn)行計(jì)算。將多個(gè)模型實(shí)例進(jìn)行模型計(jì)算，并將兩種不同計(jì)算方法的進(jìn)行對(duì)比分析。

綜合考慮數(shù)據(jù)量、反演效果及工作量之間的關(guān)系，擬布設(shè)7個(gè)測(cè)深供電點(diǎn)：A1=40 m、A2=60 m、A3=80 m、A4=100 m、A5=120 m、A6=140 m、A7=160 m。為取得有效勘探精度與深度，模擬測(cè)線(xiàn)為-50 m～250 m，實(shí)取測(cè)線(xiàn)0 m～200 m數(shù)據(jù)進(jìn)行剖面成圖。計(jì)算結(jié)果剖面圖中的藍(lán)色方框?yàn)槟Ｐ屯队皥D(圖4)。

圖4 水平地形下正方柱體模型示意圖Fig．4 Schematic diagram of affirmative column model under horizontal topography

3.1 水平地形

模型為一個(gè)20 m×20 m低阻、高極化正方柱體(見(jiàn)圖4)；中心位置(100 m,-20 m)，電阻率ρ1=50 Ω·m，極化率η1=5.0%；背景電阻率ρ0=100 Ω·m，極化率：η0=1.0%。正演模擬計(jì)算得到的各點(diǎn)源視電阻率、視極化率曲線(xiàn)如圖5、圖6所示。視電阻率、視極化率相對(duì)強(qiáng)度圖如圖7、圖8所示。視電阻率、視極化率最小二乘反演結(jié)果如圖9、圖10所示。

圖5 水平地形下正方柱體模型正演模擬各點(diǎn)源測(cè)深視電阻率曲線(xiàn)圖Fig．5 The apparent resistivity curve of affirmative column model under horizontal topography

圖6 水平地形下正方柱體模型正演模擬各點(diǎn)源測(cè)深視極化率曲線(xiàn)圖Fig．6 The polarizability curve of affirmative column model under horizontal topography

圖7 正方柱體模型正演視電阻率相對(duì)強(qiáng)度圖Fig．7 The arc intersection - relative intensity result of apparent resistivity about affirmative column model under horizontal topography

圖8 水平地形下正方柱體模型正演視極化率相對(duì)強(qiáng)度圖Fig．8 The arc intersection - relative intensity result of polarizability about affirmative column model under horizontal topography

圖9 水平地形下正方柱體模型視電阻率最小二乘反演結(jié)果Fig．9 The least square inversion result of apparent resistivity about affirmative column model under horizontal topography

圖10 水平地形下正方柱體模型視極化率最小二乘反演結(jié)果Fig．10 The least square inversion result of pola-rizability about affirmative column model under horizontal topography

由圖7、圖8可知：圓弧交匯—相對(duì)強(qiáng)度法計(jì)算結(jié)果能很好地反映出正方柱體的中心埋深。該異常體中心位置與模型中心位置吻合關(guān)系好，但該異常相對(duì)強(qiáng)度等值線(xiàn)圖不能體現(xiàn)模型的基本形態(tài)，且異常范圍較實(shí)際模型大。

對(duì)比發(fā)現(xiàn)：圖9、圖10表明最小二乘反演能較好地反映模型的基本形態(tài)和基本位置，其異常范圍與模型大小基本一致，但最小二乘反演結(jié)果未能精確反映出模型極化體中心，反演結(jié)果中極化體中心位于實(shí)體模型中上部位置。

3.2 山谷地形

模型為一個(gè)低阻、高極化長(zhǎng)方體模型(圖11)。模型中心位置，電阻率為，極化率為η1=5.0%，模型中心頂面與地表相距；背景電阻率，背景極化率η0=1.0%。正演模擬計(jì)算得到的各點(diǎn)源視電阻率、視極化率曲線(xiàn)如圖12、圖13所示。修正后的圓弧交匯—相對(duì)強(qiáng)度法[7]計(jì)算結(jié)果如圖14、圖15所示。視電阻率、視極化率最小二乘反演結(jié)果如圖16、圖17所示。

圖11 山谷地形下長(zhǎng)方體模型示意圖Fig．11 Schematic diagram of cuboid model under valley topography

圖12 山谷地形下長(zhǎng)方體模型正演模擬各點(diǎn)源測(cè)深視電阻率曲線(xiàn)圖Fig．12 The apparent resistivity curve of cuboid model under valley topography

圖13 山谷地形下長(zhǎng)方體模型正演模擬各點(diǎn)源測(cè)深視極化率曲線(xiàn)圖Fig．13 The polarizability curve of cuboid model under valley topography

圖14 山谷地形下長(zhǎng)方體模型正演視電阻率相對(duì)強(qiáng)度圖Fig．14 The arc intersection - relative intensity result of apparent resistivity about cuboid model under valley topography

由圖14、圖15可知：圓弧交匯—相對(duì)強(qiáng)度法計(jì)算結(jié)果能確定模型的基本位置，但異常體的相對(duì)強(qiáng)度中心與模型實(shí)際位置不相吻合，異常體中心位置位于實(shí)際模型底界面處，同時(shí)異常體的基本形態(tài)、大小與模型實(shí)際存在較大差距。其原因?yàn)橛捎谏焦鹊匦斡绊?，模型的異常范圍變窄?dǎo)致異常體極大值點(diǎn)與其對(duì)應(yīng)的電源點(diǎn)距離縮短；對(duì)于山脊而言剛好相反，模型的異常體范圍變大導(dǎo)致異常體極大值點(diǎn)與其對(duì)應(yīng)的電源點(diǎn)距離變大。盡管再做圓弧交匯-相對(duì)強(qiáng)度進(jìn)行成圖計(jì)算時(shí)已進(jìn)行地形校正處理，但地形影響在一定程度上仍然存在，當(dāng)極化體埋深逐漸增大時(shí)，該影響隨之減弱。

對(duì)比發(fā)現(xiàn)：圖16、圖17表明最小二乘反演能更好地反映出模型的基本形態(tài)及空間位置，異常體的基本輪廓及中心埋深位置與實(shí)際模型吻合關(guān)系一致。

3.3 山脊地形

模型為一個(gè) 低阻、高極化長(zhǎng)方體模型(圖18)，模型中心位置，電阻率 ，極化率η1=5.0%，模型中心頂面與地表相距5 m；背景電阻率 ，背景極化率η0=1.0%。正演模擬計(jì)算得到的各點(diǎn)源視電阻率、視極化率曲線(xiàn)如圖19、圖20所示。圓弧交匯-相對(duì)強(qiáng)度法計(jì)算結(jié)果如圖21、圖22所示。視電阻率、視極化率最小二乘反演結(jié)果如圖23、圖24所示。

由圖21、圖22可知：圓弧交匯-相對(duì)強(qiáng)度法計(jì)算結(jié)果能確定模型的基本位置，但異常體的相對(duì)強(qiáng)度中心與模型實(shí)際位置不相吻合，異常體中心位置位于實(shí)際模型的底界面位置。此外不能根據(jù)異常體準(zhǔn)確推斷出模型的大小及產(chǎn)狀。

圖24 山脊地形下長(zhǎng)方體模型視極化率最小二乘反演結(jié)果Fig．24 The least square inversion result of polari-zability about cuboid model under ridge topography

相比較而言：圖23、圖24表明最小二乘反演能較好地反映模型的基本形態(tài)及空間位置。反演結(jié)果中異常體中心位置及輪廓與實(shí)際模型中心位置及邊界吻合關(guān)系好，反演結(jié)果表明，最小二乘法反演能較好地適應(yīng)起伏地形進(jìn)行計(jì)算。

4 結(jié)論

(1)在水平地形條件下，圓弧交匯-相對(duì)強(qiáng)度法與最小二乘反演都能較好地反映出模型異常體的空間位置，圓弧交匯-相對(duì)強(qiáng)度法在確定異常體中心位置甚至優(yōu)于最小二乘反演計(jì)算結(jié)果。

(2)在起伏地形條件下，圓弧交匯-相對(duì)強(qiáng)度法與最小二乘反演都能較準(zhǔn)確地反映出異常體中心位置，兩種解釋方法直觀(guān)且其解釋效果與實(shí)際較為符合。

(3)采用圓弧交匯-相對(duì)強(qiáng)度法進(jìn)行作圖解釋?zhuān)涞匦斡绊懺谝欢ǔ潭壬先匀淮嬖凇４送怆y以根據(jù)異常體輪廓準(zhǔn)確推斷出模型的大小及產(chǎn)狀。

(4)最小二乘反演結(jié)果能較好地反映出模型的基本形態(tài)及空間位置。反演結(jié)果剖面圖中異常體中心位置及其輪廓模型中心位置吻合及模型邊界吻合關(guān)系好，表明最小二乘法反演能很好地適應(yīng)起伏地形進(jìn)行反演計(jì)算，其反演結(jié)果優(yōu)于圓弧交匯-相對(duì)強(qiáng)度法。

參考文獻(xiàn)：

[1] KOMAPOBB А.Теоретическиеосновыинтерпретацийн-аблюденийВметодевызваннойполяризаций[M],Л.Недра,1966.

[2] 戴前偉,肖波,馮德山，等. 基于二維高密度電阻率勘探數(shù)據(jù)的三維反演及應(yīng)用[J].中南大學(xué)學(xué)報(bào)：自然科學(xué)版,2012, 43(1): 0293 -0310.

[3] 劉國(guó)興.電法勘探原理與方法[M].北京.地質(zhì)出版社.2005.

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

[5] 李金銘,魏文博,陳本池,等.固定點(diǎn)源測(cè)深法定量解釋研究[J].物探與化探,1997,21(3):187-196.

[6] 強(qiáng)建科,羅延鐘,熊彬.固定點(diǎn)源測(cè)深激電異常研究[J].地球物理學(xué)進(jìn)展,2005,20(4):1176-1183.

[7] 熊彬,阮百堯,黃俊革. 直流電阻率測(cè)深中二維反演程序?qū)θS數(shù)據(jù)的近似解釋[J].地球科學(xué)-中國(guó)地質(zhì)大學(xué)學(xué)報(bào), 2003,28(1):0102-0106.

[8] 劉海飛.高密度數(shù)據(jù)處理方法研究[D].長(zhǎng)沙:中南大學(xué)信息物理工程學(xué)院,2003.

[9] 張大海,王興泰.二維視電阻率斷面的快速最小二乘反演[J].物探化探計(jì)算技術(shù), 1999, 21(1)：2-8.