CN101576622A - 一种超宽带电磁波的模拟方法 - Google Patents

Abstract

本发明公开了一种超宽带电磁波的模拟方法，该方法是在伸展坐标系下的PML方法和高阶时域有限差分方法FDTD(2，2M)的差分格式基础上提出了一种新的PML方法：自适应吸收参数下的分段线性卷积方式卷积完全匹配层(PCPML，Piecewise linearity CPML)吸收边界条件，该方法设定PML吸收因子是可调参数以吸收从高频到低频之间的不同类型的超宽带电磁波，并由实码遗传算法优化方法与线性迭代方法的混合方法自适应方式得到该参数，以适应超宽频带的电磁波的吸收。本发明所采用的吸收边界条件对高、低频波均能吸收处理，并可减少匹配层由网格色散和边界反射造成的反射误差，提高吸收的效能和模拟结果的准确度。

Description

一种超宽带电磁波的模拟方法

技术领域

本发明属于超宽带电磁波三维数值模拟技术范畴，是一种以模拟超宽带电磁波在同时含有空气和地下复杂有耗介质中传播的三维电磁勘探正演方法。

背景技术

用超宽带电磁波探测地下目标体或穿墙探测生命体是一项尖端技术，可用于地球物理、超宽带天线设计、雷达目标探测、穿墙生命体搜索等方面。目前，模拟超宽带电磁波在复杂介质中传播时，普遍存在着截断边界处理困难、长时反射信号精度差等问题。

超宽带电磁脉冲是一种瞬态电磁波，超宽带电磁法的模拟方法最好是直接在时间域内进行。而时域有限差分(FDTD)方法是分析瞬态电磁波的一种重要方法，也是超宽带电磁法的数值模拟方法的首选方法。目前对超宽带电磁波的模拟研究不多，集中在近天线区场的FDTD模拟上，为超宽带天线设计进行指导，但这些研究的模拟范围不大，空间介质模型较简单、反射延时不长，吸收边界易于处理。但对探地而言，需要计算的空间范围大、介质模型地电分布极其复杂、反射延时时间很长、晚时回波信号精度要求高等，这时，一般的吸收边界处理方法如Mur(G.Mur，1981)吸收边界方法、超吸收边界方法、一般的PML(Perfectly Matched Layer，1994)方法等无法处理复杂介质的超宽带电磁波的吸收问题。复杂介质中传播的超宽带电磁波的频带很宽，倏逝波、高低频的波都很丰富，特别是低频电磁波，一般的吸收方法难以统一吸收处理，长时反射波受网格色散和边界反射的影响，特别是边界反射的影响极大，导致模拟精度很低；另外，传统的PML吸收边界处理方法需要人为设置吸收介质的层数、方向电导率等吸收参数，对任意复杂有耗模型，很难人为设定能够合理地模拟超宽带电磁波的吸收参数，人们希望计算机能根据给定的地电模型自动给出优化后的吸收参数，以得到准确的模拟结果。因此研究合适的吸收边界方法对超宽带电磁法的数值模拟至关重要，必须要采用一种更为有效吸收边界处理方法和高精度的FDTD计算方法。

发明内容

本发明的目的在于提供一种超宽带电磁波的模拟方法，该方法所采用的吸收边界条件对高、低频波均能吸收处理，并可减少匹配层由网格色散和边界反射造成的反射误差，提高吸收的效能和模拟结果的准确度。

本发明利用空间2M阶的高阶FDTD方法模拟电偶极子天线激励的超宽带电磁波的传播，其吸收边界处理采用分段线性卷积方式卷积完全匹配层吸收边界方法即PCPML方法，该方法包括以下步骤：

(1)输入初始参数：给定三维模型的网格化的地电参数、背景模型的地电参数、发射电偶极天线参数、激励脉冲信号参数、时间步长、最大时间步数或模拟电磁波的最大传播时间、PML层数、PCPML吸收参数搜索范围及遗传优化算法的计算参数；

(2)计算优化PML参数的目标函数参数：计算电偶极子发射天线激发的超宽带电磁脉冲的电磁响应，作为遗传优化的目标函数中的参考响应值；

(3)使用实码混合遗传算法并结合线性反演算法、以及采用精英选择策略并结合线性反演的局部优化种群过程，优化计算PML参数，该PML参数包括各PML层内的方向电导率(σ_i)，对表面波的吸收的可调因子(κ)和对FDTD难以吸收的低频波的可调因子(α_i)；当目标函数值小于给定的误差限时，遗传算法结束；此时令时间步n＝0，开始FDTD迭代；

(4)用FDTD方法迭代计算第n个时间步的三维空间电场值：根据优化得到PML参数，用空间2M阶的高阶差分方法和分段线性递归卷积吸收边界方法得到的FDTD公式依次计算第n个时间步的三个电场分量Ex、Ey和Ez；

(5)用FDTD方法迭代计算第n+1/2个时间步的三维空间磁场值：按FDTD(2，2M)方法迭代计算第n+1/2个时间步的三个磁场分量Hx、Hy和Hz；

(6)取n＝n+1，重复第(4)、(5)步骤，直到最大时间步Nt，完成FDTD迭代计算；

(7)保存电磁场数据，按常规成图方法得到相应的波场快照。

本发明是建立在如下理论基础上的：

