在模拟点源二维电场时，为了减少计算量，避免进行三维模拟，常用的方法是用傅氏变换,在波数域中进行计算，即将解三维问题变为解多个二维问题，然后再进行傅氏反变换，获得空间中的电位.波数的选择很大程度上决定了计算的精度，同时波数的个数与计算时间是成正比的，一般来说波数越长，计算的精度越高，计算时间越长.徐世浙采用最优化方法(徐世浙，1988)，选取了5个波数，当模拟深度在100 m以内时达到了很好的精度，但是在实际勘探中，如果使用大极距大功率的电阻率法或者激电法，该组波数模拟精度不够,结果误差较大.韩思旭使用系数校正的方法(韩思旭等，2010)，采用徐世浙提出的5点波数，在数值模拟之后，对结果进行校正，在大极距条件下仍取得较好效果，但是对于不同的背景电阻率以及不同的构造，可能需要不同的校正系数，难以统一.王飞燕详细地讨论了波数的选取(王飞燕，2009)，并分别对垂直接触面模型、均匀半空间赋存单个异常体和均匀半空间赋存多个异常体三种情况提出了三组波数，取得了较好效果，同样该组波数没有考虑在大极距条件下的应用，而且三组波数难以统一，不利于反演的实现.

本文采用最优化方法，在极距为0.1 m到1000 m的范围，计算了两组波数，分别包含7个和9个，傅氏反变换的误差均在0.1%以内，通过模拟分析，该组波数精度较好.

假设地下构造为二维构造，u(x,y,z)为三维空间中的电位，U(x,y,k)是沿z(走向)方向进行傅氏变换之后在波数域中的电位.在空间中，取点源处为坐标原点，由于z是走向，根据对称性有u(x,y,z)=u(x,y,－z)，所以采用余弦形式，傅氏变换对为

ag=V ，

给定一组 k 值，在使=( Ι－V )^T( I－V )取得极小的情况可以确定一组 g ，其中 Ι = (1 … 1) ^T.

最优化方法求取波数 k－g 的过程为

(1)选取一组 r，给定一组初值k (同时也确定了波数的个数)；

(2)计算 V ，根据的极小，由 k求取对应的g；

(3)使用这一组 k－g 计算相对误差，如果精度没有达到要求则转向④，如果误差在阀值以内或者迭代次数达到设定的最大值则转向⑥；

(4)求 V 对于 k的偏导数矩阵；

(5)根据最小二乘原理，计算得到k 的修改量Δk，将Δk加入 k 中，然后执行②；

(6)保存计算得到的 k－g 退出循环.

设置r的范围为0.1～1000 m，设置阀值误差为0.1%，最大迭代次数为100次，计算了两组波数，分别包含7个波数和9个波数，误差均在0.1%以内.表 1、表 2给出了这两组波数：

表 1(Table 1)

表 1 7个波数 Table1 7 wave numbers, k and g K 0.000893 0.007744 0.028639 0.095209 0.337647 1.338806 5.766963 G 0.001869 0.007411 0.021980 0.073982 0.283508 1.210810 5.670225

表 1 7个波数 Table1 7 wave numbers, k and g

表 2(Table 2)

表 2 9个波数 Table2 9 wave numbers, k and g K 0.000704 0.005832 0.050142 0.128669 0.339511 0.018993 2.725006 0.939284 8.418160 G 0.001465 0.005285 0.030052 0.077935 0.214409 0.012413 1.879314 0.624406 6.537501

表 2 9个波数 Table2 9 wave numbers, k and g

地下电性结构为二维情况，沿走向变换到波数域中，总场法对应的变分问题为(徐世浙，1994)

使用矩形内剖分两个三角网格的方式，剖分区域，这样既方便形成地形，同时使得程序的编制变得简单，如图 1所示：

图 1

Fig. 1

图 1 区域有限元网格剖分(a) 水平地形区域剖分; (b) 起伏地形区域剖分. Fig.1 Division of region with FEM grids(a) Division of grid when flat-ground; (b) Division of grid when topography.

然后进行单元分析，得到刚度矩阵，也即方程组为

Ku=P ，

解该方程组，便可求得各个节点在波数域k中的电位值，最后进行傅氏反变换便可得到空间中的电位值，傅氏反变换的计算公式为

其中m代表点的位置，k_i，g_i即为上节中计算得到的波数.得到电位u后便可以根据具体模拟的装置，计算视电阻率.

为验证本文中的波数的有效性，首先在均匀大地与水平层状条件下，对比了文中的波数与前人提出的波数以及解析解，然后使用三维有限元模拟二维结构，对比使用本文波数与前人波数的二维模拟结果.

以下模拟计算使用的网格为：横向稀疏网格左右两边均为12个，网格大小为点距乘以系数，网格系数分别为1,1,1,1,2,2,4,4,8,16,32,64，观测区域网格为测点数，网格大小即为点距，垂向网格固定为35个，网格系数为1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,4,4,4,4,4,4,8,8,8,8,16,16,16,32,32,64,64.

3.1 均匀半空间

均匀半空间，电阻率为ρ=100 Ωm，采用对称四极测深，最大极距1000 m.采用上节中的有限元模拟方法，使用本文计算的7个以及9个波数与王飞燕计算的波数以及徐世浙计算的5个波数，以及解析解进行对比，如图 2所示：

