CN105740204A - 不规则地形下低频段大地电导率快速反演方法 - Google Patents

Abstract

本发明公开了一种不规则地形下低频段大地电导率快速反演方法，针对待确定大地电导率取值的区域，在确定测量区域后，利用长波场强/时延测量系统获取部分观测点的时延测量值；以积分方程方法作为正演算法，计算得到各观测点的时延理论预测值；采用遗传算法为优化方法，基于GPU并行技术，以时延测量值与预测值之间误差最小为优化准则，快速反演得到大地电导率分布数据。本发明方法解决了现有大地电导率反演结果中包含地形影响等效成分、无法反应地面实际导电性能的问题，本发明方法能够有效提高反演精度与反演效率。

Description

不规则地形下低频段大地电导率快速反演方法

技术领域

本发明属于电波传播及地球物理反演技术领域，具体涉及一种不规则地形下低频段大地电导率快速反演方法。

背景技术

不同的地面电特性(低频段主要是指大地电导率)会对低频地波传播性能产生非常大的影响。获取传播路径上详细、准确的大地电导率分布信息，是提高低频地波传播时延/场强理论预测精度的前提和物质基础。地波场强/时延反演方法是获取大地电导率最基本和最普遍的一种方法，其依据是：地波在传播的过程中因地表的影响而产生衰减和相位延迟，其衰减程度和相位延迟大小与信号频率、传播距离、大地电导率等因素有关。因此，可以通过测量来自地波导航台、授时台或者广播台的发播信号的场强或者时延值来反推传播路径上的大地电导率分布信息。其过程是：利用观测数据(场强或时延)，以复杂的电波传播理论为依据，以观测值与理论计算值之间的误差最小为优化准则，求取相应的地球物理模型参数(大地电导率分布)。因此，确定观测数据和地球物理模型参数之间的函数关系，实现正演计算是大地电导率反演的前提。目前，大地电导率传统反演过程中多是将传播路径假定为光滑地面，基于均匀光滑地面模型或分段均匀光滑地面模型，采用相应的算法(如平地面公式、Fock地波绕射公式、Millington经验公式等)作为正演算法来进行大地电导率的反演。

然而，对于实际的传播路径而言，地面并非光滑的，会存在高山、丘陵、峡谷等起伏变化，即不规则地面。地形的起伏变化会使电磁波发生一系列反射、绕射、散射等复杂变化，从而影响低频地波传播性能。那么，对于实际不规则地面而言，采用传统的大地电导率反演方法所反演得到的大地电导率中则包含了地形影响的等效成份，实为一“等效值”。其难以反映地面的实际电特性，且同一区域一条路径的反演结果无法推广应用到其它路径。

发明内容

本发明的目的是提供一种不规则地形下低频段大地电导率快速反演方法，该方法能有效剥离现有大地电导率中的地形影响等效成分，还原地面的实际导电特性，提高中国大地电导率数据精度。

本发明所采用的技术方案是，一种不规则地形下低频段大地电导率快速反演方法，具体步骤如下：

步骤1、针对待反演区域，选定发射台即后续测量系统所接收信号来源，确定测量区域；

步骤2、利用长波传播时延测量系统，获取各观测点的传播时延测量值T_M；

步骤3、利用先验知识，结合GIS技术获取各观测点所对应传播路径参数信息：包括电导率分区数目、各分区内电导率取值范围、各分区内的传播距离以及传播路径上的高程信息；

步骤4、根据步骤3得到的传播路径参数信息，在CPU上初始化参数；

步骤5、产生大地电导率初始种群，形成各观测点所对应模型文件；

步骤6、将步骤4初始化参数传递到GPU，完成设备初始化、数据准备，启动核函数；

步骤7、计算大地电导率初始种群适应度；

步骤8、判断是否达到结束条件，若达到，执行步骤10，若未达到，执行步骤9后返回步骤7；

步骤9、遗传算子操作，产生新个体，形成新的大地电导率种群；

步骤10、将反演结果传回CPU。

本发明的特点还在于，

步骤1中待反演区域指地球表层某一大地电导率取值未知且待反演的地理区域；

发射台的选定：选取距离待反演区域较近、信号质量较好的长波导航台或授时台为信号来源，用于后续进行接收、测量；

测量区域的确定：以所选定发射台为中心，做间隔0.01度、辐射距离为1000km的射线，则所有穿过待反演区域的射线中最外侧两条射线所夹不规则区域即为测量区域。

步骤2中观测点的选取标准：在步骤1所限定的测量区域内，综合考虑实际地形及交通便利情况，合理规划路线，进行定点测量；测量点应远离山区或峡谷等地形起伏区域；测量点后向地形应较平坦；

长波传播时延测量系统主要包括：

1)长波接收机：提供从发射台到测量点的以GPS时钟为基准的到达时刻；

2)高精度GPS接收机：提供测量点的位置信息、时间信息、1pps信号；

3)计数器：比对长波接收机与GPS接收机1pps信号的时间间隔；

4)工控机：监测长波接收机、GPS接收机、计数器的运行状态，采集和处理数据。

步骤3中获取各观测点所对应传播路径参数信息过程如下，以第i个观测点为例：

1)提取观测点所在传播路径的高程信息h_i(即地形起伏变化情况)；

2)结合高程信息和现有大地电导率地图进行大地电导率分区，并给出分区数目N和各分区内电导率取值范围l表示第l个分区；

3)计算传播路径上各段路径长度d_ij，即第i个观测点第j段大圆距离。

步骤4中初始化的参数有：积分步长Δx、大地电导率种群规模S、迭代次数P、目标函数阈值ξ、交叉概率P_c、变异概率P_m；大地电导率分区数目n、观测点个数m；各观测点时延测量数据各观测点所对应传播路径的高程数据h＝[h₁h₂…h_m]。

步骤5中产生大地电导率初始种群及建立观测点模型文件过程如下：

1)种群中第k个个体σ_k的第l个参数σ_kl在其取值范围内随机生成，式中r为随机数；

${\mathrm{\σ}}_{kl}={\mathrm{\σ}}_l^{LB}+r\×({\mathrm{\σ}}_l^{UB}-{\mathrm{\σ}}_l^{LB})$

2)按照1)的方法，形成第k个个体σ_k＝[σ_k1σ_k2...σ_kn]，此即为一种大地电导率分布模型；

3)重复步骤1)和2)，得到S个个体模型，形成大地电导率初始种群Ab|₀

4)结合步骤3中得到的各观测点所在传播路径上各段路径长度d_ij，则第i个观测点的模型为：

步骤7中定义第k个个体σ_k的适应度为

$Fit\left({\mathrm{\η}}_k\right)=\frac{max\left(\mathrm{\η}\right)-{\mathrm{\η}}_k}{S\×\[max\left(\mathrm{\η}\right)-mean\left(\mathrm{\η}\right)\]}$

式中，η_k为σ_k的目标函数，max(η)和mean(η)分别为当前种群下所有个体所对应目标函数的最大值和均值；

则目标函数η_k的计算过程分为如下步骤：

步骤7.1、采用积分方程方法并行计算第t代种群Ab|_t下每一个观测点在每一种模型下的时延预测值T_C；

