CN112285788B - 一种基于电磁波动方程的cpml吸收边界条件加载方法 - Google Patents

Abstract

本发明涉及一种基于电磁波动方程的CPML吸收边界条件加载方法，采用电磁波动方程作为控制方程并基于有限差分方法进行数值模拟，将整个计算区域分为中心区域和边界区域；在中心区域，求解三维电磁波动方程，得到中心区磁场垂直分量波场；在边界区域，将复拉伸变量代入频率域电磁波动方程，采用CPML吸收边界条件设置复拉伸变量，并将其表达式代入控制方程中，进行整理并频时变换，最后基于有限差分算法进行离散近似，得到边界区磁场垂直分量波场；将中心区和边界区波场叠加获得最终的波场。本发明目的在于克服电磁数值模拟时计算效率低及晚期反射误差大等问题，实现三维时域电磁响应的高效、高精度数值模拟。

Description

一种基于电磁波动方程的CPML吸收边界条件加载方法

技术领域

本发明涉及一种基于电磁波动方程的CPML吸收边界条件加载方法，适用于时域电磁响应数值模拟，尤其是对磁性源电磁响应的数值模拟。

背景技术

磁性源瞬变电磁探测方法，利用不接地回线向地下发射双极性脉冲电流，在发射脉冲产生的一次磁场间歇过程中，通过接收线圈与接收机采集二次感应场随时间变化的衰减曲线，对数据进行处理解释获得地下介质丰富的电性信息，由于其具有探测深度大、经济、便捷的优点，已经被广泛应用于地质资源勘探中。卷积完善匹配层(CPML，ConvolutionPerfect Matching Layer)是一种稳定高效的吸收边界条件，它是一种有限厚度的特殊媒质，包围计算空间，是基于一种虚构的本构参量来创建波匹配条件，这种匹配条件与波的频率和在边界上的入射角度无关，解决了空间截断引起数值反射并导致计算空间晚期结果不准确的问题，随着对地下介质异常特征的精细化研究，计算空间无吸收边界条件或吸收边界条件不完善已经不能满足高精度地质探测的要求。建立更符合电磁波传播规律的数学模型是提高电磁探测分辨率的关键技术之一。

电磁波动方程是由麦克斯韦基本方程组导出，表示脱离场源的电磁波在无源空间的传播规律和特点，电磁波动方程的特征是把电场与磁场分离出来，同时把电场与磁场的三个方向分量分离出来，可以只计算磁场垂直分量的响应特性，计算内存占用是麦克斯韦方程的六分之一左右，显著提高了计算效率，在实际应用中发现，电磁波动方程对地下异常体的响应分辨率明显高于麦克斯韦方程，适合对大计算空间进行高效高精度电磁探测。

三维时域有限差分(Three-dimensional Finite-difference Time-domain)方法是电磁场数值计算的重要方法之一，它的基本思想是场量对时间和空间的一阶和二阶偏导数用中心差分近似,通过在时域的迭代模拟波的传播过程,从而得出场量分布信息，与一维和二维差分方法相比，三维差分更适用于计算空间大、结构复杂的地质探测分析中，目前已被广泛应用于电磁场数值模拟计算中。

王洪华等使用二维电磁波动方程和有限元方法对探地雷达进行了数值模拟，并对比了无边界条件、Sarma边界条件和PML边界条件下均匀模型的时域反射误差和能量衰减曲线，结果表明PML边界条件的吸收效果良好，具有近似零反射系数。曹敏敏等采用隐式中心差分方法研究了具有周期性边界条件的二阶电磁波动方程的守恒性和稳定性，通过数值实验从数学的角度证明了电磁波动方程的实用性。国内外很少用三维电磁波动方程对时域电磁响应进行数值模拟，也很少给电磁波动方程加入吸收边界条件。

中国专利CN109116418B公开了一种三维波动方程混合网格有限差分数值模拟方法及装置，通过构建三维混合网格有限差分格式获得三维波动方程的有限差分离散方程；根据平面波理论计算离散方程的差分系数，对离散方程进行求解并实现对三维波动方程的数值模拟，采用混合网格的方法降低了三维波动方程数值模拟的数值频散，提高了计算精度。

