肖波，劉海飛，戴前偉

?

FPS異常的2種不同解釋方法對比

肖波1，劉海飛2，戴前偉2

(1.中國能源建設集團 廣東省電力設計研究院，廣東 廣州，510600；2.中南大學 地球科學與信息物理學院，湖南 長沙，410083)

從有限單元法正演模擬入手，采用最小二乘法進行反演計算，通過VC++程序實現(xiàn)，對多個地電模型實例進行計算，并將最小二乘法反演結果與圓弧交匯?相對強度計算結果進行對比。研究結果表明：圓弧交匯?相對強度法對簡單模型適用性高，由其反演成果所確定異常體中心位置甚至會優(yōu)于最小二乘法最優(yōu)化反演結果；最小二乘法能適應各種地電條件對多種模型進行反演計算，其計算結果優(yōu)于圓弧交匯?相對強度法計算結果。

電阻率；極化率；數(shù)值模擬；有限元法；反演；定量解釋

固定點源激電測深法(FPS)又稱固定點源雙邊三極測深法。該方法裝置在實際工作中采用多臺接受機同時測量，測量電極極距根據(jù)實際情況可靈活變化，導體極化場位置隨著供電點位置不同而變化，具有較高工作效率的同時可更好地確定異常體的位置。但該裝置在數(shù)據(jù)處理過程中，常規(guī)的相對強度?圓弧交匯解釋方法受場源的分布位置影響大，其解釋成果對模型形狀、埋深半徑比的依賴性較大；該方法對于單個模型的處理解釋通常能取得較好效果，但對組合模型的解釋通常會產(chǎn)生明顯的假異常，從而掩蓋了真實異常體的信息，為勘探解釋工作帶來了極大的困難。在此，本文作者從有限元正演模擬和最小二乘反演方法的角度入手，對多個地電模型進行實例計算，并將最小二乘反演結果與圓弧交匯?相對強度計算結果進行對比。

1 正演模擬

采用有限單元法進行正演模擬，將整個計算空間區(qū)域離散成許多相互連接的單元網(wǎng)格，通過對單元網(wǎng)格節(jié)點的電位進行求解得到計算域的電位分布，然后將電位分布值轉化為視電阻率，再通過“等效電阻率”求取視極化率。

在二維地電條件下，點電流源場各節(jié)點電位的計算可歸結為對若干個給定波數(shù)求解電位的傅里葉變換(，，)所滿足的二維偏微分方程的邊值問題：

式中：為電導率；和1為地面邊界；∞為其他邊界；I為第個供電電極的電流；為場源到測點的距離；邊界處法線向量；和分別為零階和一階修正貝塞爾函數(shù)[1?2]。

與方程式(1)等價的變分問題為

通過有限單元法數(shù)值模擬可得地表任意一點的電位，從而得到相應的視電阻率[3]：

其中：為裝置系數(shù)；為供電電流；Δ為兩節(jié)點之間的電勢差。

當?shù)叵陆橘|有極化特性時，其等效電阻率為

式中：為介質電阻率；為介質極化率。利用上面所述的電阻率二維正演方法對等效電阻率進行正演，可以得到等效視電阻率：

由式(5)可得到視極化率η為

2 解釋方法

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

電法勘探中，有2個或者2個以上的供電點，分別以供電點到視電阻率s曲線、視極化率s曲線極值點的距離為半徑畫圓弧，根據(jù)圓弧的交點確定異常體的中心位置，稱為圓弧交匯法[3?4]，如圖1所示。

圖1 圓弧交匯法示意圖