则第i个观测点对应第k个模型时的时延预测值为

$T_{C_k}^i=\frac{\arg \left(W_g\right)}{\mathrm{\ω}}\×{10}^6\mathrm{\μ}s;$

式中，

$\begin{array}{c}{W}_g\left(P\right)=1-R_0\sqrt{\frac{1}{j\mathrm{\λ}x}}\{\underset{0}{\overset{x_1}{\∫}}{W}_g\left(Q\right)\frac{F(x,l)}{\sqrt{l(x-l)}}dl+\underset{x_1}{\overset{x_{n-1}}{\∫}}{W}_g\left(Q\right)\frac{F(x,l)}{\sqrt{l(x-l)}}dl+\underset{x_{n-1}}{\overset{x_n}{\∫}}{W}_g\left(Q\right)\frac{F(x,l)}{\sqrt{l(x-l)}}dl\}\\ =1-R_0\sqrt{\frac{1}{j\mathrm{\λ}x}}(I_1+I_2+I_3)\end{array}$

$F(x,l)=\frac{l(x-l)}{R_1{R}_2}{e}^{{\mathrm{jk}}_0(R_1+R_2-R_0)}\[{\mathrm{\Δ}}_g+(1+\frac{j}{k_0{R}_2})\]\frac{\[z_Q-z_{P^{\′}}+z_Q^{\′}\·(x-l)\]}{R_2}$

${\mathrm{\Δ}}_g=\frac{\sqrt{\mathrm{\ϵ}-1+j60\mathrm{\λ}\mathrm{\σ}}}{\mathrm{\ϵ}+j60\mathrm{\λ}\mathrm{\σ}}$

其中，O为发射点的位置，P为观察点的位置，地面上动点的位置为Q，R₀表示从源点O到观察点P的直线距离；R₁和R₂分别表示地面上动点Q到源点O和观察点P之间的直线距离；z为地形高度，由步骤3提取的h_i得到；σ为动点Q处的电导率，根据AB_i判断得出；

其中I₁、I₂和I₃的计算过程如下：

采用高斯积分公式求解I₁，令则

$\begin{array}{c}{I}_1=\underset{-1}{\overset{1}{\∫}}{W}_{g}F\frac{\sqrt{x_1}}{\sqrt{x-\frac{1}{4}{x}_1{(1+y)}^2}}dy\\ =\sqrt{x_1}\underset{K=1}{\overset{5}{\Σ}}{\mathrm{\λ}}_K^{\left(5\right)}{W}_g\left(l_K^{\left(5\right)}\right)F\left(l_K^{\left(5\right)}\right)\frac{1}{\sqrt{x-l_K^{\left(5\right)}}}\end{array}$

和的取值可以通过数学手册查出；

根据n的不同取值，分别采用辛卜生和梯形公式计算I₂；

当n为偶数时

$I_2=\frac{\mathrm{\Δ}x}{3}\{\frac{W_g\left(x_1\right)F(x_n,x_1)}{\sqrt{x_1(x_n-x_1)}}+\underset{m=2}{\overset{n-2}{\Σ}}\frac{4}{M}\frac{W_g\left(x_m\right)F(x_n,x_m)}{\sqrt{x_m(x_n-x_m)}}+\frac{W_g\left(x_{n-1}\right)F(x_n,x_{n-1})}{\sqrt{x_{n-1}(x_n-x_{n-1})}}\}$

其中，

当n为奇数时

$\begin{array}{c}{I}_2=\underset{x_1}{\overset{x_{n-1}}{\∫}}{W}_g\left(Q\right)\frac{F(x,l)}{\sqrt{l(x-l)}}dl=\underset{x_1}{\overset{x_{n-2}}{\∫}}{W}_g\left(Q\right)\frac{F(x,l)}{\sqrt{l(x-l)}}dl+\underset{x_{n-2}}{\overset{x_{n-1}}{\∫}}{W}_g\left(Q\right)\frac{F(x,l)}{\sqrt{l(x-l)}}dl\\ =I_2^{\′}+I_2^{\′\′}\end{array}$

其中，I′₂仍可以用辛卜生公式求解，而I″₂可以用梯形公式进行求解，即

$I_2^{\′\′}=\frac{\mathrm{\Δ}x}{2}\{\frac{W_g\left(x_{n-2}\right)F(x_n,x_{n-2})}{\sqrt{x_{n-2}(x_n-x_{n-2})}}+\frac{W_g\left(x_{n-1}\right)F(x_n,x_{n-1})}{\sqrt{x_{n-1}(x_n-x_{n-1})}}\};$

由于I₃中不仅包含奇点l＝x，还包含待求结果W_g(x)，因此要先对其进行二次多项式拟合，而后再采用高斯积分公式，即

令对f(l)进行二次多项式拟合：

$\begin{array}{c}f\left(l\right)=\frac{+f\left(x_n\right)(l-x_{n-1})(l-x_{n-2})}{2{\mathrm{\Δx}}^2}\\ \frac{-2f\left(x_{n-1}\right)(l-x_n)(l-x_{n-2})}{2{\mathrm{\Δx}}^2}\\ \frac{+f\left(x_{n-2}\right)(l-x_n)(l-x_{n-1})}{2{\mathrm{\Δx}}^2}\end{array}$

再采用高斯积分公式，则

$I_3=\frac{-2}{15}\sqrt{\mathrm{\Δ}x}{W}_g\left(x_{n-2}\right)f\left(x_{n-2}\right)+\frac{14}{15}\sqrt{\mathrm{\Δ}x}{W}_g\left(x_{n-1}\right)f\left(x_{n-1}\right)+\frac{6}{5}\sqrt{\mathrm{\Δ}x}{W}_g\left(x_n\right)f\left(x_n\right);$

重复上述过程，得到所有观测点对应任意模型时的时延理论预测数据T_C：

步骤7.2、计算目标函数η

首先，根据步骤7.1所计算得到的时延理论预测数据T_C及各观测点时延测量数据计算时延理论预测值与时延实际测量值之间的偏差e：

其中，第i个观测点对应第k个模型时的偏差即为时延理论预测值与时延实际测量值的差值，即

$e_k^i=T_{C_k}^i-T_M^i$

其次，计算当前种群Ab|_t下的目标函数η(Ab|_t)

$\mathrm{\η}\left(Ab{|}_t\right)=\left[\begin{array}{c}{\mathrm{\η}}_1\\ {\mathrm{\η}}_2\\ .\\ .\\ .\\ {\mathrm{\η}}_S\end{array}\right]=\left[\begin{array}{c}\sqrt{\frac{{\left(e_1^1\right)}^2+{\left(e_1^2\right)}^2+\mathrm{...}+{\left(e_1^m\right)}^2}{m}}\\ \sqrt{\frac{{\left(e_2^1\right)}^2+{\left(e_2^2\right)}^2+\mathrm{...}+{\left(e_2^m\right)}^2}{m}}\\ .\\ .\\ .\\ \sqrt{\frac{{\left(e_S^1\right)}^2+{\left(e_S^2\right)}^2+\mathrm{...}+{\left(e_S^m\right)}^2}{m}}\end{array}\right]=\left[\begin{array}{c}\sqrt{\frac{\underset{i=1}{\overset{m}{\Σ}}{\left(e_1^i\right)}^2}{m}}\\ \sqrt{\frac{\underset{i=1}{\overset{m}{\Σ}}{\left(e_2^i\right)}^2}{m}}\\ .\\ .\\ .\\ \sqrt{\frac{\underset{i=1}{\overset{m}{\Σ}}{\left(e_S^i\right)}^2}{m}}\end{array}\right]$