1)自适应吸收参数下的分段线性卷积方式卷积完全匹配层(PCPML-Piecewise linearity CPML)吸收边界处理方法是对CPML(卷积PML)吸收边界条件的改进，以分段线性递归卷积方式计算伸展坐标PML中的卷积项而得到。

PCPML是建立在伸展坐标PML(Chew W.C.，Jin J.M.and MichielssenE.，1997)基础上的。

在伸展坐标PML中，直角坐标系下Maxwell的2个频域里的旋度方程的分量形式为：

$\mathrm{iw\ϵ}{\tilde{E}}_x+\mathrm{\σ}{\tilde{E}}_x=\frac{1}{s_y}\frac{\∂{\tilde{H}}_z}{\∂y}-\frac{1}{s_z}\frac{\∂{\tilde{H}}_y}{\∂z}---\left(1\right)$

$\mathrm{iw\ϵ}{\tilde{E}}_y+\mathrm{\σ}{\tilde{E}}_y=\frac{1}{s_z}\frac{\∂{\tilde{H}}_x}{\∂z}-\frac{1}{s_x}\frac{\∂{\tilde{H}}_z}{\∂x}$

$\mathrm{iw\ϵ}{\tilde{E}}_z+\mathrm{\σ}{\tilde{E}}_z=\frac{1}{s_x}\frac{\∂{\tilde{H}}_y}{\∂x}-\frac{1}{s_y}\frac{\∂{\tilde{H}}_x}{\∂y}$

$-\mathrm{iw}{\mathrm{\μ}}_0{\tilde{H}}_x=\frac{1}{s_y}\frac{\∂{\tilde{E}}_z}{\∂y}-\frac{1}{s_z}\frac{\∂{\tilde{E}}_y}{\∂z}$

$-\mathrm{iw}{\mathrm{\μ}}_0{\tilde{H}}_y=\frac{1}{s_x}\frac{\∂{\tilde{E}}_z}{\∂x}-\frac{1}{s_z}\frac{\∂{\tilde{E}}_x}{\∂z}$

$-\mathrm{iw}{\mathrm{\μ}}_0{\tilde{H}}_z=\frac{1}{s_x}\frac{\∂{\tilde{E}}_y}{\∂x}-\frac{1}{s_y}\frac{\∂{\tilde{E}}_x}{\∂y}$

其中：ε和σ分别是介质的介电常数和电导率，w为角频率，并假定计算空间磁导率为自由空间磁导率μ₀，CPML方法的伸张坐标因子s_i为：

$s_i=k_i+\frac{{\mathrm{\σ}}_i}{{\mathrm{\α}}_i+{\mathrm{iw\ϵ}}_0},(i=x,y,z)---(*)$

其中：σ_i为PML层内i方向电导率参数，是正实数，而在非PML区域，其值为0。并假定α_i取正实数，k_i≥1。k_i值的引入是为了改善PML对表面波的吸收特性，而α_i值则是为了改善PML对低频分量的吸收特性，而UWB脉冲探地雷达富含这两种波，特别是低频电磁波，因而使用CPML的效果远好于传统PML的效果。

下面以(1)为例推导CPML公式，把方程(1)变换到时域中：

${\mathrm{\ϵ}}_r{\mathrm{\ϵ}}_0\frac{{\∂E}_x}{\∂t}+{\mathrm{\σE}}_x={\hat{s}}_y\left(t\right)*\frac{{\∂H}_z}{\∂y}-{\hat{s}}_z\left(t\right)*\frac{{\∂H}_y}{\∂z}---\left(2\right)$

其中，“*”表示卷积，

是的逆拉普拉斯变换，可写为：

${\hat{s}}_i=\frac{\mathrm{\δ}\left(t\right)}{k_i}-\frac{{\mathrm{\σ}}_i}{{\mathrm{\ϵ}}_0{k}_i^2}\exp \left[(-\frac{{\mathrm{\σ}}_i}{{\mathrm{\ϵ}}_0{k}_i}+\frac{{\mathrm{\α}}_i}{{\mathrm{\ϵ}}_0})t\right]U\left(t\right)=\frac{\mathrm{\δ}\left(t\right)}{k_i}+{\mathrm{\ξ}}_i\left(t\right),(i=x,y,z)---\left(3\right)$

其中，δ(t)和U(t)分别是单位冲击函数和单位阶跃函数。

把(3)代入(2)式后得到：

${\mathrm{\ϵ}}_r{\mathrm{\ϵ}}_0\frac{{\∂E}_x}{\∂t}+{\mathrm{\σE}}_x=\frac{1}{k_y}\frac{{\∂H}_z}{\∂y}-\frac{1}{k_z}\frac{{\∂H}_y}{\∂z}+{\mathrm{\ξ}}_y\left(t\right)*\frac{{\∂H}_z}{\∂y}-{\mathrm{\ξ}}_z\left(t\right)*\frac{{\∂H}_y}{\∂z}---\left(4\right)$

对(4)式按某种差分算法进行差分离散运算后，便可得到相应分量的FDTD计算公式，对其它的方程，计算公式推导过程类似。

本发明PCPML方法是对(4)式里的卷积计算采用了近似线性卷积方式，在卷积项的各个时间段Δt内，CPML是假定磁场的偏导数项的值是恒定的，从而得到迭代公式。本发明则假定卷积项内的电磁场在时间段Δt内认为是线性变化的，同时采用空间2M阶高阶差分方法计算空间导数项，方法如下：将函数f(x)进行泰勒级数展开，不难得到适用于交错网格的f(x)的2M阶近似精度的一阶导数公式：

