CN112666612B - 基于禁忌搜索的大地电磁二维反演方法 - Google Patents

Abstract

本发明公开了一种基于禁忌搜索的大地电磁二维反演方法，在大地电磁反演迭代前，利用钻探、测井、地质调绘、区域地质图等先知信息创建电阻率禁忌搜索列表，提出禁忌搜索原则，减小电阻率参数搜索范围；根据禁忌搜索表构建反演初始模型，然后进行二维正演、残差计算，之后判断目标函数是否设定最小值，如达到，则输出地层电阻率、厚度值；如未达到，则根据灵敏度矩阵修改反演模型，不断迭代，直至目标函数达到设定最小值。在灵敏度矩阵计算中采用互换定理，将网格节点与测点互换，利用正演一次求解，获取多节点的灵敏度，减小了计算量，提高了反演结果的准确度和反演效率。

Description

基于禁忌搜索的大地电磁二维反演方法

技术领域

本发明属于铁路工程地球物理勘探技术领域，具体涉及一种基于禁忌搜索的大地电磁二维反演方法。

背景技术

目前，铁路隧道多采用大地电磁法进行地质勘察工作，但目前大地电磁对钻探、调绘资料的利用停留在反演初始模型建立和成果验证阶段，未将这些信息应用到反演迭代过程中；同时，反演灵敏度计算求解速度较慢，极大地影响了反演实施，因此需要提高灵敏度的计算速度。

发明内容

为解决上述现有技术中存在的问题，本发明提供一种能够提高灵敏度计算速度的基于禁忌搜索的大地电磁二维反演方法。

为了实现上述目的，本发明采用以下技术方案：

本发明的基于禁忌搜索的大地电磁二维反演方法，通过测井、电测深、调绘资料和分析研究，确定工区岩性分类、地下地层层序、层厚及各地层岩石电阻率，然后根据这些信息建立二维空间模型中各单元的禁忌搜索列表，列表中的参数包括：地层层序、地层厚度范围和地层电阻率范围，这些参数的下一级参数中包含二维空间模型的坐标信息。同时，建立禁忌搜索原则，逐步迭代，不断拟合实测视电阻率，其中，灵敏度矩阵计算中，采取互换源的方式，通过一次正演计算得到与一个测点对应的地下空间各节点的灵敏度，从而大幅提高反演效率。

该方法的步骤如下：

一种基于禁忌搜索的大地电磁二维反演方法，包括以下步骤：

S₁、构建二维空间模型：将地下二维空间从横、纵两个方向剖分为矩形网格，在观测目标区域内网格大小为等间隔，且网格大小与拟探测目标体大小相当；在观测区域外部，随着距离增大网格大小也逐渐线性增大；

S₂、搜集测量工区以往的钻探记录、测井成果以及地质调查资料；

S₃、依据步骤S₂搜集的资料进行综合分析，初步确定工区岩性分类、地下各地层层序、层厚度及各地层岩石电阻率值范围信息，根据这些信息建立二维空间模型中各单元的禁忌搜索列表，该禁忌搜索列表中的参数包括地层厚度范围和地层电阻率范围，这些参数的下一级参数中包含S₁中空间模型的坐标信息；

S₄、在S₃建立的禁忌搜索列表基础上，对S₁中空间模型各单元的取值及搜索范围进行约定，在建立的禁忌搜索原则中同时考虑观测数据曲线的形态，并对深浅地层结构的电阻率进行约定；

S₅、依据S₃建立的禁忌搜索列表和S₄建立的禁忌搜索原则，给S₁中构建的二维空间模型各单元赋电阻率值，建立电阻率反演初始模型ρ₀；

S₆、基于建立的模型进行大地电磁二维正演计算，计算地表观测点的电场E_x E_y、磁场H_x H_y，并根据电场、磁场值进一步计算视电阻率ρ和相位φ；

S₇、正演计算完成后，进行计算值与观测值的残差计算，并判断该残差是否达到设置值，如达到，则执行步骤S₁₀；如未达到，则执行步骤S₈；

S₈、采用互换定理，将网格节点与测点互换，利用正演一次求解，获取多节点的灵敏度，合成灵敏度矩阵；