最后，计算当前种群下目标函数η的最大值max(η)和均值mean(η)，并将当前种群下最优适应度值Fit_best对应的目标函数和个体分别作为最优目标函数η_best和最优个体σ_best。

步骤8判断条件为：最优目标函数η_best小于阈值ξ或达到最大迭代次数P，即

η_best＜ξ||t＞P。

步骤9对当前种群Ab|_t进行遗传算子操作，具体步骤如下：

步骤9.1、选择运算

首先，根据适应度越大选择概率越大的原则，定义选择概率P_s为

$P_s\left({\mathrm{\η}}_k\right)=\frac{Fit\left({\mathrm{\η}}_k\right)}{\underset{k=1}{\overset{S}{\Σ}}\[Fit\left({\mathrm{\η}}_k\right)\]};$

其次，产生S个[0,1]之间的均匀随机数，并与S个个体的选择概率进行比较，概率大的个体则被选中；

最后，重复选择过程，直至选出S个个体作为父本种群Ab_s|_t；

步骤9.2、交叉运算

对父本种群Ab_s|_t，采用随机选择的方法使种群中的个体两两配对，组成S/2对双亲，采用算数交叉以概率P_c产生子代Ab_c|_t；

以第i对双亲和为例，算数交叉产生子代的具体过程如下：

首先，产生一组随机数r＝[r₁r₂…r_n],r_l∈[0,1]；

其次，将n个随机数分别与交叉概率P_c进行比较，若第l个随机数r_l＞P_c，则将这两个父本中的第l个参数和直接传给两个子代的对应参数否则，两个父代参数按下式进行交叉，产生两个子代的对应参数和

$\left\{\begin{array}{c}{\mathrm{\σ}}_{il}^{C_1}=r_l\×{\mathrm{\σ}}_{il}^F+(1-r_l)\×{\mathrm{\σ}}_{il}^M\\ {\mathrm{\σ}}_{il}^{C_2}=r_l\×{\mathrm{\σ}}_{il}^M+(1-r_l)\×{\mathrm{\σ}}_{il}^F\end{array}\right.$

最后，按上述方法对每个参数进行运算，从而产生子代个体和

步骤9.3、变异运算

采用均匀变异的方法依照变异概率P_m对交叉运算后的种群Ab_c|_t进行变异，得到变异后的种群Ab_M|_t；

以第i个个体为例，其变异过程如下：

首先，产生一组随机数r＝[r₁r₂…r_n],r_l∈[0,1]；

其次，将n个随机数分别与变异概率P_m进行比较，以确定变异点。若第l个随机数r_l＜P_m，则第l点为变异点，将该点的参数按照下式进行变异，从而得到

${\mathrm{\σ}}_{il}^M={\mathrm{\σ}}_l^{LB}+r_l\×({\mathrm{\σ}}_l^{UB}-{\mathrm{\σ}}_l^{LB})$

最后，形成新个体

${\mathrm{\σ}}_i^M=\left[\begin{array}{cccc}{\mathrm{\σ}}_{i1}^M& {\mathrm{\σ}}_{i2}^M& \mathrm{...}& {\mathrm{\σ}}_{in}^M\end{array}\right];$

最终，通过上述遗传算子操作，产生新的种群Ab|_t+1

步骤10中回传到CPU的反演结果包括：最优大地电导率模型σ_best＝[σ₁σ₂…σ_n]即反演得到的各段大地电导率值、最优目标函数η_best即基于σ_best时的所有观测点的时延理论预测值与时延实测值之间偏差的方均根值、迭代次数t。

本发明的有益效果是，

1.本发明不规则地形下低频段大地电导率快速反演方法，不仅适用于不规则电波传播路径，同样也适用于光滑、平坦路径。与传统的以Millington经验公式为正演算法拟合路径上的大地电导率的方法相比，本发明方法能够有效剥离传统反演方法所得到的大地电导率(也是现有的大地电导率数据)中的地形影响等效成分；所获取的大地电导率数据能够更真实地反映地面的实际导电性能；便于应用于现有的更精确的低频地波场强/时延预测方法中而得到更精确的预测结果；同时为“数字地球”数据库的进一步完善打下基础。

2.本发明不规则地形下低频段大地电导率快速反演方法，是基于地波场强/时延法来进行大地电导率的反演，通过构建“实测+正演+反演”的综合反演平台，利用长波导航台/授时台获取观测点时延实测数据，以能够考虑地形起伏变化影响的积分方程方法为正演算法，结合遗传算法为优化算法，并利用基于GPU的并行技术，获取传播路径上的大地电导率分布信息。积分方程方法能够保证反演结果更接近实际地表导电性能，不包含地形影响等效成分；遗传算法能够保证反演结果真实可靠，不早熟且不陷入局部最优；GPU并行技术的引入，能够大大提高反演速度与效率，便于实际工程应用。

附图说明

图1是不规则路径下低频段大地电导率反演架构图；

图2是测量区域确定示意图；

图3是长波传播时延测量系统框图；

图4是基于GPU的并行遗传算法反演流程图；

图5是积分方程方法原理图；

图6是传播路径模型示意图；

图7是本发明方法与传统方法反演结果比较图；

图8是独立运行20次最优目标值分布情况图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

本发明不规则地形下低频段大地电导率快速反演方法，理论基础及原理是：

根据低频地波传播的理论可知，无线电波传播的场强(E)/时延(T_g)是信号频率f、传播距离d、辐射功率P_t和传播路径上媒质电参数(电导率σ、介电常数ε)的复杂函数，即 其中，媒质分为大气和地面两部分。大气环境的影响主要表现为大气折射率的复杂变化；而地面部分，在低频段主要考虑电波传播路径上的大地电导率σ在水平方向的非均匀分布。

在不考虑大气环境的复杂变化时，当频率、辐射功率和传播距离确定之后，场强/时延仅仅是大地电导率σ的函数。通常，通过接收长波授时台(如我国的BPL长波授时台)、导航台(如罗兰C台，即我国的长河二号系统)、广播台的信号，测量得到接收点的场强E和发射台到接收点的传播时延T_g。而后利用E和T_g与大地电导率σ之间的关系，反演出传播路径上的大地电导率分布。此即实测场强/时延反演大地电导率方法。综上所述，实测场强/时延反演大地电导率方法涉及三部分：(1)正演算法；(2)反演优化算法；(3)实测数据获取。如图1所示。

如果用X_i＝[(d_1i,σ₁),(d_2i,σ₂),...,(d_Ni,σ_N)]表示第i个观测点所在传播路径上的模型参数，其中，σ＝[σ₁,σ₂,…,σ_N]表示传播路径上的大地电导率分布情况，N表示电导率分区个数，d＝[d_1i,d_2i,…,d_Ni]为该传播路径中每个电导率分区所占的传播距离；用符号F表示地波传播场强/时延与X_i之间的函数关系(即为正演函数)，则F(X_i)为观测点的理论预测值；用M_i表示场强/时延实测值。此时，理论预测值和实测值之间的偏差e(X_i)可以表示为：

e(X_i)＝M_i-F(X_i)(1)

这里，将所有观测点的偏差e(X_i)的方均根值作为总的拟合误差η(σ)：

$\mathrm{\η}\left(\mathrm{\σ}\right)=\sqrt{\frac{\underset{i=1}{\overset{m}{\Σ}}{\[e\left(X_i\right)\]}^2}{m}}---\left(2\right)$

式中，m为观测点的个数。

大地电导率反演的目即为：用一定的数学方法寻找符合真实的大地电导率分布情况的模型参数。即寻找参数向量σ*＝[σ₁,σ₂,…,σ_N]，使得理论计算值F(X_i)与观测值M_i达到η(σ)最小。

在上述过程中，首先是正演算法的确定(即F的确定)，它能够准确描述电波传播过程；其次是测量系统的构建，以便获取准确的实测数据(即M_i)；最后是优化算法的确定，以期实现快速准确的反演。因此，本发明构建了“正演+实测+优化”的大地电导率反演架构。首先，采用积分方程方法作为正演算法，与传统的Millington经验公式相比，该方法能够考虑传播路径上地形起伏变化对地波传播的影响，使得反演结果中不再包含地形影响等效成分；其次，采用长波传播时延测量系统获取场强和时延实测数据；最后，采用基于GPU的并行遗传算法作为优化方法来进行大地电导率的反演，在保证反演结果的可靠性的同时提高反演速度。本发明方法能用于不规则路径下低频段大地电导率的快速高精度反演。

本发明不规则地形下低频段大地电导率快速反演方法，具体分为三个阶段共10步，详细过程如下：

阶段一：获取基本信息，完成基础数据准备

步骤1、针对待反演区域，选定发射台(即后续测量系统所接收信号来源)，确定测量区域；

待反演区域指地球表层某一大地电导率取值未知且待反演的地理区域；

发射台的选定：选取距离待反演区域较近、信号质量较好的长波导航台或授时台为信号来源，用于后续进行接收、测量；

测量区域的确定：以所选定发射台为中心，做间隔0.01度、辐射距离为1000km的射线，则所有穿过待反演区域的射线中最外侧两条射线所夹不规则区域(如图2中阴影区域)即为测量区域。

步骤2、利用长波传播时延测量系统，获取各观测点的传播时延测量值T_M；

观测点的选取标准：在步骤1所限定的测量区域内，综合考虑实际地形及交通便利情况，合理规划路线，进行定点测量；测量点应远离山区或峡谷等地形起伏区域；测量点后向地形应较平坦。

长波传播时延测量系统如图3所示，主要包括：

1)长波接收机：提供从发射台到测量点的以GPS时钟为基准的到达时刻；