$f^{\′}\left(x\right)=\frac{1}{\mathrm{\ΔL}}\underset{i=-M}{\overset{M-1}{\mathrm{\Σ}}}c\left(i\right)f(x+i+0.5)---\left(5\right)$

其中， $c\left(i\right)=\frac{{(-1)}^i{[(2M-1)!!]}^2}{2{(i+0.5)}^2(2M-2-2i)!!(2M+2i)!!}---\left(6\right)$

现在用(5)式的差分方法，对(4)式进行离散化，但由于直接计算其中的卷积项过于费时，本发明沿用FDTD计算色散介质的线性卷积处理方法，以卷积项为例，它的离散化近似计算公式为：

${\mathrm{\ξ}}_y\left(t\right)*\frac{{\∂H}_z}{\∂y}={\∫}_0^{\mathrm{n\Δt}}{\mathrm{\ξ}}_y\left(\mathrm{\τ}\right)\frac{1}{\mathrm{\Δy}}\underset{i=-M}{\overset{M-1}{\mathrm{\Σ}}}c\left(i\right)H_z(y+i+0.5|\mathrm{n\Δt}-\mathrm{\τ})\mathrm{d\τ}---\left(7\right)$

$=\underset{m=0}{\overset{n-1}{\mathrm{\Σ}}}\frac{1}{\mathrm{\Δy}}\underset{i=-M}{\overset{M-1}{\mathrm{\Σ}}}c\left(i\right)H_z(y+i+0.5|(n-m)\mathrm{\Δt}){\∫}_{\mathrm{m\Δt}}^{(m+1)\mathrm{\Δt}}{\mathrm{\ξ}}_y\left(\mathrm{\τ}\right)\mathrm{d\τ}$

其中，H_z(y+i+0.5|nΔt-τ)和H_z(y+i+0.5|(n-m)Δt)分别表示卷积和项在nΔt-τ时间变量及其离散时刻点(n-m)Δt时，Hz分量对应y方向网格节点y+i+0.5处的场值。根据ξ_y(t)的定义，有：

$Z_y\left(m\right)={\∫}_{\mathrm{n\Δt}}^{(n+1)\mathrm{\Δt}}{\mathrm{\ξ}}_y\left(\mathrm{\τ}\right)\mathrm{d\τ}=a_y\exp [-(\frac{{\mathrm{\σ}}_y}{k_y}+{\mathrm{\α}}_y)\frac{\mathrm{m\Δt}}{{\mathrm{\ϵ}}_0}]---\left(8\right)$

其中， $a_i=\frac{{\mathrm{\σ}}_i}{{\mathrm{\σ}}_i{k}_i+k_i^2{\mathrm{\α}}_i}\left\{\exp \right[-(\frac{{\mathrm{\σ}}_i}{k_i}+{\mathrm{\α}}_i)\frac{\mathrm{\Δt}}{{\mathrm{\ϵ}}_0}]-1.0\},(i=x,y,z)$

类似地写出Z_z(m)的表达式。

这样，由PCPML方法和高阶差分方法的计算方式则可获得更高的计算精度，进而可减少匹配层由网格离散造成的反射误差，提高吸收的效能和模拟结果的准确度，降低差分计算的低通滤波的影响，使得计算的频带更宽，符合超宽带数值模拟的要求。

在时间段[iΔt，(i+1)Δt]内，连续的电磁场偏导数项A(t)(如式(7)里的

)可由节点场值Aⁱ和Aⁱ⁺¹近似表出：

$A\left(t\right)=A^i+\frac{A^{i+1}-A^i}{\mathrm{\Δt}}(t-\mathrm{i\Δt})---\left(9\right)$

显然比用在时间段Δt内场值函数取为恒定的CPML方法的计算结果要精确，如图1、图2所示，图1为场函数分段线性近似(实线)与CPML计算方法(虚线)的对比图，图2显示了卷积项的分段线性近似计算方法，可写为：

$A(\mathrm{n\Δt}-\mathrm{\τ})=A^{n-m}+\frac{A^{n-m-1}-A^{n-m}}{\mathrm{\Δt}}(\mathrm{\τ}-\mathrm{m\Δt})---\left(10\right)$

则卷积(7)可表示为：

其中， $X\left(m\right)={\∫}_{\mathrm{m\Δt}}^{(m+1)\mathrm{\Δt}}{\mathrm{\ξ}}_y\left(\mathrm{\τ}\right)\mathrm{d\τ}=a_y{e}^{-\mathrm{sm\Δt}},$ $s=\frac{{\mathrm{\σ}}_y}{k_y{\mathrm{\ϵ}}_0}+\frac{{\mathrm{\α}}_y}{{\mathrm{\ϵ}}_0}.$

$Y\left(m\right)=\frac{1}{\mathrm{\Δt}}{\∫}_{\mathrm{m\Δt}}^{(m+1)\mathrm{\Δt}}(\mathrm{\τ}-\mathrm{m\Δt}){\mathrm{\ξ}}_y\left(\mathrm{\τ}\right)\mathrm{d\τ}=\frac{1}{\mathrm{\Δt}}{\∫}_{\mathrm{m\Δt}}^{(m+1)\mathrm{\Δt}}\mathrm{\τ}{\mathrm{\ξ}}_y\left(\mathrm{\τ}\right)\mathrm{d\τ}$