相對強度法是將二維地電斷面網(wǎng)格化，以各網(wǎng)格節(jié)點到點源的距離為半徑作圓弧，取圓弧與水平地表交點位置對應的視電阻率、視極化率為該網(wǎng)格節(jié)點的值，位于點源坐標左邊的網(wǎng)格節(jié)點由視電阻率、視極化率左支實測值確定，位于點源坐標右邊的網(wǎng)格節(jié)點由視電阻率、視極化率右支實測值確定。對于多個點源，先對各節(jié)點上求和，進而計算出他們的平均值找出斷面上最大的平均值，用各網(wǎng)格節(jié)點的平均值去除以，最后得到了相對強度；對于低阻體，提取斷面上最小平均值，然后用除以各網(wǎng)格節(jié)點的平均值，從而突出斷面異常體的中心位置。相對強度法的作圖示意圖如圖2所示。

圖2 相對強度法的成圖示意圖

2.2 最小二乘反演法

采用基于圓滑約束最小二乘反演方法，利用正演模型和實測數(shù)據(jù)構造一目標函數(shù)，并使其達到極小。圓滑約束最小二乘法是基于以下方程[5?8]：

對方程組(7)進行求解，得到模型參數(shù)修正矢量Δ，將其代入式(8)便得到新的預測模型參數(shù)矢量(k)：

如此重復，直至模擬數(shù)據(jù)和實測數(shù)據(jù)之間的平均均方誤差滿足要求為止[9?14]。平均均方誤差為

在反演過程中，為提高求解大型、病態(tài)線性方程組的計算速度和反演效率，采用廣義共軛梯度算法(GCG)[15?16]對偏導數(shù)矩陣進行求解，其基本公式為

其實現(xiàn)步驟如下。

1) 初始向量(0)=0，(0)=，(0)=T(0)=(0)；誤差限，阻尼因子；

如果||(0)||≤，則終止，否則轉入步驟3)；對于=2，3，…，max，計算到步驟10)；

8) 如果||(j+1)||≤終止, 否則轉入步驟9) ；

3 模型計算

采用有限單元法數(shù)值計算進行正演模擬，分別用圓弧交匯?相對強度法成圖反演和最小二乘最優(yōu)化反演進行反演計算。對多個模型實例的計算結果進行對比分析。

綜合考慮數(shù)據(jù)量、反演效果及工作量之間的關系，擬設7個測深供電點：1=40 m，2=60 m，3=80 m，4=100 m，5=120 m，6=140 m，7=160 m。為達到有效勘探深度，勘探區(qū)域為?50~250 m，取0~200 m進行剖面成圖，反演結果剖面圖中方框內(nèi)為模型投影圖。

3.1 正方柱體模型Ⅰ

模型Ⅰ為長×寬為20 m×20 m低阻、高極化正方柱體(見圖3)，中心位置為(100 m，?20 m)，電阻率為1=50 Ω?m，極化率為1=5.0%；背景電阻率為0=100 Ω?m，極化率為0=1.0%。正演模擬各點源對應的視極化率曲線如圖4所示。視極化率相對強度如圖5所示。視極化率最小二乘反演結果如圖6所示。

圖3 正方柱體模型示意圖

1—A1；2—A2；3—A3；4—A4；5—A5；6—A6；7—A7

圖5 正方柱體模型視極化率相對強度

圖6 正方柱體模型視極化率最小二乘反演結果

對比圖5和6可知：最小二乘反演、圓弧交匯?相對強度法能通過視極化率反演成圖能較好的體現(xiàn)出正方柱體的位置。視極化率相對強度法計算結果能較準確地反映出正方柱體的中心埋深位置；而最小二乘反演結果未能體現(xiàn)模型極化體中心，反演極化體中心位于實際模型中上部，但異常體處于模型投影范圍內(nèi)。

3.2 正方柱體模型Ⅱ

模型Ⅱ為長×寬為20 m×20 m正方低阻、高極化柱狀體(見圖7)，中心位置為(40 m，?17 m)，電阻率為1=50 Ω?m，極化率為1=5.0%；背景電阻率為0=100 Ω?m，極化率為0=1.0%。正演模擬各點源對應的視極化率曲線如圖8所示。視極化率相對強度如圖9所示。視極化率最小二乘反演結果如圖10所示。