2)高精度GPS接收机：提供测量点的位置信息、时间信息、1pps信号

3)计数器：比对长波接收机与GPS接收机1pps信号的时间间隔；

4)工控机：监测长波接收机、GPS接收机、计数器的运行状态，采集和处理数据。

步骤3、利用先验知识，结合GIS技术获取各观测点所对应传播路径参数信息：包括电导率分区数目、各分区内电导率取值范围、各分区内的传播距离以及传播路径上的高程信息；

获取各观测点所对应传播路径参数信息过程如下(以第i个观测点为例)：

1)提取观测点所在传播路径的高程信息h_i(即地形起伏变化情况)；

2)结合高程信息和现有大地电导率地图进行大地电导率分区，并给出分区数目n和各分区内电导率取值范围(l表示第l个分区)；

3)计算传播路径上各段路径长度(大圆距离)d_ij(即第i个观测点第j段大圆距离)。

阶段二：并行反演(流程如图4所示)

步骤4、根据步骤3得到的传播路径参数信息，在CPU上初始化参数；

需初始化的参数主要有：积分步长Δx、大地电导率种群规模S、迭代次数P、目标函数阈值ξ、交叉概率P_c、变异概率P_m；大地电导率分区数目N、观测点个数m；各观测点时延测量数据各观测点所对应传播路径的高程数据h＝[h₁h₂…h_m]。

步骤5、产生大地电导率初始种群，形成各观测点所对应模型文件；

产生初始种群及建立观测点模型文件过程如下：

1)种群中第k个个体σ_k的第l个参数σ_kl在其取值范围内随机生成(式中r为随机数)；

${\mathrm{\σ}}_{kl}={\mathrm{\σ}}_l^{LB}+r\×({\mathrm{\σ}}_l^{UB}-{\mathrm{\σ}}_l^{LB})$

2)按照1)的方法，形成第k个个体σ_k＝[σ_k1σ_k2...σ_kn]，此即为一种大地电导率分布模型；

3)重复步骤1)和2)，得到S个个体(模型)，形成大地电导率初始种群Ab|₀

4)结合步骤3中得到的各观测点所在传播路径上各段路径长度d_ij，则第i个观测点的模型为

步骤6、将初始化的数据传递到GPU，完成设备初始化、数据准备，启动核函数；

步骤7、计算大地电导率初始种群适应度；

定义第k个个体σ_k的适应度为

$Fit\left({\mathrm{\η}}_k\right)=\frac{max\left(\mathrm{\η}\right)-{\mathrm{\η}}_k}{S\×\[max\left(\mathrm{\η}\right)-mean\left(\mathrm{\η}\right)\]}$

式中，η_k为σ_k的目标函数，max(η)和mean(η)分别为当前种群下所有个体所对应目标函数的最大值和均值。

则目标函数η_k的计算过程分为如下步骤：

步骤7.1、采用积分方程方法并行计算第t代种群Ab|_t下每一个观测点在每一种模型下的时延预测值T_C。

则第i个观测点对应第k个模型时的时延预测值为

$T_{C_k}^i=\frac{\arg \left(W_g\right)}{\mathrm{\ω}}\×{10}^6\mathrm{\μ}s.$

式中，

$\begin{array}{c}{W}_g\left(P\right)=1-R_0\sqrt{\frac{1}{j\mathrm{\λ}x}}\{\underset{0}{\overset{x_1}{\∫}}{W}_g\left(Q\right)\frac{F(x,l)}{\sqrt{l(x-l)}}dl+\underset{x_1}{\overset{x_{n-1}}{\∫}}{W}_g\left(Q\right)\frac{F(x,l)}{\sqrt{l(x-l)}}dl+\underset{x_{n-1}}{\overset{x_n}{\∫}}{W}_g\left(Q\right)\frac{F(x,l)}{\sqrt{l(x-l)}}dl\}\\ =1-R_0\sqrt{\frac{1}{j\mathrm{\λ}x}}(I_1+I_2+I_3)\end{array}$

$F(x,l)=\frac{l(x-l)}{R_1{R}_2}{e}^{{\mathrm{jk}}_0(R_1+R_2-R_0)}\[{\mathrm{\Δ}}_g+(1+\frac{j}{k_0{R}_2})\]\frac{\[z_Q-z_{P^{\′}}+z_Q^{\′}\·(x-l)\]}{R_2}$

