2.5维直流电法正演中Fourier逆变换离散波数的最优化选取

潘克家^{1, 2}，汤井田¹

(1. 中南大学 有色金属成矿预测教育部重点实验室，地球科学与信息物理学院，湖南 长沙，410083

2. 中南大学 数学与统计学院，湖南 长沙，410083)

摘要：从优化模型和计算方法2方面改进点源2.5维直流电法正演中Fourier逆变换离散波数的最优化选取方法。首先利用均匀半空间点源电位的精确解，基于最小二乘法给出计算离散波数的改进非线性最优化问题；然后，利用差分进化(DE)算法进行求解；最后，研究参与计算的电极距数对模拟精度的影响。对具有解析解的典型地电模型，通过与已有文献计算结果进行比较，验证方法的可行性。研究结果表明：增加参与计算的电极距数可有效提高视电阻率曲线近源处的计算精度，并能保证较大电性差异情形下的计算精度；与现有离散波数相比，本文方法得到的波数具有更高的精度和更大的适用范围。

关键词：直流电法；有限单元法；Fourier逆变换；离散波数；差分进化算法

中图分类号：P631           文献标志码：A         文章编号：1672-7207(2013)07-2819-08

Optimized selection of discrete wavenumbers for inverse Fourier transform in 2.5D DC resistivity modeling

PAN Kejia^{1, 2}, TANG Jingtian¹

(1. Key Laboratory of Metallogenic Prediction of Nonferrous Metals, Ministry of Education,

School of Geosciences and Info-Physics, Central South University, Changsha 410083, China;

2. School of Mathematics and Statistics, Central South University, Changsha 410083, China)

Abstract: The optimized selection method was improved in terms of both the optimization model and computational method for selecting the discrete wavenumbers of inverse Fourier transform, which arise in the 2.5-dimensional direct current resistivity modeling with a point current source. Firstly, with an analytical solution of electric potential produced by a point source in semi-infinite homogeneous medium, an improved nonlinear optimization problem for computing these wavenumbers based on the method of least squares was presented, then it was resolved using the differential evolution (DE) algorithm. Finally, the computational accuracy against the impact of the number of electrode spacings included in the simulation was discussed. For a number of typical geoelectric models with analytical solutions, the method was justified through comparing the computed results to those in literatures. The results show that increasing the number of electrode spacings not only dramatically reduces the error of the apparent resistivity near the point source, but also ensures the computational accuracy for models with a large conductivity contrast. The discrete wavenumbers obtained from the proposed method have a higher accuracy and wider range of applications compared with existing ones.

Key words: direct current modeling; finite element method; inverse Fourier transform; discrete wavenumbers; differential evolution algorithm

利用有限元法(FEM)或者有限差分法(FDM)求解2.5维点源电阻率法正演问题时，需要求Fourier逆变换，而通过数值方法计算Fourier逆变换时，波数的选取是保证计算精度和节省计算时间的关键问题^[1-8]。罗延钟等^[9]根据零阶修正贝塞尔函数的特点，构造线性和负指数函数，对其进行分段积分来完成Fourier逆变换。该方法选择的波数呈等比数列，称之为等比离散波数。由于该方法仅考虑最大和最小电极距，所以，用等比离散波数进行Fourier逆变换时，需要较多的波数才能满足模拟精度，显然正演的计算量随着波数个数的增加而线性增加。为此，徐世浙等^[10-11]采用最优化方法计算离散波数，优点在于将电极距序列的每个电极都参与计算，通过迭代法求解优化问题使得Fourier逆变换的计算结果达到最优，由此计算出来的离散波数称之为最优化离散波数。该方法是一种非常有效的Fourier逆变换方法，已在2.5维电法正演中得到广泛运用^[12-20]。柳建新等^[21]在最优化波数的基础上，研究了离散波数初值的选定以及偏导数矩阵的计算问题，使得该方法在计算量和计算精度方面有了进一步的改善。近年来，Li等^[22-23]将计算离散最优化波数中的均匀地层参考模型改为水平分层介质模型，提出增强最优化波数的计算方法，进一步提高了离散波数的精度。综合现有文献，基于最优化波数的2.5维电法正演整体精度较高，但在点电源附近精度仍然比较差^{[17, 23]}。尽管采用二次场(背景场+异常场)的方法可以避开电源点的奇性^[24]，提高点源附近的计算精度，但是，此方法难以推广到考虑起伏地表的直流电法正演问题中。因此，进一步提高2.5维直流电法正演的模拟精度具有重要的意义。本文作者从模型和算法两方面改进最优化波数的计算方法，解决点源2.5维直流电法模拟精度(尤其源点附近)不高的问题。对具有解析解的水平地电模型和三层垂直分层Dike模型，与文献[6]和文献[17]中离散波数的模拟结果进行比较与分析，验证本文改进最优化方法得到的波数具有更高的精度和更大的适用范围，且适用于大电性差异地电模型。