圖7 正方柱體模型示意圖

1—A1；2—A2；3—A3；4—A4；5—A5；6—A6；7—A7

圖9 正方柱體模型正演視極化率相對強度圖

圖10 正方柱體模型視極化率最小二乘反演結果

由圖9可知：視極化率的相對強度中心與模型實際位置中心較吻合，但相對強度等值線不能確定模型的基本形態(tài)。相對強度中心周邊等值線向右成圓弧形傾斜，其原因是采用圓弧交匯法，與點源距離相等的各網(wǎng)格節(jié)點視電阻率、視極化率的取值相同，當點源分布相對模型中心不對稱時，各網(wǎng)格節(jié)點視電阻率、視極化率難以均衡，從而導致相對強度中心周邊等值線呈圓弧形發(fā)生明顯傾斜。

圖10表明：最小二乘反演能基本反映了模型的形態(tài)和位置；反演結果異常體頂面與實際模型頂面較吻合。但異常體中心位于實際模型中的中上部，模型的底界面無法通過異常體形態(tài)確定。圖10還表明：最小二乘反演結果受點源位置分布影響小。

3.3 水平板狀體模型

模型為長×寬為40 m×40 m水平低阻、高極化板狀體(見圖11)，中心位置為(100 m，?10 m)，電阻率1=50 Ω?m，極化率1=5.0%；背景電阻率0=100 Ω?m，極化率0=1.0%。正演模擬各點源對應的視極化率曲線如圖12所示。視極化率相對強度如圖13所示。視極化率最小二乘反演結果如圖14所示。

圖11 水平板狀體模型示意圖

1—A1；2—A2；3—A3；4—A4；5—A5；6—A6；7—A7

圖13 水平板狀體模型正演視極化率相對強度

圖14 水平板狀體模型視極化率最小二乘反演結果

對比圖13和圖14可知：最小二乘反演能很好的反映出水平低阻、極化板狀體的埋深、基本形態(tài)，異常體位置與模型實際位置對應關系較好；圓弧交匯?相對強度法能基本確定水平低阻、高極化板狀體位置，但視極化率相對強度中心位于低阻極化體模型位置正下部邊緣處，從視極化率相對強度等值線圖無法確定模型的基本形態(tài)。

3.4 組合模型

模型為2個長×寬為20 m×20 m正方柱形低阻、高極化模型(見圖15)，中心位置分別為(80 m，?18 m)、(120 m，?18 m)，電阻率1=50 Ω?m，極化率1=5.0%；背景電阻率0=100 Ω?m，極化率0=1.0%。正演模擬各點源對應的視極化率曲線如圖16所示，視極化率相對強度如圖17所示，視極化率最小二乘反演結果如圖18所示。

圖15 組合模型示意圖

1—A1；2—A2；3—A3；4—A4；5—A5；6—A6；7—A7

圖17 組合模型正演視極化率相對強度

圖18 組合模型視極化率最小二乘反演結果

圖17表明：組合模型視極化率相對強度中心與模型的實際位置中心相吻合；但模型的中垂線下方位置存在1個明顯的假異常，其原因是圓弧交匯法和相對強度法是建立在單個模型之上。對于組合模型，由于同一點源計算的視極化率曲線會產(chǎn)生2個極值點，從而進行圓弧交匯時就會產(chǎn)生其他相對強度，即為假異常，假異常通常隨著模型間的距離增大而變深。

圖18表明：最小二乘視極化率反演較好地反映了組合模型的實際位置，反演結果與實際模型較吻合，無明顯假異常存在，但反演結果難以反映出模型的底界面。

后期對高阻、高極化模型進行模型計算對比，其對比結論與低阻、高極化模型計算對比結果一致：最小二乘法反演結果優(yōu)于圓弧交匯?相對強度法。

4 結論

1) 固定點源激電測深具有簡單、高效的特點，其解釋方法圓弧交匯?相對強度法與最小二乘反演都能較準確地反映出異常體中心位置，這2種解釋方法直觀且其解釋結果與實際結果較符合。