$-m{\∫}_{\mathrm{m\Δt}}^{(m+1)\mathrm{\Δt}}{\mathrm{\ξ}}_y\left(\mathrm{\τ}\right)\mathrm{d\τ}=\frac{1}{\mathrm{\Δt}}{\∫}_{\mathrm{m\Δt}}^{(m+1)\mathrm{\Δt}}\mathrm{\τ}{\mathrm{\ξ}}_y\left(\mathrm{\τ}\right)\mathrm{d\τ}-\mathrm{mX}\left(m\right).$

$=-p_y{e}^{-\mathrm{sm\Δt}}-{\mathrm{ma}}_y{e}^{-\mathrm{sm\Δt}}$

其中 $p_y=\frac{{\mathrm{\σ}}_y{\mathrm{\ϵ}}_0(1-e^{-\mathrm{s\Δt}})}{({\mathrm{\σ}}_y{k}_y+{\mathrm{\α}}_y{k}_y^2)\mathrm{\Δt}}.$

改写(11)得：

$A^n{a}_y-[A^{n-1}-A^n]p_y$

$=\underset{m=0}{\overset{n-2}{\mathrm{\Σ}}}\left\{\right[(m+1)A^{n-1-m}-{\mathrm{mA}}^{n-1-m-1}+A^{n-1-m}-A^{n-1-m-1}]a_y{e}^{-\mathrm{sm\Δt}}---\left(12\right)$

$-[A^{n-1-m-1}-A^{n-1-m}]p_y{e}^{-\mathrm{sm\Δt}}\}{e}^{-\mathrm{s\Δt}}+A^n{a}_y-[A^{n-1}-A^n]p_y$

并记 $u^n=\underset{m=0}{\overset{n-2}{\mathrm{\Σ}}}[A^{n-m-1}-A^{n-m-2}]e^{-\mathrm{sm\Δt}},$ 它也通过递归方法计算得到：

$u^n=\left\{\begin{array}{cc}{A}^{n-1}-A^{n-2}+u^{n-1}{e}^{-\mathrm{s\&Delta;t}}& n\&GreaterEqual;2\\ 0& n<2\end{array}\right.---\left(13\right)$

这样，(7)式里的分段线性卷积的迭代计算公式为：

迭代时，的初值为0。(4)式里的另一个卷积项

可类似地求出，代入(4)后，即可得到Ex分量的迭代公式，分别对(1)式表示的Maxwell旋度方程其它各个式子进行类似的推导。

2)实数遗传算法与线性反演方法的原理

一般的PML方法对各PML参数(PML层数、各PML层内的方向电导率σ_i、对表面波的吸收的可调因子κ和对FDTD一般难以吸收的低频波的可调因子α_i)是人为给定的，具有盲目性，没有统一的选择规律可以遵循。本发明对PML参数采用混合线性反演的实数遗传算法进行优化得到，可得到的更好的参数选择结果。

先选择与模型背景较为接近的均匀或层状模型，对电偶极子源或磁偶极子源，解析计算背景模型电磁响应，并设为目标函数中的参考值。这样，便可计算各边界节点上的反射误差值，而目标函数取为各边界各节点上的反射误差之和的倒数：

个体与群体生成：给定待优化的各PML参数的变化范围和离散个数，在对应搜索范围内均匀分割，在各小区间端点序列中随机生成一个初始个体。假定其中一个参数的给定参数范围为[X_imin，X_imax]，进行均匀分割获得个体。在[0，Ni]中产生一个随机整数，Ni为第i个参数所在区间被分割的个数，即为个体的元素。把待优化参数对应的依次连在一起构成问题解的编码形式为(k₁，k₂，…，k_n)。在计算适应值时将k_i做映射X_i＝X_i min+k_i*Δ_i，其中 ${\mathrm{\&Delta;}}_i=\frac{X_{i\max }-X_{i\min }}{\mathrm{dvi}},$ X_i∈[X_i min，X_i max]，(i＝1，2，…，n)，便可得到各优化参数X_i，用于计算个体适应值。

选择：从遗传算法整个选择策略来讲，精英选择策略是使群体收敛到最优解的一种保障，故本文采用精英选择策略。如果当前群体的最佳个体适应值大于下一代群体最佳个体的适应值，则将当前群体最佳个体或适应值大于下一代最佳个体适应值的多个个体直接复制到下一代，随机替换或者替换适应值最差者。

杂交：设X_i ^t和X_j ^t分别为第t代要进行交叉的两个个体，则交叉后产生的子代个体为：