1  点源2.5维电场边值问题及其Fourier逆变换

考虑具有一定走向的二维构造地电模型，记u(x, y, z)为点源二维介质电场问题中的三维电位，U(x, y, k)为其在z (走向)方向进行Fourier变换后的二维电位，即

       (1)

其中：k为波数。实际上，式(1)为余弦变换，这是因为二维地电构造关于z=0对称，电位u(x, y, z)关于z为偶函数。二维电位U(x, y, k)满足如下带参数(波数k)的二维变系数非齐次亥姆霍茨(Helmholtz)方程边值问题^{[6, 17]}:

       (2)

其中：为地表以下求解区域；为地表-空气界面；为截断无穷远边界；r为点电源到无穷远边界的距离；为电导率；I为电流；n为边界的外法向；为狄拉克函数，K₀和K₁分别为第二类零阶、一阶修正贝塞尔函数。

由Fourier逆变换知：

      (3)

理论上，先求解偏微分方程边值问题(2)，再由式(3)进行Fourier逆变换即可得点源二维介质电场问题中的三维电位u。然而，边值问题(2)只能通过有限元(FEM)等方法数值求解，即只能得到一系列离散波数{k_i, i=1, 2, …, N}对应的U(x, y, k_i)的近似解U_h(x, y, k_i)。因此，如何选取最优的离散波数{k_i}是非常值得研究的一个问题。

2  离散波数选取非线性优化问题

文献[10-11]给出了用最优化方法确定离散波数的基本原理。大多数情形下，只需考虑通过点电源剖面(主剖面)上的电位，此时z=0，由式(3)可知：

             (4)

将式(4)写成如下形式：

               (5)

其中：，为主剖面上的点至电源点的距离；N为离散波数的个数；k_i为离散波数；g_i为相应的权系数。离散波数k_i及相应的权系数g_i需要适当选取使得式(5)对一定范围内的r尽可能精确成立。然而，一般地电模型得不到u和U的精确表达式，故无法直接通过式(4)确定k_i和g_i。

在均匀半空间情形下， ，将其代入式(1)得：

          (6)

将式(6)代入式(4)得：

      (7)

可利用式(7)确定k_i及g_i，使得式(7)对一定范围内的r尽可能精确地成立。

选取电极距r_j (j=1, 2, …, M)，根据非线性最小二乘法，可通过求解如下带非负约束的非线性优化问题确定离散波数k_i及相应系数g_i：

             (8)

其中：M和N分别为电极距数和波数个数。

非线性优化模型(8)与文献[11]的不同，文献[11]中采用的办法是将式(7)改写为

                  (9)

转化为如下非线性优化问题：

      (10)

值得注意的是，式(7)和式(9)在数学上是等价的，但由这2式分别导出的优化问题(8)和(10)不等价，计算结果也完全不同。考虑电位的误差，电位为O(1/r)的量级，因此优化模型(8)更加合理。而优化问题式(10)过多地考虑了较大r处的精度，忽略了近源处的计算精度。

求解优化问题式(8)的数值方法与文献[11]中的迭代求解法不同。本文采用的方法是将离散波数k_i及相应系数g_i当成同等地位待定参数，直接利用差分进化算法(differential evolution，DE)算法^[25-27]求解带非负约束的非线性优化问题(8)，同时确定k_i和g_i。DE算法本质是一种基于实数编码的具有保优思想的贪婪遗传算法，具有很好的全局寻优能力，不需要设定波数初值，也避免了繁琐的计算；并且非负约束条件自动满足，同时确定k_i及g_i具有更高的精度。