${\mathrm{\Δ}}_g=\frac{\sqrt{\mathrm{\ϵ}-1+j60\mathrm{\λ}\mathrm{\σ}}}{\mathrm{\ϵ}+j60\mathrm{\λ}\mathrm{\σ}}$

如图5所示，O为发射点的位置，P为观察点的位置，地面上动点的位置为Q，R₀表示从源点O到观察点P的直线距离；R₁和R₂分别表示地面上动点Q到源点O和观察点P之间的直线距离；z为地形高度，由步骤3提取的h_i得到；σ为动点Q处的电导率，根据AB_i判断得出。

其中I₁、I₂和I₃的计算过程如下：

采用高斯积分公式求解I₁，令则

$\begin{array}{c}{I}_1=\underset{-1}{\overset{1}{\∫}}{W}_{g}F\frac{\sqrt{x_1}}{\sqrt{x-\frac{1}{4}{x}_1{(1+y)}^2}}dy\\ =\sqrt{x_1}\underset{K=1}{\overset{5}{\Σ}}{\mathrm{\λ}}_K^{\left(5\right)}{W}_g\left(l_K^{\left(5\right)}\right)F\left(l_K^{\left(5\right)}\right)\frac{1}{\sqrt{x-l_K^{\left(5\right)}}}\end{array}$

和的取值可以通过数学手册查出。

根据n的不同取值，分别采用辛卜生和梯形公式计算I₂，

当n为偶数时

$I_2=\frac{\mathrm{\Δ}x}{3}\{\frac{W_g\left(x_1\right)F(x_n,x_1)}{\sqrt{x_1(x_n-x_1)}}+\underset{m=2}{\overset{n-2}{\Σ}}\frac{4}{M}\frac{W_g\left(x_m\right)F(x_n,x_m)}{\sqrt{x_m(x_n-x_m)}}+\frac{W_g\left(x_{n-1}\right)F(x_n,x_{n-1})}{\sqrt{x_{n-1}(x_n-x_{n-1})}}\}$

其中，

当n为奇数时

$\begin{array}{c}{I}_2=\underset{x_1}{\overset{x_{n-1}}{\∫}}{W}_g\left(Q\right)\frac{F(x,l)}{\sqrt{l(x-l)}}dl=\underset{x_1}{\overset{x_{n-2}}{\∫}}{W}_g\left(Q\right)\frac{F(x,l)}{\sqrt{l(x-l)}}dl+\underset{x_{n-2}}{\overset{x_{n-1}}{\∫}}{W}_g\left(Q\right)\frac{F(x,l)}{\sqrt{l(x-l)}}dl\\ =I_2^{\′}+I_2^{\′\′}\end{array}$

其中，I′₂仍可以用辛卜生公式求解，而I″₂可以用梯形公式进行求解，即

$I_2^{\′\′}=\frac{\mathrm{\Δ}x}{2}\{\frac{W_g\left(x_{n-2}\right)F(x_n,x_{n-2})}{\sqrt{x_{n-2}(x_n-x_{n-2})}}+\frac{W_g\left(x_{n-1}\right)F(x_n,x_{n-1})}{\sqrt{x_{n-1}(x_n-x_{n-1})}}\}$

由于I₃中不仅包含奇点l＝x，还包含待求结果W_g(x)，因此要先对其进行二次多项式拟合，而后再采用高斯积分公式，即

令对f(l)进行二次多项式拟合：

$\begin{array}{c}f\left(l\right)=\frac{+f\left(x_n\right)(l-x_{n-1})(l-x_{n-2})}{2{\mathrm{\Δx}}^2}\\ \frac{-2f\left(x_{n-1}\right)(l-x_n)(l-x_{n-2})}{2{\mathrm{\Δx}}^2}\\ \frac{+f\left(x_{n-2}\right)(l-x_n)(l-x_{n-1})}{2{\mathrm{\Δx}}^2}\end{array}$

再采用高斯积分公式，则

$I_3=\frac{-2}{15}\sqrt{\mathrm{\Δ}x}{W}_g\left(x_{n-2}\right)f\left(x_{n-2}\right)+\frac{14}{15}\sqrt{\mathrm{\Δ}x}{W}_g\left(x_{n-1}\right)f\left(x_{n-1}\right)+\frac{6}{5}\sqrt{\mathrm{\Δ}x}{W}_g\left(x_n\right)f\left(x_n\right).$

重复上述过程，得到所有观测点对应任意模型时的时延理论预测数据T_C：

步骤7.2、计算目标函数η。

首先，根据步骤7.1所计算得到的时延理论预测数据T_C及各观测点时延测量数据计算时延理论预测值与时延实际测量值之间的偏差e：

其中，第i个观测点对应第k个模型时的偏差即为时延理论预测值与时延实际测量值的差值，即

$e_k^i=T_{C_k}^i-T_M^i$

其次，计算当前种群Ab|_t下的目标函数η(Ab|_t)

$\mathrm{\η}\left(Ab{|}_t\right)=\left[\begin{array}{c}{\mathrm{\η}}_1\\ {\mathrm{\η}}_2\\ .\\ .\\ .\\ {\mathrm{\η}}_S\end{array}\right]=\left[\begin{array}{c}\sqrt{\frac{{\left(e_1^1\right)}^2+{\left(e_1^2\right)}^2+\mathrm{...}+{\left(e_1^m\right)}^2}{m}}\\ \sqrt{\frac{{\left(e_2^1\right)}^2+{\left(e_2^2\right)}^2+\mathrm{...}+{\left(e_2^m\right)}^2}{m}}\\ .\\ .\\ .\\ \sqrt{\frac{{\left(e_S^1\right)}^2+{\left(e_S^2\right)}^2+\mathrm{...}+{\left(e_S^m\right)}^2}{m}}\end{array}\right]=\left[\begin{array}{c}\sqrt{\frac{\underset{i=1}{\overset{m}{\Σ}}{\left(e_1^i\right)}^2}{m}}\\ \sqrt{\frac{\underset{i=1}{\overset{m}{\Σ}}{\left(e_2^i\right)}^2}{m}}\\ .\\ .\\ .\\ \sqrt{\frac{\underset{i=1}{\overset{m}{\Σ}}{\left(e_S^i\right)}^2}{m}}\end{array}\right]$

最后，计算当前种群下目标函数η的最大值max(η)和均值mean(η)，并将当前种群下最优适应度值Fit_best对应的目标函数和个体分别作为最优目标函数η_best和最优个体σ_best。

步骤8、判断是否达到结束条件，若达到，执行步骤10，若未达到，执行步骤9后返回步骤7；

判断条件为：最优目标函数η_best小于阈值ξ或达到最大迭代次数P，即

η_best＜ξ||t＞P。

步骤9、遗传算子操作，产生新个体，形成新的大地电导率种群；

对当前种群Ab|_t进行遗传算子操作，具体步骤如下：

步骤9.1、选择运算

首先，根据适应度越大选择概率越大的原则，定义选择概率P_s为