$\left\{\begin{array}{c}{X}_i^{t+1}=X_i^t+{\mathrm{\&tau;}}_1(X_i^t-X_j^t)\\ X_j^{t+1}=X_j^t+{\mathrm{\&tau;}}_2(X_i^t-X_j^t)\end{array}\right.$

其中τ₁和τ₂为[0，1]上均匀分布的随机数。从加快收敛速度和全局搜索性能两方面考虑，受自然界中家庭内兄弟间竞争现象的启发，在算法中加入小范围竞争、择优操作。方法是：对将要进行交叉的父母对A、B进行n(3-6)次交叉，生成2n个体，选出最优的两个个体进入子代。该方法的实质是在相同父母的情况下，预先加入兄弟间的小范围的竞争择优机制。变异：采用高斯变异法，按照如下的变异概率对选中的个体进行变异操作：

$P_m=\left\{\begin{array}{cc}k(f_{\max }-f)/(f_{\max }-f_{\mathrm{avg}})& 0<(f_{\max }-f)/(f_{\max }-f_{\mathrm{avg}})\&le;1\\ 0.75& (f_{\max }-f)/(f_{\max }-f_{\mathrm{avg}})<0\end{array}\right.,k\&Element;\left[\mathrm{0,1}\right]$

其中，Pm为个体的变异概率；fmax为种群的最大适应值；favg为种群的平均适应值；f为个体适应值。

线性反演方法：由于遗传算法优化速度慢，优化结果可能在最优值附近而非真实的最优值，本发明采用实数编码遗传算法与线性反演方法相结合的方法予以解决，进一步对各个体(PML参数)再进行局部优化，把进化一定代数后的种群中的个体作为最速下降反演方法的初始模型，将线性反演方法迭代一定次数后的结果返回种群中，更新种群。本发明选取最速下降法作为补充优化方法，迭代一定次数，即可快速得到本代种群的PML参数局部最优解。

由以上三个步骤完成实数遗传算法后，当目标函数小于给定误差限时，遗传算法结束，再进行交错网格的FDTD(2，2M)迭代计算三维空间的电磁场步进计算过程，便完成了超宽带电磁法的三维波场模拟。

1.本发明与现有模拟方法相比，其核心在于：(1)使用分段线性递归卷积法计算CPML方法中的卷积项，取代CPML法中对前后两次迭代时间段内的电磁场以常数进行计算的方法；(2)对各个空间微分项(含PCPML的人工边界部分和计算空间区域部分)均采用2M阶高阶差分进行计算；(3)使用实数遗传算法结合线性反演方法，获得最优PML参数，用于FDTD模拟。基于上述特点，本发明具有如下优点：

(1)能够较准确地模拟超宽带电磁波的传播，通过调整参数因子不仅能象一般的PML方法那样吸收高频的波，也能吸收低频的波，同时还能较高精度地模拟超宽带电磁波的长时反射，符合超宽带电磁波对模拟的基本要求；

(2)吸收边界处理效率高，由于其递推公式不再与模拟物体的类型有关，从而它可模拟各种复杂三维介质模型，不仅适用于空气介质，也可用于任意复杂的兼有有耗介质和空气介质的三维模型，以及各种色散、非线性介质等，有较宽的模型适用范围；

(3)所采用的用分段线性递归卷积法计算卷积项和2M阶高阶差分方法大大提高了计算精度、降低了网格色散和边界反射造成的反射误差；

(4)采用改进的遗传优化算法与线性反演方法相结合的办法获取吸收边界参数，让计算机智能化选择方法代替人为盲目设定某个具体值，从而降低了模拟参数选择的难度，吸收效果也大大提高；

(5)所采用的高阶时域有限差分计算方式能较高精度地模拟空间电磁场，并能以直观的方式给出各时间步进的电磁波在复杂介质中的传播。

下面以实例说明本发明与现有模拟方法相比所具有的效果：

(1)图3为用CPML-FDTD(2，4)计算的y＝0m断面上的Ex分量在39.75ns时的分布图，图4为用PCPML-FDTD(2，4)计算的y＝0m断面上的Ex分量在39.75ns时的分布图，图3中Ex分量的长时边界反射分量误差污染了波场数据图，对比图3、图4可见，超宽带电磁波的PCPML吸收边界方法的长时反射的吸收效果明显好于CPML吸收边界方法，仅在空气部分的两侧边界处有弱的边界反射。

(2)图5为PCPML-FDTD(2，4)方法在原点处的Ex分量波形曲线，图6为CPML-FDTD(2，4)方法在原点处的Ex分量波形曲线，图7为解析解在原点处的Ex分量波形曲线，可见直达波和地面波同解析解，而在39ns左右时，用CPML吸收边界处理方法模拟结果出现多余峰值波，这就是超宽带的长时反射误差；对比两种方法可看出，PCPML方法对长时反射的良好吸收效果与计算精度的大幅提高，与解析解很接近，仅在27纳秒处有弱的误差反射，而CPML的长时反射吸收效果很不好，而且精度没有PCPML-FDTD(2，4)方法的精度高，充分说明了本发明的优势。

(3)图8-图18为模拟超宽带电磁波穿墙探测生命体模型的例子。图8-图9为含生存者的地震倒塌房屋模型，设计算区域为[-0.3m，1.8m]×[-0.9m，0.9m]×[-1.6m，0.3m]，人体位于区域[0.6m，0.9m]×[-0.5m，0.5m]×[-0.4m，0m]，相对介电常数为12，电导率为0.5S/m。如图(8)、图(9)分别对应的xz垂直断面(y＝0m)和xy水平剖面(z＝-0.2m)。图(8)中z＜0m部分为相对介电常数为6、电导率为0.05S/m的地下介质，x在0m～0.3m之间部分为墙体(相对介电常数为4，电导率为0.04S/m)，人体上方为倒塌的天花板(相对介电常数为8，电导率为0.0125S/m)，厚度为0.2m。假定发射天线是4cm×4cm×1cm的均匀体天线，中心位于墙体上点(0m，0m，-0.2m)处，天线上馈入y方向的均匀截断三正弦电流信号，脉冲宽度为1.25ns。遗传算法的PML参数初始值设定为(层数8～17层，层电导率参数变化范围为0.01～1S/m，因子κ_i的范围均取为1～20，因子α_i均取为(i＝x，y，z)，其中λ的变化范围为0.01～2.0，其它参数按下面的具体实施步骤中进行设置。

用本发明PCPML-FDTD(2，4)方法迭代计算三维空间各节点处的场值。图10-图17为图8-图9模型在2.5ns、5ns、6ns和8ns时y＝0m和z＝-0.2m对应的两个断面的波场快照。图片显示，在2.5ns时，脉冲电磁波在墙外空气、地下和障碍墙体内传播，尚未进入房屋内；在5ns时，电磁波穿过障碍墙体，进入人体，同时，人体的回波进入障碍墙体，跟在内墙面的反射波之后向墙体外方向传播；在6ns时，接收点开始收到人体的回波信号，而房屋内的电磁波的一部分继续在人体内部传播，另一部分进入倒塌的天花板内，而天花板的反射波又进入人体；8ns时，电磁波透射出天花板，向房屋外空中传播，而在房屋和墙体内，由天花板引起的反射波占优势地位，并与墙体和人体之间相互作用使得波场分布极为复杂，回波信号的波形也非常复杂。

假定人体微动区域为以x＝0.6m为中心、在y方向[-0.1m，0.1m]和z方向[-0.3m，-0.1m]区域沿x方向向左(扩张微动)或向右(收缩微动)运动。根据上述方法和模拟参数，对人体在常态(dL＝0cm)、收缩4cm(dL＝4cm)和扩张4cm(dL＝-4cm)3种情况下进行仿真，以探讨连续发射超宽带电磁波探测人体微动情况。图18为点(0m，0m，-0.2m)处的Hy分量的波形曲线，其中在5ns前，人体回波还没到达接收点，波形重合；在6～7ns间，开始出现人体微动异常，在7ns～8.5ns间，回波信号主要为人体的反射波，人体微动导致空间介质电性分布变化，使得接收信号波形不同，人体扩张微动时，波的旅行时间短，反之则旅行时间长，与实际情况是一致的；而在8.5ns后，倒塌天花板和墙体内的回波和它们的多次波等对回波信号的影响占优势地位，波形差异相对减弱，这与上述的波场快照结果是一致的。

附图说明

图1为场函数A(t)的分段线性近似(实线)与CPML计算方法(虚线)的对比图。

图2为卷积项内场函数A(nΔt-τ)的分段线性近似计算方法的示意图。

图3为用CPML-FDTD(2，4)计算的y＝0m断面上的Ex分量在39.75ns时的分布图。

图4为用PCPML-FDTD(2，4)计算的y＝0m断面上的Ex分量在39.75ns时的分布图。

图5为PCPML-FDTD(2，4)方法在原点处的Ex分量波形曲线。

图6为CPML-FDTD(2，4)方法在原点处的Ex分量波形曲线。

图7为解析解在原点处的Ex分量波形曲线。

图8为含生存者的地震倒塌房屋模型的xz垂直断面(y＝0m)。

图9为含生存者的地震倒塌房屋模型的xy水平剖面(z＝-0.2m)。

图10为图8、图9所示模型在2.5ns时y＝0m对应的断面的波场快照。

图11为图8、图9所示模型在2.5ns时z＝-0.2m对应的断面的波场快照。

图12为图8、图9所示模型在5ns时y＝0m对应的断面的波场快照。

图13为图8、图9所示模型在5ns时z＝-0.2m对应的断面的波场快照。

图14为图8、图9所示模型在6ns时y＝0m对应的断面的波场快照。

图15为图8、图9所示模型在6ns时z＝-0.2m对应的断面的波场快照。

图16为图8、图9所示模型在8ns时y＝0m对应的断面的波场快照。

图17为图8、图9所示模型在8ns时z＝-0.2m对应的断面的波场快照。

图18为点(0m，0m，-0.2m)处的Hy分量的波形曲线。

具体实施方式

本发明是基于递归计算方式的分段线性离散近似与高阶差分离散CPML的卷积项和PML参数的自适应获取方法。具体包括如下步骤：

(1)给定三维模型的网格化的地电参数、背景模型的地电参数、电偶极子发射天线参数、激励脉冲信号参数、时间步长、最大时间步数Nt(或最大传播时间)、PML层数、各PCPML参数搜索范围及遗传算法计算参数；

(2)根据遗传算法优化PML参数时所用的均匀半空间模型或均匀层状模型，计算遗传优化的目标函数中的参考响应值。遗传优化的目标函数为： $\mathrm{Error}=20\underset{P}{\mathrm{\&Sigma;}}\frac{1}{\mathrm{Np}}\underset{y=1}{\overset{3}{\mathrm{\&Sigma;}}}\underset{t=0}{\overset{\mathrm{Nt}}{\mathrm{\&Sigma;}}}{\log }_{10}\frac{|E_y\left(t\right)-E_y^{\mathrm{ref}}\left(t\right)|}{\left|E_{y\_\max }^{\mathrm{ref}}\left(t\right)\right|}$ 的倒数为遗传优化的目标函数，其中，Nt为最大时间步数，y为空间三个方向指标，E_y(t)为边界面上点P的场Ey信号，E_y ^ref(t)为P点的场Ey信号的参考值，E_{y_max} ^ref(t)为P点的场Ey参考值的最大值，

表示对所有边界面上的点P求和，Np表示边界节点个数；

(3)使用改进的实码混合遗传算法并结合线性反演算法，优化计算PML参数，下面是遗传算法优化计算的步骤：

(a)初始种群生成：把各待优化参数的搜索范围均匀分割，在各小区间端点序列中随机生成一个初始个体。根据各参数的给定参数范围[X_i min，X_i max]和离散个数进行均匀分割。在[0，Ni]中产生一个随机整数k_i，Ni为第i个参数所在区间被分割的个数，k_i即为个体的元素。把待优化参数对应的k_i依次连在一起构成问题解的编码形式为(k₁，k₂，…，k_n)。在计算适应值时将k_i做映射X_i＝X_i min+k_i*Δ_i，其中 ${\mathrm{\&Delta;}}_i=\frac{X_{i\max }-X_{i\min }}{\mathrm{dvi}},$ X_i∈[X_i min，X_i max]，(i＝1，2，…，n)，解编的X_i用于计算个体适应值；

(b)选择计算：采用精英选择策略，如果当前群体的最佳个体适应值大于下一代群体最佳个体的适应值，则将当前群体最佳个体或适应值大于下一代最佳个体适应值的多个个体直接复制到下一代，随机替换或者替换适应值最差者；

(c)交叉计算：设X_i ^t和X_j ^t分别为第t代要进行交叉的两个个体，则交叉后产生的子代个体为： $\left\{\begin{array}{c}{X}_i^{t+1}=X_i^t+{\mathrm{\&tau;}}_1(X_i^t-X_j^t)\\ X_j^{t+1}=X_j^t+{\mathrm{\&tau;}}_2(X_i^t-X_j^t)\end{array}\right.,$ 其中τ₁和τ₂为[0，1]上均匀分布的随机数，对将要进行交叉的父母对A、B进行n(3-6)次交叉，生成2n个个体，选出最优的两个个体进入子代依次对选中交叉的父母对都进行6次交叉，在产生子代和对应的父母对中进行优选，选出两个最优的个体进入子代；

(d)变异计算：分别计算种群的最大适应值fmax、种群的平均适应值，设f为个体适应值，按如下的变异概率Pm对选中的个体进行变异操作：

$P_m=\left\{\begin{array}{cc}k(f_{\max }-f)/(f_{\max }-f_{\mathrm{avg}})& 0<(f_{\max }-f)/(f_{\max }-f_{\mathrm{avg}})\&le;1\\ 0.75& (f_{\max }-f)/(f_{\max }-f_{\mathrm{avg}})<0\end{array}\right.$

其中，k取为0～1的某个值；

(e)线性反演：当遗传到初始设定的代数或保持较大交叉概率、变异概率达20代以上时，则退出遗传，得出初步优化结果。再把得到优化后PML参数结果作为最速下降方法的初始模型X⁰，迭代反演一定次数(取为10次左右)后的结果返回种群中，更新种群，得到最终优化的PML参数。具体迭代过程如下：首先计算 $S^k=-\&dtri;f\left(X^k\right);$ 其次，循环计算(k＝0～10)如下式：计算第k次迭代的雅可比矩阵

更新PML参数

k＝k+1，参数β取为1/2或1/4；循环计算完成后，线性反演得到的结果更新了本代遗传优化的种群，进行下一代的遗传优化反演。当遗传反演的目标函数小于某个给定误差范围后，便得到PML参数最终优化结果，退出遗传优化过程。令n＝0，进入步骤(4)；

(4)根据优化得到PML参数和FDTD的迭代初值一系列参数，用分段线性递归卷积法计算如下吸收边界迭代公式(以x方向分量为例)：

${\mathrm{\&epsiv;}}_r{\mathrm{\&epsiv;}}_0\frac{{\&PartialD;E}_x}{\&PartialD;t}+{\mathrm{\&sigma;E}}_x=\frac{1}{k_y}\frac{{\&PartialD;H}_z}{\&PartialD;y}-\frac{1}{k_z}\frac{{\&PartialD;H}_y}{\&PartialD;z}+{\mathrm{\&xi;}}_y\left(t\right)*\frac{{\&PartialD;H}_z}{\&PartialD;y}-{\mathrm{\&xi;}}_z\left(t\right)*\frac{{\&PartialD;H}_y}{\&PartialD;z}$

里的卷积项，其中的差分计算采用空间2M阶的高阶差分方法：

$f^{\&prime;}\left(x\right)=\frac{1}{\mathrm{\&Delta;L}}\underset{i=-M}{\overset{M-1}{\mathrm{\&Sigma;}}}c\left(i\right)f(x+i+0.5)$ $c\left(i\right)=\frac{{(-1)}^i{[(2M-1)!!]}^2}{2{(i+0.5)}^2(2M-2-2i)!!(2M+2i)!!}$

若记第n个时间步的电磁场偏导数项为A(t)，则卷积项

可按如下公式迭代计算：

其中uⁿ按如下递归方式计算：

$u^n=\left\{\begin{array}{cc}{A}^{n-1}-A^{n-2}+u^{n-1}{e}^{-\mathrm{s\&Delta;t}}& n\&GreaterEqual;2\\ 0& n<2\end{array}\right.$

另一个卷积项也同样进行，在非PML区域的三维空间中，按如下迭代公式计算磁网格下的x方向分量电场：

$E_x^{n+1}(i,j+0.5,k+0.5)=p1E_x^{n+1}(i,j+0.5,k+0.5)+p2[27H_z^{n+0.5}(i,j+1,k+$

$0.5)-27H_z^{n+0.5}(i,j,k+0.5)-H_z^{n+0.5}(i,j+2,k+0.5)+H_z^{n+0.5}(i,j+0.5,k-1)]$

$+p3{\mathrm{\&psi;}}_{\mathrm{exy}}^{n+0.5}(i,j+0.5,k+0.5)-p3{\mathrm{\&psi;}}_{\mathrm{exz}}^{n+0.5}(i,j+0.5,k+0.5)$

其中， $p1=\frac{2\mathrm{\&epsiv;}-\mathrm{\&sigma;\&Delta;t}}{2\mathrm{\&epsiv;}+\mathrm{\&sigma;\&Delta;t}},$ $p2=\frac{\mathrm{\&Delta;t}}{12(3\mathrm{\&epsiv;}+\mathrm{\&sigma;\&Delta;t})k_y\mathrm{\&Delta;y}},$ $p3=\frac{2\mathrm{\&Delta;t}}{2\mathrm{\&epsiv;}+\mathrm{\&sigma;\&Delta;t}},$

${\mathrm{\&psi;}}_{\mathrm{exy}}^{n+0.5}(i,j+0.5,k+0.5)=B_y{\mathrm{\&psi;}}_{\mathrm{exy}}^{n-0.5}(i,j+0.5,k+0.5)+A_y[27H_z^{n+0.5}(i,j+1,k+,$

$0.5)-27H_z^{n+0.5}(i,j,k+0.5)-H_z^{n+0.5}(i,j+2,k+0.5)+H_z^{n+0.5}(i,j-1,k+0.5)]$

${\mathrm{\&psi;}}_{\mathrm{exz}}^{n+0.5}(i,j+0.5,k+0.5)=B_2{\mathrm{\&psi;}}_{\mathrm{exy}}^{n-0.5}(i,j+0.5,k+0.5)+A_z[27H_z^{n+0.5}(i,j+0.5,k,$

$+1)-27H_z^{n+0.5}(i,j+0.5,k)-H_z^{n+0.5}(i,j+0.5,k+2)+H_z^{n+0.5}(i,j+0.5,k-1)]$

$A_s=\frac{a_s}{24\mathrm{\&Delta;s}},$ $B_s=\exp [-(\frac{{\mathrm{\&sigma;}}_s}{k_s}+{\mathrm{\&alpha;}}_s)\frac{\mathrm{\&Delta;t}}{{\mathrm{\&epsiv;}}_0}],$ (s＝x，y，z)，

$a_i=\frac{{\mathrm{\&sigma;}}_i}{{\mathrm{\&sigma;}}_i{k}_i+k_i^2{\mathrm{\&alpha;}}_i}\left\{\exp \right[-(\frac{{\mathrm{\&sigma;}}_i}{k_i}+{\mathrm{\&alpha;}}_i)\frac{\mathrm{\&Delta;t}}{{\mathrm{\&epsiv;}}_0}]-1.0\},$ (i＝x，y，z)。i，j，k为磁网格的空间采样节点编号，n为时间采样序号，Δs为s方向上网格步长，Δt为时间步长。根据以上公式计算三维空间的电场x分量；

(5)按第(4)步骤方法依次迭代计算y、z方向电场分量；

(6)按FDTD(2，2M)方法迭代计算第n+1/2个时间步的三个磁场分量；

(7)取n＝n+1，重复第(4)、(5)、(6)步骤，直到最大时间步Nt，完成FDTD迭代计算；

(8)保存电磁场数据，按常规成图方法得到相应的波场快照。

Claims (1)

1.一种超宽带电磁波的模拟方法，其特征在于该方法是利用空间2M阶的高阶FDTD方法模拟电偶极子天线激励的超宽带电磁波的传播，其吸收边界处理采用分段线性卷积方式卷积完全匹配层吸收边界方法即PCPML方法，该方法包括以下步骤：

(1)输入初始参数：给定三维模型的网格化的地电参数、背景模型的地电参数、发射电偶极天线参数、激励脉冲信号参数、时间步长、最大时间步数或模拟电磁波的最大传播时间、PML层数、PCPML吸收参数搜索范围及遗传优化算法的计算参数；

(2)计算优化PML参数的目标函数参数：计算电偶极子发射天线激发的超宽带电磁脉冲的电磁响应，作为遗传优化的目标函数中的参考响应值；

(3)使用实码混合遗传算法并结合线性反演算法、以及采用精英选择策略并结合线性反演的局部优化种群过程，优化计算PML参数，该PML参数包括各PML层内的方向电导率(σ_i)，对表面波的吸收的可调因子(κ)和对FDTD难以吸收的低频波的可调因子(α_i)；当目标函数值小于给定的误差限时，遗传算法结束；此时令时间步n＝0，开始FDTD迭代；

(4)用FDTD方法迭代计算第n个时间步的三维空间电场值：根据优化得到PML参数，用空间2M阶的高阶差分方法和分段线性递归卷积吸收边界方法得到的FDTD公式依次计算第n个时间步的三个电场分量Ex、Ey和Ez；

(5)用FDTD方法迭代计算第n+1/2个时间步的三维空间磁场值：按FDTD(2，2M)方法迭代计算第n+1/2个时间步的三个磁场分量Hx、Hy和Hz；

(6)取n＝n+1，重复第(4)、(5)步骤，直到最大时间步Nt，完成FDTD迭代计算；

(7)保存电磁场数据，按常规成图方法得到相应的波场快照。

Images