式(6)采用均匀半空间模型作为参考模型，事实上，根据实际问题地下电性结构，可采用同样具有解析解的水平层状模型或纵向垂直分层模型作为参考模型，将解析解代入式(4)和(6)得到类似的非线性优化问题，进一步提高离散波数的精度^[22-23]。

3  数值计算结果

本文的数值测试平台为CPU i5 760 2.8G，RAM 2G；测试系统为Compaq Visual Fortran Professional Edition 6.6，自带IMSL国际数学和统计链接库，调用贝塞尔函数比较方便。

为了说明本文方法所得的离散波数k_i及相应权系数g_i的可靠性，首先研究了电极距个数对离散波数精度的影响，然后，对具有解析解的水平地电模型以及3层垂直分层Dike模型，并与文献[6]中最优化波数(见表1)，以及文献[17]中最优化波数(见表2)的计算结果进行比较与分析。

3.1  电极距数的影响

首先以对称四极测深常用极距r_j序列为例：r_j=1.5，2.5，4，6，9，15，25，40，65，100，150和220 m，共12个极距。为方便与已有文献比较，取离散波数个数N=5，利用DE算法求得非线性优化问题式(8)，求得离散波数及权系数如表3所示。

取电极距r_j=j (1≤j≤220)，共220个极距，保持与对称四极测深常用极距序列的范围不变。同样取离散波数个数N=5，利用DE算法求解非线性优化问题式(8)，求得最优化离散波数及相应系数如表4所示。

图1所示为4组不同波数的相对误差。相对误差e为

          (11)

表1  文献[6]中离散波数及其系数

Table 1  Wavenumbers and weights in Ref.[6]

表2 文献[17]中离散波数及其系数

Table 2  Wavenumbers and weights in Ref.[17]

表3  12极距得到的波数及其系数

Table 3  Wavenumbers and weights obtained from 12 polar distances

表4  220极距得到的波数及其系数

Table 3  Wavenumbers and weights obtained from 220 polar distances

              (12)

从图1可看出，在2＜r＜60范围内，4组波数对应的计算值的相对误差均比较小，控制在0.5%以内。但在近源处，文献[6]中的波数的相对误差比较大，尤其在r=1时，相对误差达到6.04%；当r＞100时，文献[6]中的波数的相对误差随着r的增大迅速增大，至r=200时可达10.10%。文献[17]中的波数在近源处精度比较高，但整体范围内误差波动比较大，平均相对误差达到0.19%，文献[6]中的波数、12极距波数和220极距波数平均相对误差分别为0.17%，0.04%和0.03%；并且当r＞60时，误差亦随着r的增大激增，当r=200时，达到14.84%。本文方法得到的12极距波数的误差相对于文献[6]中的波数、文献[17]中的波数有明显减小，平均相对误差减小为原来的1/5，但近源处相对误差仍比较大，达到2.24%；而220极距波数的误差，在1＜r＜200范围内均非常小，最大相对误差仅为0.5%。从图1可以看出，本文方法显著提高了近源处的精度，同时对远源处精度的提高更加明显，这主要是因为本文方法采用了比较多的电极距，且设置其范围比较宽(1≤r≤220)。总之，本文改进优化方法得到的离散波数精度非常高，适用范围较已有波数更大，并且增加极距个数及其范围有望获得更加精确、稳定的离散波数。

图1  不同波数的相对误差

Fig.1  Relative errors of different wavenumbers

3.2  水平层状地电模型

为了进一步验证本文方法的正确性，将最优化波数应用于求解如图2和图3所示的水平2层和3层地电模型。其中2层模型取自文献[6]，3层模型取自文献[28]。采用单极供电(电流强度为I=1 A)，供电电极A位于原点。

图2  二层地电断面示意图

Fig.2  Schematic diagram of two-layered geoelectric section

图3  三层地电断面示意图

Fig.3  Schematic diagram of three-layered geoelectric section