2) 圓弧交匯?相對強度法對于單個模型能準確地反映出模型的中心埋深位置；但該解釋方法受場源位置影響大；對于稍復雜模型未能精確的體現(xiàn)出異常體的位置及形狀，對于組合模型通常會產(chǎn)生較明顯的假異常。

3) 最小二乘反演不受場源位置影響，能較好地反映出復雜模型的實際位置及輪廓，在對組合模型進行反演解釋時候，不會產(chǎn)生假異常；但該解釋方法難以體現(xiàn)模型的下底面位置。

[1] 李小康. 時間域激電二維正演研究[D]. 北京: 中國地質大學地球物理與信息技術學院, 2008: 27?47.
LI Xiaokang. 2D forward numerical simulation of ip method in time domain[D]. Beijing: China University of Geosciences. School of Geophysics and Geoinformation Systems, 2008: 27?47.

[2] 戴前偉, 肖波, 馮德山, 等. 基于二維高密度電阻率勘探數(shù)據(jù)的三維反演及應用[J]. 中南大學學報(自然科學版), 2012, 43(1): 293?310.
DAI Qianwei, XIAO Bo, FENG Deshan, et al. 3-D inversion of the high density resistivity method based on 2-D exploration data and its application[J]. Journal of Central South University (Science and Technology), 2012, 43(1): 293?310.

[3] 劉國興. 電法勘探原理與方法[M]. 北京: 地質出版社, 2005: 111?128.
LIU Guoxing, Electrical prospecting principle and method[M]. Beijing: Geological Publishing House, 2005: 111?128.

[4] 李金銘. 地電場與電法勘探[M]. 北京: 地質出版社, 2005: 258?260.
LI Jinming. Electric field and electrical prospecting[M]. Beijing: Geological Publishing House, 2005: 258?260.

[5] 李金銘, 魏文博, 陳本池, 等. 固定點源測深法定量解釋研究[J]. 物探與化探, 1997, 21(3): 187?196.
LI Jinming, WEI Wenbo, CHEN Benchi, et al. The quantitative interpretation of the fixed point source sounding methods[J]. Geophysical & Geochemical Exploration, 1997, 21(3): 187?196.

[6] 陳本池, 李金銘, 魏文博. 起伏地形條件下固定點源測深法定量解釋研究[J]. 現(xiàn)代地質, 1998, 12(1): 123?129.
CHEN Benchi, LI Jinming, WEI Wenbo. A study of quantitative interpretation for fixed point source sounding method under undulatory terrai[J]. Geoscience, 1998, 12(1): 123?129.

[7] 強建科, 阮百堯. 用FPS定量解釋法實現(xiàn)三維近似成像的探討[J]. CT理論與應用研究, 2003, 12(1): 13?16.
QIANG Jianke, RUAN Bairao. The study of 3-D IP anomaly by fixed point source sounding[J]. Study on the theory and Application of CT, 2003, 12(1): 13?16.

[8] 強建科, 羅延鐘, 熊彬. 固定點源測深激電異常研究[J]. 地球物理學進展, 2005, 20(4): 1176?1183.
QIANG Jianke, LUO Yanzhong, XIONG Bin. Study on the abnormal fixed point source sounding IP[J]. Progress in Geophysics, 2005, 20(4): 1176?1183.

[9] 熊彬, 阮百堯, 黃俊革. 直流電阻率測深中二維反演程序對三維數(shù)據(jù)的近似解釋[J]. 地球科學(中國地質大學學報), 2003, 28(1): 102?106.
XIONG Bin, RUAN Baiyao, HUANG Junge. Approximate interpretation of 3-D data using 2-D inversion program in the DC resistivity sounding[J]. Earth Science (Journal of China University of Geosciences), 2003, 28(1): 102?106.