S₉、依据S₈合成的灵敏度矩阵对S₅建立的电阻率反演初始模型中各单元的电阻率值进行修改，使各单元电阻率值ρ＝ρ₀+δρ，δρ是地下介质电阻率的修正量，模型参数修改后，返回步骤S₆，重新进行大地电磁二维正演计算；

S₁₀、输出反演结果：输出观测区域下方各单元的电阻率值ρ以及各单元的实际坐标值x,y。

其中，步骤S₃建立禁忌搜索列表的步骤为：

S₃₁、依据搜集的资料进行统计分析，初步确定工区岩性分类，并按工区岩性类型将整个工区划分为5-8个子区域；S₃₂、确定各子区域的地层层序和层厚度；S₃₃、各子区域以测井成果中的电阻率为基础，进行综合分析，确定各地层岩石电阻率，给出岩石电阻率参考范围；S₃₄、按照前述综合分析结果创建禁忌搜素列表，表内包括区域编号、岩性、平面坐标范围、地层层序、层厚度和地层岩石电阻率范围。

其中，步骤S₄建立禁忌搜索原则的步骤为：S₄₁、根据平面坐标值确定计算单元所处位置，并在禁忌搜索列表中查找对应信息；S₄₂、根据禁忌搜索列表给计算单元的电阻率属性赋值；S₄₃、根据大地电磁观测数据曲线形态，在纵向上给各计算单元的电阻率值乘以加权系数。优选的是，对于曲线起伏较大区段，所述加权系数取0.8-0.9；对于曲线起伏较小区段，所述加权系数取1.0-1.1。

其中，步骤S₆中大地电磁二维正演的方法为：

TE模式下，在S₅构建的电阻率反演初始模型基础上求解下述方程，未知数为E_y，求解地表或地下空间各节点的电场值：

TM模式下，在S₅构建的电阻率反演初始模型基础上求解下述方程，未知数为H_y，求解地表或地下空间各节点的磁场值：

其中，

σ为电导率，μ为磁场率，i为虚数单位，圆频率ω＝2πf，π为圆周率，f为频率；

上述方程求解中，采用矩形网格剖分，8个节点，双二次插值建立形函数，进而有限元求解：

形函数N_i构造如下：

式中，ε,η为规则变换网格中的坐标轴；

将以上形函数代入TE、TM方程，对方程组求解，得到电场E和磁场H，并根据电场、磁场值进一步计算视电阻率ρ和相位φ：

视电阻率

相位

步骤S8中，灵敏度矩阵的合成方法为：

TE模式下灵敏度矩阵方程如下：

TM模式下灵敏度矩阵方程如下：

式中

σ为电导率，μ为磁场率，i为虚数单位，圆频率ω＝2πf，π为圆周率，f为频率；

上述方程中

和

为走向电场和磁场对各单元电导率的偏导数；

采用互换定理求解该方程，将方程右端置于地表，与地下单元节点互换，此时，灵敏度矩阵方程与正演方程相似，通过m次正演计算得到地下节点的场

然后利用互换定理可得灵敏度矩阵

其中F′为

s_j为灵敏度矩阵方程中右端项，最后，将F′合成，即为灵敏度矩阵，其中m为地表测量点数量，m远小于n。

本发明的方法，在大地电磁反演迭代前，先利用钻探、测井、地质调绘、区域地质图等已知信息创建电阻率禁忌搜索列表，提出电阻率搜索原则，减小电阻率参数搜索范围，并分析研究区域地质信息，根据禁忌搜索列表构建反演初始模型，进行二维正演、残差计算，之后判断目标函数是否设定最小值，如达到，则输出地层电阻率、厚度值；如未达到，则根据灵敏度矩阵修改反演模型，不断迭代，直至目标函数达到设定的最小值。并且，在灵敏度矩阵计算中，采取互换源的方式，通过一次正演计算得到与一个测点对应的地下空间各节点的灵敏度。

本发明具有以下有益效果：

1、本发明利用钻探、测井、地质调绘、区域地质图等先知信息，分析确定了工区岩性分类、地下各地层层序、层厚度及各地层岩石电阻率值范围大概信息，构建了电阻率禁忌搜索列表，通过电阻率禁忌搜索列表以及禁忌搜索原则的措施缩小了反演迭代中电阻率参数的取值范围，避免了无目标、纯数学式搜索，提高了计算效率；