图4和图5所示分别为两层地电模型和三层地电模型采用不同离散波数得到的数值模拟结果。从图5可以看出：220波数对应的视电阻率与解析解基本吻合，相对误差控制在0.73%以内；而采用文献[6]和[17]中的波数时，点源附近和r＞100 m处误差比较大，如r=200 m时，相对误差分别达到11.27%和16.36%，r= 1 m时相对误差分别达到3.15%和0.92%。三层水平层状地电模型与二层模型相比，可以得到完全类似的结论。值得指出的是：文献[17]波数的相对误差波动较大。

进一步分析可知直流电法2.5维有限元正演的误差主要来自2个方面：一是模型本身(如截断边界条件等)及其有限元离散带来的误差，二是离散Fourier逆变换，即用式(5)近似代替式(3)带来的误差。从图4(b) 和图5(b)可以看出：误差曲线的形状与图1中曲线形状基本相同，变化趋势也一样。这说明2个算例中有限元逼近较准确，引起的误差相对较小，计算误差主要来自离散Fourier逆变换。因此，图4(b)和图5(b)中相对误差曲线反映了3组最优化波数的精确程度及其适用范围。

图4  二层地电模型数值解与解析解的对比

Fig.4  Comparison of analytical and numerical solutions for two-layered geoelectric model

图5  三层地电模型数值解与解析解的对比

Fig.5  Comparison of analytical and numerical solutions for three-layered geoelectric model

3.3  大电性差异垂直分层模型

Dike模型^[29]如图6所示。该模型为垂直的三层介质，中间层厚度很薄。考虑大电性差异Dike模型：取中间层电阻率为ρ₂=1Ω·m，厚度为5 m，左边电阻率ρ₁=100Ω·m，右边电阻率ρ₃=50Ω·m。单点电源A位于坐标原点，采用二极AM测深装置。3组最优化波数下的视电阻率及其相对误差如图7所示。从图7可以看出：对大电性差异的Dike模型，利用本文最优化波数的模拟结果与精确解仍然基本重合，相对误差控制在5.58%以内。而另外2组离散波数的模拟结果仅在2~20 m的范围内相对误差比较小，经过低阻薄层(x＞25 m)后，相对误差迅速增大，到100 m时相对误差分别达到13%和17%左右，根本不能满足实际的精度需求。

图6  Dike模型示意图

Fig.6  Schematic diagram of the Dike model

图7  垂直Dike模型数值解与解析解的对比

Fig.7  Comparison of analytical and numerical solutions for vertical Dike model

4  结论

(1) DE算法求解离散波数的最优化问题，无需设定初值、也避免繁琐的计算，且同时确定离散波数及其权系数，精度比较高。

(2) 增加用于计算离散波数的极距数能有效提高离散波数的精度，增大极距的范围能扩大所得离散波数的适用范围。

(3) 与已有离散波数相比，本文改进方法得到的波数具有更高的精度和更大适用范围，尤其在近源处和远源处精度明显提高。

(4) 本文改进最优化方法确定的离散波数适用于大电性差异模型，而已有文献中波数对大电性差异模型误差比较大，不能满足实际的精度需求。

参考文献：

[1] Coggon J H. Electromagnetic and electrical modeling by the finite element method[J]. Geophysics, 1971, 36(1): 132-145.

[2] Dey A, Morrison H F. Resistivity modeling for arbitrarily shaped two dimensional structures[J]. Geophysical Prospecting, 1979, 27(1): 106-136.

[3] Mundry E. Geoelectrical model calculations for two-dimensional resistivity distributions[J]. Geophysical Prospecting, 1984, 32(1): 124-131.

[4] Mufti I R. Finite-difference resistivity modeling for arbitrarily two-dimensional structures[J]. Geophysics, 1976, 41(1): 62-78.

[5] 罗延钟, 张桂青. 电子计算机在电法勘探中的应用[M]. 武汉: 武汉地质学术出版社, 1987: 10-15.

LUO Yanzhong, ZHANG Guiqing. Application of electric computers in electrical exploration[M]. Wuhan: Wuhan College of Geology Press, 1987:10-15.

[6] 徐世浙. 地球物理中的有限元法[M]. 北京: 科学出版社, 1994: 178-188.