中国专利CN103616721A公开了一种基于二阶偏微分波动方程的PML吸收边界条件的方法，通过加载某一采样时刻震源的纵波波场，计算纵波波场在三维空间某一方向的二阶偏导数；根据二阶偏导数对此方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减，继续对不同方向、不同时刻的纵波波场进行吸收衰减，计算结果证明了吸收边界条件的实用性和有效性。

中国专利CN109188517A公开了一种基于Higdon余弦型加权的混合吸收边界条件方法，对地震勘探进行正演模拟，通过采用二维标量声波方程进行有限差分数值模拟，将研究区域划分为中心波场区域、过渡区域及边界区域，对不同区域分别求解二维标量声波双程波方程及二阶Higdon吸收边界方程；在过渡区域，将双程波波场和单程波波场进行余弦型加权叠加，得到最终的波场，这种方法使得内边界与中心波场区域耦合更好，有效压制了内边界反射。

以上所述方法公布了国内外关于吸收边界条件加载的数值模拟方法。但对三维电磁波动方程中的边界条件加载，目前几乎没有关于在时间域进行三维有限差分运算的研究，且对于电磁数值模拟在运算时具有占用内存大导致计算区域有限及分辨率小等问题，如何在电磁场正演模拟领域中高效、准确的进行高分辨率有限差分运算，是本领域技术人员迫切解决的一个技术问题。

发明内容

本发明所要解决的技术问题在于提供一种基于电磁波动方程的CPML吸收边界条件加载方法，采用电磁波动方程作为控制方程并基于有限差分方法进行数值模拟，将整个计算区域分为中心区域和边界区域；在中心区域，求解三维电磁波动方程，得到中心区磁场垂直分量波场；在边界区域，将复拉伸变量代入频率域电磁波动方程，采用CPML吸收边界条件设置复拉伸变量，并将其表达式代入控制方程中，进行整理并频时变换，最后基于有限差分算法进行离散近似，得到边界区磁场垂直分量波场；将中心区和边界区波场叠加获得最终的波场。本发明目的在于克服电磁数值模拟时计算效率低及晚期反射误差大等问题，实现三维时域电磁响应的高效、高精度数值模拟。

本发明是这样实现的，一种基于电磁波动方程的CPML吸收边界条件加载方法：

1)、对Maxwell方程组进行变换得到电磁波动方程作为磁场垂直分量迭代控制方程，电磁波动方程将电场与磁场分离开来，且电场与磁场的各个分量可以独立运算，在计算过程中可以大幅降低计算量就可以获得感应电动势的衰减曲线；

2)、对磁场垂直分量控制方程中的时间和空间偏导项用中心差分形式进行离散，并推导出磁场垂直分量迭代方程；

3)、对电磁波动方程进行时-频变换，将复拉伸变量代入电磁波动方程，依据CPML吸收边界条件选取复拉伸变量，将其表达式代入方程中，对控制方程进行整理并进行频-时变换,完成电磁波动方程吸收边界条件的加载；

4)、采用有限差分方法对控制方程时间和空间偏导项进行离散，递归得出含CPML吸收边界条件时磁场垂直分量的迭代方程；

5)、采用非均匀三维Yee氏网格对计算区域进行剖分，设置计算域电导率、磁导率、人工介电常数和CPML参数，计算初始场；

6)、使用GPU加速，并在观测时间内开展磁场垂直分量迭代运算，在边界层应用含CPML吸收边界条件时磁场垂直分量的控制方程；

7)、迭代计算结束后，提取磁场垂直分量响应并进行成图，对结果进行分析处理；

其中步骤1)中，电磁波动方程磁场垂直分量表达式为：

式(1)中▽²为拉普拉斯算子，H_z为磁场垂直分量，μ为磁导率，ε为介电常数，σ为电导率，t为时间；电磁波动方程可以表示磁场的垂直分量控制方程，大幅降低计算量，提高计算效率并增大计算空间；

进一步地，步骤2)中将控制方程中的空间和时间偏导项进行离散可得：

其中Δt_n+1＝t_n+1-t_n，Δt_n＝t_n-t_n-1，Δx₂,Δx₁,Δy₂,Δy₁,Δz₂,Δz₁为各方向非均匀网格步长；

将控制方程进行展开，进行差分离散近似，得到磁场垂直分量迭代方程：