图 2

Fig. 2

图 2 均匀半空间结果对比 Fig. 2 Homogenous half-space results contrast

图中analytical代表解析解结果，5_kg代表使用徐世浙5点波数的结果，wang_kg为使用王飞燕计算的波数的模拟结果，7_kg与9_kg分别为使用本文计算的7个与9个波数的模拟结果(以下图示相同).

如图可知，当极距大于200 m时，使用前人学者的波数其数值解与解析解误差较大，最大相对误差可达到30%，而使用本文计算的7个波数以及9个波数，极距到了1000 m时，仍然与解析解吻合较好.在均匀半空间中，使用该组波数，解析解与数值解均方误差在3%以内，同时7个波数与9个波数的模拟结果一致，它们之间的均方误差在1%以内.

3.2 K型地电模型

K型地电模型，如图 3所示，其中，第一层电阻率ρ₁=100 Ωm，层厚度H₁=30 m，第二层电阻率ρ₂=1000 Ωm，层厚度H₂=10 m，基底电阻率ρ₃=50 Ωm，采用对称四极测深，模拟结果如图 4所示：

图 3

Fig. 3

图 3 K型地电模型示意 Fig. 3 schematic of K-type geoelectric model

图 4

Fig. 4

图 4 K型地电模型模拟结果对比 Fig. 4 K-type geoelectric model results contrast

如图所示，使用本文计算的两组波数均与解析解吻合较好，相对误差在3%以内，从图中可以看出使用前人学者计算的波数，虽然曲线整体形态一致，但是误差较大.

3.3 H型地电模型

H型地电模型，如图 5所示，其中，第一层电阻率ρ₁=1000 Ωm，层厚度H₁=30 m，第二层电阻率ρ₂=50 Ωm，层厚度H₂=10 m，基底电阻率ρ₃=2000 Ωm，采用对称四极测深，模拟结果如图 6所示：

图 5

Fig. 5

图 5 H型地电模型示意 Fig. 5 schematic of H-type geoelectric model

图 6

Fig. 6

图 6 H型地电模型模拟结果对比 Fig. 6 H-type geoelectric model results contrast

如图所示，使用本文计算的波数与解析解吻合较好，同时发现使用7个波数与9个波数得到的结果基本一致，与解析解对比均能够控制误差在3%以内，所以建议采用7个波数进行模拟.

3.4 含单个异常体二维模型

背景电阻率ρ=100 Ωm，包含一个低阻异常体，模型具体参数如图 7所示：

图 7

Fig. 7

图 7 含低阻异常体模型示意 Fig. 7 schematic of model with a low resistivity abnormal body

采用对称四极装置，极距为20 m，50个电极，最大电极距为800 m，取模拟结果中，极距为600 m的测点，视电阻率曲线图如图 8所示，视电阻率剖面图如图 9所示：

图 8

Fig. 8

图 8 极距600 m的视电阻率曲线图 Fig. 8 Apparent resistivity curves when polar distance is 600 m

图 9

Fig. 9

图 9 含异常体模型模拟结果剖面图 Fig. 9 The analog results profile map of model with a low resistivity abnormal body

图 8，图 9中3D为使用有限单元法六面体剖分三维模拟的结果(模拟二维模型)，与使用不同波数的二维模拟结果进行对比.通过观察图 8，可以发现使用文中计算的波数与三维模拟结果吻合较好，而前人学者提出的波数均出现一定的偏离.剖面图(图 9)中，本文的波数与三维模拟结果在形态和数值上均吻合(使用9个波数得到的结果与7个波数是一致的，在此仅给出7个波数模拟的结果，下同).

3.5 含两个异常体二维模型

背景电阻率ρ=100 Ωm，包含两个异常体，模型具体参数如图 10所示：

图 10

Fig. 10

图 10 含两个异常体模型示意 Fig. 10 schematic of model with two abnormal bodies

采用对称四极装置，极距为20 m，一共50个电极，最大电极距为800 m，取模拟结果中，极距为600 m的测点，视电阻率曲线图如图 11所示，视电阻率剖面图如图 12所示：

图 11

Fig. 11

图 11 极距600 m的视电阻率曲线图 Fig. 11 Apparent resistivity curves when polar distance is 600 m

图 12

Fig. 12

图 12 含异常体模型模拟结果剖面图 Fig. 12 The analog results profile maps of model with two abnormal bodies

从图 12可以看出，在整体形态上，使用文中的波数，与三维模拟结果一致，而使用前人学者计算的波数，视电阻率在深部误差较大.

本文通过最优化方法，选取点源二维电场傅氏反变换中的波数以适应极距达到1000 m.通过计算得到了两组波数，分别包含7个波数和9个波数.采用文中波数，模拟层状大地，两组波数的精度均较高，与解析解的相对误差在3%以内，而前人学者提出的波数在极距大于200 m时与解析解对比，误差较大.同时得出7个波数得到的结果与9个波数得到的结果基本一致，所以实际模拟中建议采用7个波数.最后本文通过模拟二维模型与三维模拟结果对比，再次验证文中波数的正确性.