XU Shizhe. The finite element method in geophysics[M]. Beijing: Science Press, 1994: 178-188.

[7] Xu S Z. The boundary element method in geophysics[M]. Geophysical Monograph Series. Tulsa: Society of Exploration, 2001: 12-20.

[8] Xu S Z, Zhao S, Ni Y. A boundary element method for 2-D dc resistivity modeling with a point current source[J]. Geophysics, 1998, 63(2): 399-404.

[9] 罗延钟, 孟永良. 关于用有限元法对二维构造作电阻率法模拟的几个问题[J]. 地球物理学报, 1986, 29(6): 613-621.

LUO Yanzhong, MENG Yongliang. Some problems on resistivity modeling for two-dimensional structures by the finite element method[J]. Chinese J Geophys, 1986, 29(6): 613-621.

[10] 徐世浙. 点电源二维电场问题中付氏变换的波数k的选择[J]. 物探化探计算技术, 1988, 10(3): 235-239.

XU Shizhe. Selection of wave number k in Fourier inverse transformation for point source and 2-D electric field problems[J]. Computing Techniques for Geophysical and Geochemical Exploration, 1988, 10(3): 235-239.

[11] Xu S Z, Duan B C, Zhang D H. Selection of the wavenumbers k using an optimization method for the inverse Fourier transform in 2.5D electrical modeling[J]. Geophysical Prospecting, 2000, 48(5): 789-796.

[12] 徐世浙. 电导率分段线性变化的水平层的点电源电场的数值解[J]. 地球物理学报, 1986, 29(1): 84-90.

XU Shizhe. A numerical method for solving electric field of point source on a layered model with linear change of conductivity in each layer[J]. Chinese J Geophys, 1986, 29(1): 84-90.

[13] 阮百尧, 徐世浙. 电导率分块线性变化二维地电断面电阻率测深有限元数值模拟[J]. 地球科学: 中国地质大学学报, 1998, 23(3): 303-307.

RUAN Baiyao, XU Shizhe. FEM for modeling resistivity sounding on 2D geoelectric model with line variation of conductivity within each block[J]. Earth Science: Journal of China University of Geosciences, 1998, 23(3): 303-307.

[14] 周熙襄, 钟本善, 严忠琼, 等. 电法勘探正演数值模拟的若干结果[J]. 地球物理学报, 1983, 26(5): 479-491.

ZHOU Xixiang, ZHONG Benshan, YAN Zhongqiong, et al. Some results of resistivity modelling[J]. Chinese J Geophys, 1983, 26(5): 479-491.

[15] 熊彬, 阮百尧. 电位双二次变化二维地电断面电阻率测深有限元数值模拟[J]. 地球物理学报, 2002, 45(2): 285-295.

XIONG Bin, RUAN Baiyao. A numerical simulation of 2-D geoelectric section with biquadratic change of potential for resistivity sounding by the finite element method[J]. Chinese J Geophys, 2002, 45(2): 285-295.

[16] 底青云, 王妙月. 稳定电流场有限元法模拟研究[J]. 地球物理学报, 1998, 41(2): 252-260.

DI Qingyun, WANG Miaoyue. The real-like 2D FEM modeling research on the field characteristics of direct electric current field[J]. Chinese J Geophys, 1998, 41(2): 252-260.

[17] 汤井田, 王飞燕, 任政勇. 基于非结构化网格的2.5-D直流电阻率自适应有限元数值模拟[J]. 地球物理学报, 2010, 53(3): 708-716.

TANG Jingtian, WANG Feiyan, REN Zhengyong. 2.5-D DC resistivity modeling by adaptive finite-element method with unstructured triangulation[J]. Chinese J Geophys, 2010, 53(3): 708-716.

[18] 郑洲顺, 黄辉, 张继锋, 等. 2.5-D稳定电流场h-自适应有限元的后验误差分析[J]. 地球物理学进展, 2009, 24(2): 663-667.

ZHEN Zhoushun, HUANG Hui, ZHANG Jifeng, et al. Posteriori error estimation analysis of the h-adaptive finite element in the 2.5-dimensional direct current field[J]. Progress in Geophysics, 2009, 24(2): 663-667.