$P_s\left({\mathrm{\η}}_k\right)=\frac{Fit\left({\mathrm{\η}}_k\right)}{\underset{k=1}{\overset{S}{\Σ}}\[Fit\left({\mathrm{\η}}_k\right)\]}.$

其次，产生S个[0,1]之间的均匀随机数，并与S个个体的选择概率进行比较，概率大的个体则被选中。

最后，重复选择过程，直至选出S个个体作为父本种群Ab_s|_t。

步骤9.2、交叉运算

对父本种群Ab_s|_t，采用随机选择的方法使种群中的个体两两配对，组成S/2对双亲，采用算数交叉以概率P_c产生子代Ab_c|_t。

以第i对双亲和为例，算数交叉产生子代的具体过程如下：

首先，产生一组随机数r＝[r₁r₂…r_n],r_l∈[0,1]。

其次，将n个随机数分别与交叉概率P_c进行比较，若第l个随机数r_l＞P_c，则将这两个父本中的第l个参数和直接传给两个子代的对应参数否则，两个父代参数按下式进行交叉，产生两个子代的对应参数和

$\left\{\begin{array}{c}{\mathrm{\σ}}_{il}^{C_1}=r_l\×{\mathrm{\σ}}_{il}^F+(1-r_l)\×{\mathrm{\σ}}_{il}^M\\ {\mathrm{\σ}}_{il}^{C_2}=r_l\×{\mathrm{\σ}}_{il}^M+(1-r_l)\×{\mathrm{\σ}}_{il}^F\end{array}\right.$

最后，按上述方法对每个参数进行运算，从而产生子代个体和

步骤9.3、变异运算

采用均匀变异的方法依照变异概率P_m对交叉运算后的种群Ab_c|_t进行变异，得到变异后的种群Ab_M|_t。

以第i个个体为例，其变异过程如下：

首先，产生一组随机数r＝[r₁r₂…r_n],r_l∈[0,1]。

其次，将n个随机数分别与变异概率P_m进行比较，以确定变异点。若第l个随机数r_l＜P_m，则第l点为变异点，将该点的参数按照下式进行变异，从而得到

${\mathrm{\σ}}_{il}^M={\mathrm{\σ}}_l^{LB}+r_l\×({\mathrm{\σ}}_l^{UB}-{\mathrm{\σ}}_l^{LB})$

最后，形成新个体

${\mathrm{\σ}}_i^M=\left[\begin{array}{cccc}{\mathrm{\σ}}_{i1}^M& {\mathrm{\σ}}_{i2}^M& \mathrm{...}& {\mathrm{\σ}}_{in}^M\end{array}\right].$

最终，通过上述遗传算子操作，产生新的种群Ab|_t+1

阶段三：得到反演结果。

步骤10、将反演结果传回CPU。

回传到CPU的反演结果包括：最优大地电导率模型σ_best＝[σ₁σ₂…σ_n](即反演得到的各段大地电导率值)、最优目标函数η_best(即基于σ_best时的所有观测点的时延理论预测值与时延实测值之间偏差的方均根值)、迭代次数t。

实施例1

本发明方法有效性及可靠性验证。

1)仿真验证模型

建立如图6所示的地面模型，其中，传播路径总长度为100km；整个路径分为三段，第一段和第三段的电导率分别为σ₁＝σ₃＝3×10^-3S/m，第二段为σ₂＝1×10^-3S/m；地形起伏为高斯曲线取山的峰值高度H＝2km，山的宽度l＝20km，山的中心位置与发射天线之间的距离ρ₀＝50km。

取传播路径上的四个点(分别是20km、40km、60km和80km处)，以FDTD计算得到的场强值作为实测值。

2)仿真验证方法

分别采用本发明方法和以Millington经验公式为正演算法(即在本发明方法中适应度计算时采用Millington经验公式计算时延理论预测值)进行仿真模型中大地电导率的反演，将反演结果与实际模型进行比较。其中遗传算法的结束条件为：η_best＜0.001||t＞1000(即达到最大迭代次数1000或达到最优目标函数小于阈值0.001μs(即所有观测点的时延理论预测值与时延实测值之间的偏差的方均根值小于1纳秒))。

3)仿真结果

将两种反演方法各独立运行20次，反演结果如图7和图8所示。由反演结果可见，与传统的基于Millington经验公式的大地电导率拟合方法相比，本发明方法可以得到更接近仿真模型的大地电导率分布数据，有效剔除传统反演结果中的地形影响等效成分，具有较高的精度。

实施例2

本发明方法并行加速比验证。

仿真模型如实例1中模型，分别采用基于GPU并行和基于CPU串行的算法程序进行大地电导率反演，将二者的运行时间和加速比进行比较，以验证本发明方法的加速性能。其中仿真实验环境为：CPU为Inter(R)Core(TM)2DuoCPUE7300，显卡为NVIDIAGeForce9800GT。仿真结果如表1所示。

表1基于GPU并行和基于CPU串行的算法运行时间和加速比

仿真结果显示，随着种群规模的增加，本发明方法反演速度更快、效率更高，能够很好的解决大规模低频段大地电导率反演问题。

Claims (10)

1.一种不规则地形下低频段大地电导率快速反演方法，其特征在于，具体步骤如下：

步骤1、针对待反演区域，选定发射台即后续测量系统所接收信号来源，确定测量区域；

步骤2、利用长波传播时延测量系统，获取各观测点的传播时延测量值T_M；

步骤3、利用先验知识，结合GIS技术获取各观测点所对应传播路径参数信息：包括电导率分区数目、各分区内电导率取值范围、各分区内的传播距离以及传播路径上的高程信息；

步骤4、根据步骤3得到的传播路径参数信息，在CPU上初始化参数；

步骤5、产生大地电导率初始种群，形成各观测点所对应模型文件；

步骤6、将步骤4初始化参数传递到GPU，完成设备初始化、数据准备，启动核函数；

步骤7、计算大地电导率初始种群适应度；

步骤8、判断是否达到结束条件，若达到，执行步骤10，若未达到，执行步骤9后返回步骤7；

步骤9、遗传算子操作，产生新个体，形成新的大地电导率种群；

步骤10、将反演结果传回CPU。

2.根据权利要求1所述的不规则地形下低频段大地电导率快速反演方法，其特征在于，步骤1中待反演区域指地球表层某一大地电导率取值未知且待反演的地理区域；

发射台的选定：选取距离待反演区域较近、信号质量较好的长波导航台或授时台为信号来源，用于后续进行接收、测量；

测量区域的确定：以所选定发射台为中心，做间隔0.01度、辐射距离为1000km的射线，则所有穿过待反演区域的射线中最外侧两条射线所夹不规则区域即为测量区域。

3.根据权利要求1中所述的不规则地形下低频段大地电导率快速反演方法，其特征在于，步骤2中观测点的选取标准：在步骤1所限定的测量区域内，综合考虑实际地形及交通便利情况，合理规划路线，进行定点测量；测量点应远离山区或峡谷等地形起伏区域；测量点后向地形应较平坦；

长波传播时延测量系统主要包括：

1)长波接收机：提供从发射台到测量点的以GPS时钟为基准的到达时刻；

2)高精度GPS接收机：提供测量点的位置信息、时间信息、1pps信号；