2、本发明将先知条件应用至大地电磁反演迭代过程中，相当于利用已知条件进行约束，提高了反演结果的准确度；

3、在反演灵敏度矩阵计算中，将源置于测量点位置，采用互换定理，与地下单元节点互换，通过一次正演计算得到地下单元对一个测点处的偏导数，从而大幅提高了反演迭代速度和反演效率。

附图说明

图1是本发明的基于禁忌搜索的大地电磁二维反演方法的流程图；

图2是本发明的大地电磁正演中形函数代号示意图；

图3a是本发明计算灵敏度矩阵中接收点与地下节点场源互换前的示意图；

图3b是本发明计算灵敏度矩阵中接收点与地下节点场源互换后的示意图。

具体实施方式

下面结合附图和具体实施例对本发明的基于禁忌搜索的大地电磁二维反演方法进行详细说明。

下述各实例仅用于说明本发明而并非对本发明的限制。

图1是本发明的基于禁忌搜索的大地电磁二维反演方法的流程图。如图1所示，本发明的基于禁忌搜索的大地电磁二维反演方法，包括以下步骤：

S₁、构建二维空间模型：

将地下二维空间从横、纵两个方向剖分为矩形网格，在观测目标区域内网格大小为等间隔，根据观测目标，设定地下空间单元网格大小，网格大小与拟探测目标体大小相当；在观测区域外部，随着距离增大网格大小也逐渐线性增大。

设定剖分模型左上角节点为坐标系原点(0,0)，横向向右为正、纵向向下为正，然后计算各单元节点的坐标值。

S₂、搜集测井、钻探、地质调查资料：搜集测量工区以往的钻探记录、测井成果以及地质调查资料。

S₃、创建禁忌搜索列表，包括：

依据前期搜集的资料进行综合分析，初步确定工区岩性分类及占比情况；

分析确定地下各地层层序、层厚度；

建立各岩层岩石电阻率参考范围：以测井电阻率成果为基础，综合考虑工区岩性和各区域地层情况，初步确定各地层岩石电阻率，给出岩石电阻率参考范围，然后根据这些信息建立二维空间模型中各单元的禁忌搜索列表，列表中的参数包括地层厚度范围、地层电阻率范围以及这些参数的禁忌范围，这些参数的下一级参数中包含S₁中空间模型的坐标信息。

S₄、建立禁忌搜索原则：在S₃创建的禁忌搜索列表基础上，对S₁中空间模型各单元的取值及搜索范围进行约定，由于观测数据曲线为随频率变化的曲线，高频代表浅地表、低频反映深地层的电阻率信息，所以，在建立禁忌搜索原则中同时考虑观测数据曲线的形态，对深浅地层结构的电阻率进行约定。具体如下：

S₄₁、判定单元所处岩性区域：根据平面坐标值确定计算单元所处位置，并在禁忌搜索列表中查找对应信息，进行岩性区域判定。

S₄₂、从禁忌搜索表取值：根据单元所处位置、区域岩性，在S₃₃中建立的禁忌搜索表内取值，给单元电阻率赋值。

S₄₃、电阻率值加权：大地电磁观测数据曲线为随频率变化的曲线，高频代表浅地表、低频反映深地层的电阻率信息，从地表点位给纵向各单元的电阻率值乘以加权系数，加权系数的取值主要依据观测数据曲线的形态，曲线起伏较大区段，加权系数取0.8-0.9，曲线起伏较小区段取1.0-1.1。

S₅、依据S₃建立的禁忌搜索列表和S₄建立的禁忌搜索原则，给S₁中构建的二维空间模型各单元赋电阻率值，建立电阻率反演初始模型ρ₀。

S₆、大地电磁二维正演：基于所述电阻率反演初始模型进行大地电磁二维正演计算，计算地表观测点的电场E_x E_y、磁场H_x H_y：

(1)TE模式下，在S₅构建模型的基础上，求解下述方程，未知数为E_y，求解地表或地下空间各节点的电场值。

(2)TM模式下，在S₅构建模型的基础上，求解下述方程，未知数为H_y，求解地表或地下空间各节点的磁场值。

其中，

σ为电导率，μ为磁场率，i为虚数单位，圆频率ω＝2πf，π为圆周率，f为频率。