其中步骤3)中对电磁波动方程进行时-频变换及代入复拉伸变量过程为：

其中

为H_z关于时间的傅里叶变换，ω为角频率，S_x、S_y、S_z分别为计算空间各个方向的复拉伸变量，其表达式为：

其中σ_pei为CPML层电导率，ε₀为真空中介电常数，κ_ei、α_ei为CPML吸收边界条件引入的参数，κ_ei的取值大于等于1，α_ei的取值大于等于0；

式(8)两边同乘以S_xS_yS_z，并整理可以得到：

将式(9)代入式(10)，整理后可得：

其中A_i,B_i,C_i,D_i,(i＝0,1,2,3,4,5,6,7)为公式整理后的复拉伸变量中的参数组合，其中一组参数的表达式为；

对式(11)进行通分并省略三阶及三阶以上高阶项得到：

其中P_i,Q_i,R_i,Q₄,R₄,(i＝0,1,2,3)为公式整理后的系数，其中一组系数表达式为：

对式(13)进行频-时变换得到含CPML吸收边界条件时磁场垂直分量的控制方程：

其中步骤4)中采用有限差分方法对式(15)中偏导项进行离散，将磁场控制方程进行展开并整理可得：

其中K,L_i,M_i,N_i,(i＝0,1,2,3)为公式整理后的系数，其中一组系数表达式为：

移项后得到磁场垂直分量迭代方程：

本发明与现有技术相比，有益效果在于：在时域电磁三维数值模拟中使用电磁波动方程，能够明显提高计算效率并增大计算空间；针对电磁数值模拟过程中空间截断引起数值反射并导致计算空间晚期结果不准确的问题，引入CPML吸收边界条件，可以最大限度的消除截断边界对数值模拟结果的影响，提高数值模拟的精确度，同时电磁波动方程与CPML吸收边界条件的组合显著提高了对异常体的电磁响应分辨率，达到了对三维时域电磁响应的高效、高精度、大分辨率数值模拟的目的。

附图说明

图1是一种基于电磁波动方程的CPML吸收边界条件加载方法示意图；

图2是基于CPML吸收边界条件的三维时域电磁波动方程在均匀半空间中(电导率为0.01S/m)数值模拟得到的接收线圈感应电动势衰减曲线与数值积分解对比图及误差曲线；

图3是不含吸收边界条件的三维时域电磁波动方程在均匀半空间中数值模拟得到的地面二维电磁响应平面图；

图4是含CPML吸收边界条件的三维时域电磁波动方程在均匀半空间中数值模拟得到的地面二维电磁响应平面图；

图5是基于CPML吸收边界条件的三维时域电磁波动方程对不同电导率的异常体(异常体大小为300m×300m×300m，中心分布，上表面距地面高度230m)数值模拟得到的接收线圈感应电动势衰减曲线对比图；

图6是基于CPML吸收边界条件的三维时域电磁波动方程对异常体数值模拟得到的电磁响应切面图；

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

实施例

参见图1，一种基于电磁波动方程的CPML吸收边界条件加载方法，包括：

1)、根据电磁波动方程推导磁场垂直分量迭代方程，结合非均匀Yee式网格和DuFort-Frankel方法构建磁场垂直分量显式差分格式，建立7390m×7390m×4175m的计算区域，网格数目为161×161×85，其中x、y方向上的网格数均为161个，z方向上网格数85个，最小网格步长为10m，最大网格步长为120m，计算区域中除地面外其余5个面的最外层6个网格加载CPML边界条件。

2)、在整个计算区域内设置电导率、磁导率、人工介电常数、CPML系数等参数，大地电导率设置为σ₀＝0.01S/m，磁导率设置为真空磁导率。

3)、将计算区域的电性参数代入到发射和接收部分均在地面的电场响应表达式中，计算t₀、t₁时刻的x，y方向的电场响应作为初始场。

E_z＝0 (21)

其中I为发射电流，μ和σ为均匀半空间的磁导率和电导率，W和L分别为x和y方向的积分范围，J₀为0阶贝塞尔函数，iω为拉普拉斯变量，λ为汉克尔变换积分变量。

4)、将参与迭代的CPU序列转化为GPU序列，加速运算。