3)计数器：比对长波接收机与GPS接收机1pps信号的时间间隔；

4)工控机：监测长波接收机、GPS接收机、计数器的运行状态，采集和处理数据。

4.根据权利要求1所述的不规则地形下低频段大地电导率快速反演方法，其特征在于，步骤3中获取各观测点所对应传播路径参数信息过程如下，以第i个观测点为例：

1)提取观测点所在传播路径的高程信息h_i即地形起伏变化情况；

2)结合高程信息和现有大地电导率地图进行大地电导率分区，并给出分区数目N和各分区内电导率取值范围l表示第l个分区；

3)计算传播路径上各段路径长度d_ij，即第i个观测点第j段大圆距离。

5.根据权利要求1所述的不规则地形下低频段大地电导率快速反演方法，其特征在于步骤4中初始化的参数有：积分步长Δx、大地电导率种群规模S、迭代次数P、目标函数阈值ξ、交叉概率P_c、变异概率P_m；大地电导率分区数目n、观测点个数m；各观测点时延测量数据各观测点所对应传播路径的高程数据h＝[h₁h₂…h_m]。

6.根据权利要求4所述的不规则地形下低频段大地电导率快速反演方法，其特征在于，步骤5中产生大地电导率初始种群及建立观测点模型文件过程如下：

1)种群中第k个个体σ_k的第l个参数σ_kl在其取值范围内随机生成，式中r为随机数；

${\mathrm{\σ}}_{kl}={\mathrm{\σ}}_l^{LB}+r\×({\mathrm{\σ}}_l^{UB}-{\mathrm{\σ}}_l^{LB})$

2)按照1)的方法，形成第k个个体σ_k＝[σ_k1σ_k2...σ_kn]，此即为一种大地电导率分布模型；

3)重复步骤1)和2)，得到S个个体模型，形成大地电导率初始种群Ab|₀

4)结合步骤3中得到的各观测点所在传播路径上各段路径长度d_ij，则第i个观测点的模型为：

7.根据权利要求1所述的不规则地形下低频段大地电导率快速反演方法，其特征在于，步骤7中定义第k个个体σ_k的适应度为

$Fit\left({\mathrm{\η}}_k\right)=\frac{max\left(\mathrm{\η}\right)-{\mathrm{\η}}_k}{S\×\[\max \left(\mathrm{\η}\right)-mean\left(\mathrm{\η}\right)\]}$

式中，η_k为σ_k的目标函数，max(η)和mean(η)分别为当前种群下所有个体所对应目标函数的最大值和均值；

则目标函数η_k的计算过程分为如下步骤：

步骤7.1、采用积分方程方法并行计算第t代种群Ab|_t下每一个观测点在每一种模型下的时延预测值T_C；

则第i个观测点对应第k个模型时的时延预测值为

$T_{C_k}^i=\frac{\arg \left(W_g\right)}{\mathrm{\ω}}\×{10}^6\mathrm{\μ}s;$

式中，

$\begin{array}{c}{W}_g\left(P\right)=1-R_0\sqrt{\frac{1}{j\mathrm{\λ}x}}\{\underset{0}{\overset{x_1}{\∫}}{W}_g\left(Q\right)\frac{F(x,l)}{\sqrt{l(x-l)}}dl+\underset{x_1}{\overset{x_{n-1}}{\∫}}{W}_g\left(Q\right)\frac{F(x,l)}{\sqrt{l(x-l)}}dl+\underset{x_{n-1}}{\overset{x_n}{\∫}}{W}_g\left(Q\right)\frac{F(x,l)}{\sqrt{l(x-l)}}dl\}\\ =1-R_0\sqrt{\frac{1}{j\mathrm{\λ}x}}(I_1+I_2+I_3)\end{array}$

$F(x,l)=\frac{l(x-l)}{R_1{R}_2}{w}^{{\mathrm{jk}}_0(R_1+R_2-R_0)}\[{\mathrm{\Δ}}_g+(1+\frac{j}{k_0{R}_0})\]\frac{\[z_Q-z_{P^{\′}}+z_Q^{\′}\·(x-l)\]}{R_2}$

${\mathrm{\Δ}}_g=\frac{\sqrt{\mathrm{\ϵ}-1+j60\mathrm{\λ}\mathrm{\σ}}}{\mathrm{\ϵ}+j60\mathrm{\λ}\mathrm{\σ}}$

其中，O为发射点的位置，P为观察点的位置，地面上动点的位置为Q，R₀表示从源点O到观察点P的直线距离；R₁和R₂分别表示地面上动点Q到源点O和观察点P之间的直线距离；z为地形高度，由步骤3提取的h_i得到；σ为动点Q处的电导率，根据AB_i判断得出；

其中I₁、I₂和I₃的计算过程如下：

采用高斯积分公式求解I₁，令则

$\begin{array}{c}{I}_1=\underset{-1}{\overset{1}{\∫}}{W}_{g}F\frac{\sqrt{x_1}}{\sqrt{x-\frac{1}{4}{x}_1{(1+y)}^2}}dy\\ =\sqrt{x_1}\underset{K=1}{\overset{5}{\Σ}}{\mathrm{\λ}}_K^{\left(5\right)}{W}_g\left(l_K^{\left(5\right)}\right)F\left(l_K^{\left(5\right)}\right)\frac{1}{\sqrt{x-l_K^{\left(5\right)}}}\end{array}$

和的取值可以通过数学手册查出；

根据n的不同取值，分别采用辛卜生和梯形公式计算I₂；

当n为偶数时

$I_2=\frac{\mathrm{\Δ}x}{3}\left\{\frac{W_g\left(x_1\right)F(x_n,x_1)}{\sqrt{x_1(x_n-x_1)}}+\underset{m=2}{\overset{n-2}{\Σ}}\frac{4}{M}\frac{W_g\left(x_m\right)F(x_n,x_m)}{\sqrt{x_m(x_n-x_m)}}+\frac{W_g\left(x_{n-1}\right)F(x_n,x_{n-1})}{\sqrt{x_{n-1}(x_n-x_{n-1})}}\right\}$

其中，

当n为奇数时

$\begin{array}{c}{I}_2=\underset{x_1}{\overset{x_{n-1}}{\∫}}{W}_g\left(Q\right)\frac{F(x,l)}{\sqrt{l(x-l)}}dl=\underset{x_1}{\overset{x_{n-2}}{\∫}}{W}_g\left(Q\right)\frac{F(x,l)}{\sqrt{l(x-l)}}dl+\underset{x_{n-2}}{\overset{x_{n-1}}{\∫}}{W}_g\left(Q\right)\frac{F(x,l)}{\sqrt{l(x-l)}}dl\\ =I_2^{\′}+I_2^{\′\′}\end{array}$

其中，I′₂仍可以用辛卜生公式求解，而I″₂可以用梯形公式进行求解，即

$I_2^{\′\′}=\frac{\mathrm{\Δ}x}{2}\left\{\frac{W_g\left(x_{n-2}\right)F(x_n,x_{n-2})}{\sqrt{x_{n-2}(x_n-x_{n-2})}}+\frac{W_g\left(x_{n-1}\right)F(x_n,x_{n-1})}{\sqrt{x_{n-1}(x_n-x_{n-1})}}\right\};$