(3)上述方程求解中，采用矩形网格剖分，8节点，如图2所示，图2为本发明的大地电磁正演中形函数编号示意图，其中，图2(a)的局部自然坐标系中的正方形单元称为计算单元；图2(b)中在整体坐标系中的任意四边形单元称为实际单元。两种单元内的参变量函数采用8个节点参数和相同的形函数进行变换。本发明采用双二次插值建立形函数，

进行等参变换。然后可以在局部自然坐标系中的计算单元上进行单元分析，然后再映射到实际直角坐标系中的单元上进行计算，这种等参变化具有计算精度高和适用性好的特点。

双二次插值形函数N_i构造如下：

式中，ε,η为规则变换网格中的坐标轴。

将以上形函数代入TE、TM方程，对方程组求解，即可得到电场E和磁场H，并根据电场、磁场值进一步计算视电阻率ρ和相位φ。

视电阻率

和相位

S₇、残差计算：正演计算完成后，进行计算值与观测值的残差计算，具体计算公式如下：

e＝ρ_obs-ρ_res

其中：e为残差，ρ_obs为视电阻率观测值，ρ_res为视电阻率正演计算值。

然后，判断该残差是否达到设置值1.0×10^-5，该值可根据数据质量灵活选择，数据质量较好时可向更小值选择，如达到，则执行步骤S₁₀；如未达到，则执行步骤S₈，进行灵敏度矩阵计算；

S₈、采用互换定理，将网格节点与测点互换，利用正演一次求解，获取多节点的灵敏度，减小计算量。

图3a和图3b是本发明计算灵敏度矩阵中接收点与地下节点场源互换前后的示意图。图3a中K单元的S₁节点变化对地表测点产生的响应F_s1为灵敏度，F_s1为走向电场和磁场对各单元电导率的偏导数

和

即灵敏度矩阵求解方程。

TE模式下灵敏度矩阵方程如下：

TM模式下灵敏度矩阵方程如下：

式中

σ为电导率，μ为磁场率，i为虚数单位，圆频率ω＝2πf，π为圆周率，f为频率。

上述方程中

和

为走向电场和磁场对各单元电导率的偏导数，该方程为灵敏度矩阵，其中的未知数数量为地下空间单元个数，每变换一次单元区域，就需要求解一次方程组，如果直接求解，计算量极大。但该方程与正演方程形式相似，因此，可以采用互换定理进行计算，即将源置于测量点位置，与地下单元节点互换，如图3b中S₁测点变化对地下节点产生的响应

通过互换定理可得

(F′为

s_j为灵敏度矩阵方程中右端项)，最后，将F′合成，即为灵敏度矩阵。这样，通过一次正演计算就可以得到地下单元对一个测点处的偏导数F′，如果有m个测点则只需做m次正演计算，就可以形成灵敏度矩阵P，如此计算将会大幅提升计算效率。

S₉、依据S₈合成的灵敏度矩阵对S₅建立的电阻率反演初始模型中各单元的电阻率值进行修改，使各单元电阻率值ρ＝ρ₀+δρ，δρ是地下介质电阻率的修正量，模型参数修改后，返回步骤S₆，重新进行大地电磁二维正演计算；

S₁₀、输出反演结果：输出观测区域下方各单元的电阻率值ρ以及各单元的实际坐标值x,y。

Claims (5)

1.一种基于禁忌搜索的大地电磁二维反演方法，包括以下步骤：

S₁、构建二维空间模型：将地下二维空间从横、纵两个方向剖分为矩形网格，在观测目标区域内网格大小为等间隔，且网格大小与拟探测目标体大小相当；在观测区域外部，随着距离增大网格大小也逐渐线性增大；

S₂、搜集测量工区以往的钻探记录、测井成果以及地质调查资料；

S₃、依据步骤S₂搜集的资料进行综合分析，初步确定工区岩性分类、地下各地层层序、层厚度及各地层岩石电阻率值范围信息，根据这些信息建立二维空间模型中各单元的禁忌搜索列表，该禁忌搜索列表中的参数包括地层厚度范围和地层电阻率范围，这些参数的下一级参数中包含S₁中空间模型的坐标信息；

S₄、在S₃建立的禁忌搜索列表基础上，对S₁中空间模型各单元的取值及搜索范围进行约定，在建立的禁忌搜索原则中同时考虑观测数据曲线的形态，并对深浅地层结构的电阻率进行约定；