5)、将电场初始场值代入麦克斯韦方程中磁场控制方程，在整个计算区域更新磁场值H_z，并保存t₀、t₁时刻磁场值。

6)、计算当前时刻迭代时间步长。

Δ_min为最小网格步长，a取值为0.1。

7)、将t₀、t₁时刻磁场值代入磁场垂直分量迭代方程中，在计算区域中心层更新磁场值H_z，并保存当前时刻及前一时刻磁场值。

8)、将t₀、t₁时刻磁场值代入含CPML吸收边界条件时磁场垂直分量迭代方程中，在计算区域边界层更新磁场值H_z，并保存当前时刻及前一时刻磁场值。

9)、判断是否完成全部迭代，若没有完成，重复6-8步骤，如果完成全部迭代，将计算结果从GPU释放出来，并对计算结果进行显示、成图，完成全部计算。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种基于电磁波动方程的CPML吸收边界条件加载方法，其特征在于，包括如下步骤：

1)、对Maxwell方程组进行变换得到电磁波动方程作为磁场垂直分量迭代控制方程，电磁波动方程将电场与磁场分离开来，且电场与磁场的各个分量可以独立运算，在计算过程中可以大幅降低计算量就可以获得感应电动势的衰减曲线；

2)、对磁场垂直分量控制方程中的时间和空间偏导项用中心差分形式进行离散，并推导出磁场垂直分量迭代方程；

3)、对电磁波动方程进行时-频变换，将复拉伸变量代入电磁波动方程，依据CPML吸收边界条件选取复拉伸变量，将其表达式代入方程中，对控制方程进行整理并进行频-时变换,完成电磁波动方程吸收边界条件的加载；

4)、采用有限差分方法对控制方程时间和空间偏导项进行离散，递归得出含CPML吸收边界条件时磁场垂直分量的控制方程；

5)、采用非均匀三维Yee氏网格对计算区域进行剖分，设置计算域电导率、磁导率、人工介电常数和CPML参数，计算初始场；

6)、使用GPU加速，并在观测时间内开展磁场垂直分量迭代运算，在边界层应用含CPML吸收边界条件时磁场垂直分量的迭代方程；

7)、迭代计算结束后，提取磁场垂直分量响应并进行成图，对结果进行分析处理；

其中步骤1)中，电磁波动方程磁场垂直分量表达式为：

式(1)中

为拉普拉斯算子，H_z为磁场垂直分量，μ为磁导率，ε为介电常数，σ为电导率，t为时间；电磁波动方程可以表示磁场的垂直分量控制方程，大幅降低计算量，提高计算效率并增大计算空间；

进一步地，步骤2)中将控制方程中的空间和时间偏导项进行差分离散近似，得到磁场垂直分量迭代方程:

其中步骤3)中对电磁波动方程进行时-频变换及代入复拉伸变量过程为：

其中

为H_z关于时间的傅里叶变换，ω为角频率，S_x、S_y、S_z分别为计算空间各个方向的复拉伸变量，其表达式为：

其中i为x、y或z，表示空间不同的方向，σ_pei为CPML层电导率，ε₀为真空中介电常数，κ_ei、α_ei为CPML吸收边界条件引入的参数，κ_ei的取值大于等于1，α_ei的取值大于等于0；

式(3)两边同乘以S_xS_yS_z，并整理可以得到：

将式(4)代入式(5)，整理后可得：

其中A_i，B_i，C_i，D_i为公式整理后的复拉伸变量中的参数组合，其中i为0、1、2、3、4、5、6、7，表示不同的参数组合，其中一组参数的表达式为：

对式(6)进行通分并省略三阶及三阶以上高阶项得到：

其中P_i,Q_i,R_i,Q₄,R₄为公式整理后的系数，其中i为0、1、2、3，表示不同的系数，其中一组系数表达式为：

对式(8)进行频-时变换得到含CPML吸收边界条件时磁场垂直分量的控制方程：

其中步骤4)中采用有限差分方法对式(10)中偏导项进行离散差分近似，得到含CPML吸收边界条件磁场垂直分量迭代方程：

其中K,L_i,M_i,N_i为公式整理后的系数，其中i为0、1、2、3，表示不同的系数，其中一组系数表达式为：

Images