由于I₃中不仅包含奇点l＝x，还包含待求结果W_g(x)，因此要先对其进行二次多项式拟合，而后再采用高斯积分公式，即

令对f(l)进行二次多项式拟合：

$\begin{array}{c}f\left(l\right)=\frac{+f\left(x_n\right)(l-x_{n-1})(l-x_{n-2})}{2{\mathrm{\Δx}}^2}\\ \frac{-2f\left(x_{n-1}\right)(l-x_n)(l-x_{n-2})}{2{\mathrm{\Δx}}^2}\\ \frac{+f\left(x_{n-2}\right)(l-x_n)(l-x_{n-1})}{2{\mathrm{\Δx}}^2}\end{array}$

再采用高斯积分公式，则

$I_3=\frac{-2}{15}\sqrt{\mathrm{\Δ}x}{W}_g\left(x_{n-2}\right)f\left(x_{n-2}\right)+\frac{14}{15}\sqrt{\mathrm{\Δ}x}{W}_g\left(x_{n-1}\right)f\left(x_{n-1}\right)+\frac{6}{5}\sqrt{\mathrm{\Δ}x}{W}_g\left(x_n\right)f\left(x_n\right);$

重复上述过程，得到所有观测点对应任意模型时的时延理论预测数据T_C：

步骤7.2、计算目标函数η

首先，根据步骤7.1所计算得到的时延理论预测数据T_C及各观测点时延测量数据计算时延理论预测值与时延实际测量值之间的偏差e：

其中，第i个观测点对应第k个模型时的偏差即为时延理论预测值与时延实际测量值的差值，即

$e_k^i=T_{C_k}^i-T_M^i$

其次，计算当前种群Ab|_t下的目标函数η(Ab|_t)

$\mathrm{\η}\left(Ab{|}_t\right)=\left[\begin{array}{c}{\mathrm{\η}}_1\\ {\mathrm{\η}}_2\\ .\\ .\\ .\\ {\mathrm{\η}}_S\end{array}\right]=\left[\begin{array}{c}\sqrt{\frac{{\left(e_1^1\right)}^2+{\left(e_1^2\right)}^2+\mathrm{...}+{\left(e_1^m\right)}^2}{m}}\\ \sqrt{\frac{{\left(e_2^1\right)}^2+{\left(e_2^2\right)}^2+\mathrm{...}+{\left(e_2^m\right)}^2}{m}}\\ .\\ .\\ .\\ \sqrt{\frac{{\left(e_S^1\right)}^2+{\left(e_S^2\right)}^2+\mathrm{...}+{\left(e_S^m\right)}^2}{m}}\end{array}\right]=\left[\begin{array}{c}\sqrt{\frac{\underset{i=1}{\overset{m}{\Σ}}{\left(e_1^i\right)}^2}{m}}\\ \sqrt{\frac{\underset{i=1}{\overset{m}{\Σ}}{\left(e_2^i\right)}^2}{m}}\\ .\\ .\\ .\\ \sqrt{\frac{\underset{i=1}{\overset{m}{\Σ}}{\left(e_S^i\right)}^2}{m}}\end{array}\right]$

最后，计算当前种群下目标函数η的最大值max(η)和均值mean(η)，并将当前种群下最优适应度值Fit_best对应的目标函数和个体分别作为最优目标函数η_best和最优个体σ_best。

8.根据权利要求1所述的不规则地形下低频段大地电导率快速反演方法，其特征在于，步骤8判断条件为：最优目标函数η_best小于阈值ξ或达到最大迭代次数P，即

η_best<ξ||t>P。

9.根据权利要求7所述的不规则地形下低频段大地电导率快速反演方法，其特征在于，步骤9对当前种群Ab|_t进行遗传算子操作，具体步骤如下：

步骤9.1、选择运算

首先，根据适应度越大选择概率越大的原则，定义选择概率P_s为

$P_s\left({\mathrm{\&eta;}}_k\right)=\frac{Fit\left({\mathrm{\&eta;}}_k\right)}{\underset{k=1}{\overset{S}{\&Sigma;}}\&lsqb;Fit\left({\mathrm{\&eta;}}_k\right)\&rsqb;};$

其次，产生S个[0,1]之间的均匀随机数，并与S个个体的选择概率进行比较，概率大的个体则被选中；

最后，重复选择过程，直至选出S个个体作为父本种群Ab_s|_t；

步骤9.2、交叉运算

对父本种群Ab_s|_t，采用随机选择的方法使种群中的个体两两配对，组成S/2对双亲，采用算数交叉以概率P_c产生子代Ab_c|_t；

以第i对双亲和为例，算数交叉产生子代的具体过程如下：

首先，产生一组随机数r＝[r₁r₂…r_n],r_l∈[0,1]；

其次，将n个随机数分别与交叉概率P_c进行比较，若第l个随机数r_l>P_c，则将这两个父本中的第l个参数和直接传给两个子代的对应参数否则，两个父代参数按下式进行交叉，产生两个子代的对应参数和

$\left\{\begin{array}{c}{\mathrm{\&sigma;}}_{il}^{C_1}=r_l\&times;{\mathrm{\&sigma;}}_{il}^F+(1-r_l)\&times;{\mathrm{\&sigma;}}_{il}^M\\ {\mathrm{\&sigma;}}_{il}^{C_2}=r_l\&times;{\mathrm{\&sigma;}}_{il}^M+(1-r_l)\&times;{\mathrm{\&sigma;}}_{il}^F\end{array}\right.$

最后，按上述方法对每个参数进行运算，从而产生子代个体和

步骤9.3、变异运算

采用均匀变异的方法依照变异概率P_m对交叉运算后的种群Ab_c|_t进行变异，得到变异后的种群Ab_M|_t；

以第i个个体为例，其变异过程如下：

首先，产生一组随机数r＝[r₁r₂…r_n],r_l∈[0,1]；

其次，将n个随机数分别与变异概率P_m进行比较，以确定变异点，若第l个随机数r_l<P_m，则第l点为变异点，将该点的参数按照下式进行变异，从而得到

${\mathrm{\&sigma;}}_{il}^M={\mathrm{\&sigma;}}_l^{LB}+r_l\&times;({\mathrm{\&sigma;}}_l^{UB}-{\mathrm{\&sigma;}}_l^{LB})$

最后，形成新个体

${\mathrm{\&sigma;}}_i^M=\left[\begin{array}{cccc}{\mathrm{\&sigma;}}_{i1}^M& {\mathrm{\&sigma;}}_{i2}^M& ...& {\mathrm{\&sigma;}}_{in}^M\end{array}\right];$

最终，通过上述遗传算子操作，产生新的种群Ab|_t+1

10.根据权利要求1所述的不规则地形下低频段大地电导率快速反演方法，其特征在于，步骤10中回传到CPU的反演结果包括：最优大地电导率模型σ_best＝[σ₁σ₂…σ_n]即反演得到的各段大地电导率值、最优目标函数η_best即基于σ_best时的所有观测点的时延理论预测值与时延实测值之间偏差的方均根值、迭代次数t。