[10] 黃俊革, 阮百堯. 直流電阻率測深中二維與三維反演結果的對比與分析[J]. 物探與化探, 2004, 28(5): 447?450.
HUANG Junge, RUAN Baiyao. An analytical comparison between 2D and 3D inversions in dc resistivity sounding[J]. Geophysical & Geochemical Explopation, 2004, 28(5): 447?450.

[11] 黃俊革, 王家林, 阮百堯. 三維高密度電阻率E-SCAN法有限元模擬異常特征研究[J]. 地球物理學報, 2006, 49(4): 1206?1214.
HUANG Junge, WANG Jialin, RUAN Baiyao. A study on FEM modeling of anomalies of 3-D high-density E-SCAN resistivity survey[J]. Chinese Journal of Geophysics, 2006, 49(4): 1206?1214.

[12] 劉海飛. 高密度數(shù)據(jù)處理方法研究[D]. 長沙: 中南大學信息物理工程學院, 2003: 30?60.
LIU Haifei. The processing method of High-density data[D]. Changsha: Central South University. School of Info-physics and Geomatics Engineering, 2003: 30?60.

[13] 黃俊革. 三維電阻率/極化率有限元正演模擬與反演成像[D]. 長沙: 中南大學信息物理工程學院, 2003: 20?60.
HUANG Junge. 3-Dresistivity/IP modeling and inversion based on FEM[D]. Changsha: Central South University. School of Info-physics and Geomatics Engineering, 2003: 20?60.

[14] 張大海, 王興泰. 二維視電阻率斷面的快速最小二乘反演[J]. 物探化探計算技術, 1999, 21(1): 1?12.
ZHANG Dahai, WANG Xingtai. Rapd least-squares inversion of 2-D apparent resistivity pseudosection[J]. Computing Techniques for Geophysical And Geochemical Exploration, 1999, 21(1): 1?12.

[15] 吳小平, 徐果明. 利用共軛梯度法的電阻率三維反演研究[J]. 地球物理學報, 2000, 43(3): 420?427.
WU Xiaoping, XU Guoming. Study 3-D resistivity inversion using conjugate gradient method[J]. Chinese Journal of Geophysics, 2000, 43(3): 420?427.

[16] 周竹生. 廣義共軛梯度算法[J]. 物探與化探, 1996, 20(5): 351?358.
ZHOU Zhusheng. The generalized conjugate gradient algorithm[J]. Geophysical and Geochemical Exploration, 1996, 20(5): 351?358.

Two different interpretation methods about FPS anomaly

XIAO Bo1, LIU Haifei2, DAI Qianwei2

(1. Guangdong Electric Power Design Institute, China Energy Engineering Group, Guangzhou 510600, China;2. School of Geosciences and Info-Physics, Central South University, Changsha 410083, China)

The plurality of geoelectric model instance was calculated based on the forward modeling of the finite element method, using the least square method and the program with VC++. Then the least square inversion results and arc intersection?relative intensity calculation results were compared. The results show that arc intersection?relative intensity method has high applicability in simple model, the inverse mapping effect is even superior to the least square optimization inversion results in determining the center of abnormal body position. The least square inversion method can adapt various geoelectric conditions on a variety of complex model inversion, and the inversion effects of the complex model and combined model are superior to the calculation results of the arc intersection?relative intensity method.

resistivity; polarizability; numerical simulation; finite element method; inversion; quantitative interpretation

P631

A

1672?7207(2015)02?0595?08

2014?02?13；

2014?04?20

國家自然科學基金資助項目(41374118，41174102，41074085)；教育部博士點基金資助項目(20120162110015)(Project (41374118, 41174102, 41074085) supported by the National Natural Science Foundation of China; Project (20120162110015) supported by the Doctoral Fund of Ministry of Education of China)

戴前偉，博士，教授，從事電磁法方法及理論、工程地球物理勘探等研究；E-mail：qwdai@mail.csu.edu.cn

10.11817/j.issn.1672-7207.2015.02.030

(編輯 趙俊)