[19] 王武, 赵林, 汤井田. 基于ANSYS的2.5维直流电法正演模拟[J]. 地球物理学进展, 2010, 25(2): 516-524.

WANG Wu, ZHAO Lin, TANG Jingtian. 2.5D forward modelling of the direct current method using ANSYS[J]. Progress in Geophysics, 2010, 25(2): 516-524.

[20] 张东良, 孙建国, 孙章庆. 2维和2.5维起伏地表直流电法有限差分数值模拟[J]. 地球物理学报, 2011, 54(1): 234-244.

ZHANG Dongliang, SUN Jianguo, SUN Zhangqing. Finite- difference DC electrical field modeling on 2D and 2.5D undulate topography[J]. Chinese J Geophys, 2011, 54(1): 234-244.

[21] 柳建新, 刘海飞. 计算最优化离散波数的优化算法[J]. 物探化探计算技术, 2005, 27(1): 34-38.

LIU Jianxin, LIU Haifei. The optimum algorithm for the calculation of the optimized discrete wave-number[J]. Computing Techniques for Geophysical and Geochemical Exploration, 2005, 27(1): 34-38.

[22] Li Y G, Spitzer K. Three-dimensional DC resistivity forward modelling using finite elements in comparison with finite- difference solutions[J]. Geophys J Int, 2002, 151: 924-934.

[23] Tang J T, Wang F Y, Xiao X, et al. 2.5-D DC resistivity modeling considering flexibility and accuracy[J]. Journal of Earth Science, 2011, 22(1): 124-130.

[24] 徐世浙, 刘斌, 阮百尧. 电阻率法中求解异常电位的有限元法[J]. 地球物理学报, 1994, 37(S2): 511-515.

XU Shizhe, LIU Bin, RUAN Baiyao. The finite element method for solving anomalous potential for resistivity surveys[J]. Chinese J Geophys, 1994, 37(S2): 511-515.

[25] Storn R, Price K. Differential evolution: A simple and efficient heuristic for global optimization over continuous spaces[J]. Journal of Global Optimization, 1997, 11: 341-359.

[26] 潘克家, 王文娟, 谭永基, 等. 基于混合差分进化算法的地球物理线性反演[J]. 地球物理学报, 2009, 52(12): 3083-3090.

PAN Kejia, WANG Wenjuan, TAN Yongji, et al. Geophysical linear inversion based on hybrid differential evolution algorithm[J]. Chinese J Geophys, 2009, 52(12): 3083-3090.

[27] 潘克家, 陈华, 谭永基. 基于差分进化算法的核磁共振T2谱多指数反演[J]. 物理学报, 2008, 57(9): 5956-5961.

PAN Kejia, CHEN Hua, TAN Yongji. Multi-exponential inversion of T2 spectrum in NMR based on differential evolution algorithm[J]. Acta Physica Sinica, 2008, 57(9): 5956-5961.

[28] 李长伟, 阮百尧, 吕玉增, 等. 点源场井-地电位测量三维有限元模拟[J]. 地球物理学报, 2010, 53(3): 729-736.

LI Changwei, RUAN Baiyao, L Yuzeng, et al. Three- dimensional FEM modeling of point source borehole-ground electrical potential measurements[J]. Chinese J Geophys, 2010, 53(3): 729-736.

[29] 潘克家, 汤井田, 胡宏伶, 等. 直流电阻率法2.5维正演的外推瀑布式多重网格法[J]. 地球物理学报, 2012, 55(8): 2769-2778.

PAN Kejia, TANG Jingtian, HU Hongling, et al. Extrapolation cascadic multigrid method for 2.5D direct current resistivity modeling[J]. Chinese J Geophys, 2012, 55(8): 2769-2778.

(编辑  赵俊)

收稿日期：2012-06-18；修回日期：2012-09-20

基金项目：国家自然科学基金资助项目(41204082)；国家科技专项(Sino Probe-03)；高等学校博士学科点科研基金资助项目(20120162120036)；中国博士后科学基金资助项目(2011M501295)

通信作者：汤井田(1965-)，男，江苏连云港人，博士，教授，从事电磁场理论及应用、数字信号处理等研究；电话：13507317396；E-mail: jttang@csu.edu.cn