S₅、依据S₃建立的禁忌搜索列表和S₄建立的禁忌搜索原则，给S₁中构建的二维空间模型各单元赋电阻率值，建立电阻率反演初始模型ρ₀；

S₆、基于建立的模型进行大地电磁二维正演计算，计算地表观测点的电场E_xE_y、磁场H_xH_y，并根据电场、磁场值进一步计算视电阻率ρ和相位φ；

S₇、正演计算完成后，进行计算值与观测值的残差计算，并判断该残差是否达到设置值，如达到，则执行步骤S₁₀；如未达到，则执行步骤S₈；

S₈、采用互换定理，将网格节点与测点互换，利用正演一次求解，获取多节点的灵敏度，合成灵敏度矩阵；

S₉、依据S₈合成的灵敏度矩阵对S₅建立的电阻率反演初始模型中各单元的电阻率值进行修改，使各单元电阻率值ρ＝ρ₀+δρ，δρ是地下介质电阻率的修正量，模型参数修改后，返回步骤S₆，重新进行大地电磁二维正演计算；

S₁₀、输出反演结果：输出观测区域下方各单元的电阻率值ρ以及各单元的实际坐标值x,y，

其中，步骤S8中灵敏度矩阵的合成方法为：

TE模式下灵敏度矩阵方程如下：

TM模式下灵敏度矩阵方程如下：

式中

σ为电导率，μ为磁场率，i为虚数单位，圆频率ω＝2πf，π为圆周率，f为频率；

上述方程中

和

为走向电场和磁场对各单元电导率的偏导数；

采用互换定理求解该方程，将方程右端置于地表，与地下单元节点互换，此时，灵敏度矩阵方程与正演方程相似，通过m次正演计算得到地下节点的场

然后利用互换定理可得灵敏度矩阵

其中F′为

s_j为灵敏度矩阵方程中右端项，最后，将F′合成，即为灵敏度矩阵，其中m为地表测量点数量，m远小于n。

2.根据权利要求1所述的大地电磁二维反演方法，其特征在于，步骤S₃建立禁忌搜索列表的步骤为：

S₃₁、依据搜集的资料进行统计分析，初步确定工区岩性分类，并按工区岩性类型将整个工区划分为5-8个子区域；

S₃₂、确定各子区域的地层层序和层厚度；

S₃₃、各子区域以测井成果中的电阻率为基础，进行综合分析，确定各地层岩石电阻率，给出岩石电阻率参考范围；

S₃₄、按照前述综合分析结果创建禁忌搜素列表，表内包括区域编号、岩性、平面坐标范围、地层层序、层厚度和地层岩石电阻率范围。

3.根据权利要求1所述的大地电磁二维反演方法，其特征在于，步骤S₄建立禁忌搜索原则的步骤为：

S₄₁、根据平面坐标值确定计算单元所处位置，并在禁忌搜索列表中查找对应信息；

S₄₂、根据禁忌搜索列表给计算单元的电阻率属性赋值；

S₄₃、根据大地电磁观测数据曲线形态，在纵向上给各计算单元的电阻率值乘以加权系数。

4.根据权利要求3所述的大地电磁二维反演方法，其特征在于：对于曲线起伏较大区段，所述加权系数取0.8-0.9；对于曲线起伏较小区段，所述加权系数取1.0-1.1。

5.根据权利要求1所述的大地电磁二维反演方法，其特征在于，步骤S₆中大地电磁二维正演的方法为：

TE模式下，在S₅构建的电阻率反演初始模型基础上求解下述方程，未知数为E_y，求解地表或地下空间各节点的电场值：

TM模式下，在S₅构建的电阻率反演初始模型基础上求解下述方程，未知数为H_y，求解地表或地下空间各节点的磁场值：

其中，

σ为电导率，μ为磁场率，i为虚数单位，圆频率ω＝2πf，π为圆周率，f为频率；

上述方程求解中，采用矩形网格剖分，8个节点，双二次插值建立形函数，进而有限元求解：

形函数N_i构造如下：

式中，ε,η为规则变换网格中的坐标轴；

将以上形函数代入TE、TM方程，对方程组求解，得到电场E和磁场H，并根据电场、磁场值进一步计算视电阻率ρ和相位φ：

视电阻率

相位

Images