一种地下各向同性介质球电磁散射的计算方法与流程

本发明属于电磁散射理论计算领域，具体涉及一种地下各向同性介质球电磁散射的计算方法。

背景技术

任何一个电磁问题的分析计算通常都需要求解麦克斯韦方程组和边界值。经过多年的发展，一些经典的电磁散射问题已经得出准确的解析。而实际的物理环境是非常复杂的，需要采用数值计算得到具体的电磁波传播特性。而采用解析方法和数值计算相结合的方法是常用的研究方法。根据研究目标和工作波长之间的关系，电磁场的数值计算方法可以分为低频数值方法和高频近似方法两类。低频数值法有：矩量法(mom)、有限元法(fem)、时域有限差分法(fdtd)等，这类数值方法能准确的解决几何结构和组成材料复杂的电磁问题，但低频数值方法受限于计算机物理内存、cpu的工作速度。高频近似方法主要有：几何光学法(go)、物理光学法(po)、几何绕射理论(gto)等。这类方法具有物理概念清晰、快速便捷的特点。

自然界中地表下分布着各种各样的媒质，利用低频电磁波的大穿透能力可以获取地下土质信息，但也会受到地下掩埋物的影响，因此对地下三维目标电磁散射的研究受到广泛关注。对自由空间中的介质体解析解的研究主要有各向同性介质，单轴、双轴、铁氧体介质、等离子体等。对地下三维目标的解析方法的研究相对比较少。对地下导体球电磁散射解析方法的研究，如专利号cn201710160218.7的发明专利公开的计算地下导体球电磁散射的解析方法，步骤1.利用均匀平面波对无界均匀媒质中的传播特性得到反射系数r和透射系数t,以及在z＝-d处的反射电场e^r(z)和透射电场e^t(z)；步骤2.将透射波作为导体球的入射波,然后根据本征矢量和平面波因子乘积的解析,将透射波展开成地下各向同性的球矢量波函数,利用导体球表面的边界条件得到散射电场展开系数；步骤3.导体球的散射波作为入射波对地平面进行反射和透射,散射波的反射波又一次作为导体球的入射波,循环反复；在此循环过程中根据分界面的边界条件得到总体坐标系下的透射场的参数值与局部坐标系下反射场的参数值并利用球矢量波函数的加法定理来得到在总体坐标系下反射场的参数值。对双各向异性介质球解析解的研究，如专利号cn201310156058.0的发明专利公开的内容：步骤1.利用无源麦克斯韦方程组和双各向异性媒介的本征方程推导磁感应强度b的微分方程；步骤2.将微分方程中和b相关的因子以球矢量波函数的形式表达,然后利用球矢量波函数m,n的正交性质得出一个含参的矩阵方程,先利用矩阵方程满足非零解的条件计算出该矩阵方程的参数,再将参数代回到含参数的矩阵方程中得到矩阵方程的非零解；步骤3.构造一个新的函数,用新函数重新表示磁感应强度b,进而求出介质球内部的电磁场,然后把介质球内的电磁场和球外的入射电磁场、散射电磁场代入到边界条件中,得出散射矩阵。现有技术中对于地下各向同性介质球的散射特性的解析方法几乎没有。

技术实现要素：

本发明的目的在于针对地下三维目标电磁散射解析解理论的不足，提出一种地下各向同性介质球电磁散射的计算方法，该方法基于递归t矩阵算法，具有快速、高效地求解地下三维目标的电磁散射特性的优点。

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

一种地下各向同性介质球电磁散射的计算方法，其具体步骤包括：

步骤1.一均匀平面波从自由空间沿坐标系正z方向垂直入射到距离坐标系原点d处，即z＝-d的无限大的有耗介质分界平面上，利用均匀平面波在无界均匀媒质中的传播特性，根据边界条件建立方程组得出z＝-d处的反射电场和透射电场

步骤2.根据本征矢量和平面波因子乘积的解析展开式，将各向同性介质球的入射波的电磁场用各向同性球矢量波函数展开，然后利用均匀平面波的传播特性及各向同性介质球表面的边界条件，得出介质球的反射波、透射波的电场和磁场

步骤3.入射到无限大的有耗介质分界平面上进行反射和透射，而反射波又一次作为入射波入射到各向同性介质球上，以此循环反复；在此循环过程中根据介质平面的镜像原理得出反射场和透射场的表达式，再利用边界条件联立方程组得出全局坐标系下的透射场与镜像坐标系下的反射场的具体值。

步骤4.利用球矢量波函数的叠加原理求得在全局坐标系下对应反射场然后采用递归t矩阵算法和传统理论方法得出矩阵变换系数，得出的具体值。

进一步，通过最后联立方程组得出某点的电场值，在解析过程中建立数据库，提高运算效率与精度。

进一步，步骤1中在z＝-d处的反射电场和透射电场的具体计算如下：

在z＜-d空间是自由空间，参数ε0,μ0,σ0；自由空间是无损耗空间其σ0＝0，ε0＝ε0εr0,μ0＝μ0μr0，其中σ0、ε0、εr0、μ0、μr0分别表示自由空间的电导率、介电常数、相对介电常数、磁导率、相对磁导率，μ0＝4π×10^-7h/m，εr0＝1,μr0＝1；在此空间一均匀平面波从自由空间垂直入射到z＝-d的无限大的有耗介质分界平面上，该入射波是均匀平面波，其电场的幅度为1，其沿正z方向垂直入射，其电场的极化方向为正x方向；故在自由空间中沿正z方向传播的平面波的入射电场和磁场可以表示为：

其中表示自由空间中的本征阻抗、波数；黑色加粗字体表示该值是矢量，表示直角坐标系的单位矢量，i表示虚数单位；在z＝-d的无限大的有耗介质分界平面上，利用均匀平面波在无界均匀媒质中的传播特性，得出反射波和透射波的电场、磁场分别为：

透射波所在的空间即z＞-d，空间介质是媒质1，该空间是有耗介质空间，电导率σ1≠0，ε1＝ε0εr1,μ1＝μ0μr1，其中σ1、ε1、εr1、μ1、μr1分别表示媒质1的电导率、介电常数、相对介电常数、磁导率、相对磁导率，；该空间掩埋有一半径为a的各向同性介质球，埋藏深度为d，且d＞a，各向同性介质球的参数为ε2,μ2,σ2，其中其中σ2、ε2、εr2、μ2、μr2分别表示介质球的电导率、介电常数、相对介电常数、磁导率、相对磁导率；r,t为反射波和透射波的反射和透射系数；η1、k1为媒质1的本征阻抗、波数，η2、k2为介质球的本征阻抗、波数；0表示自由空间，1表示媒质1空间，2表示各向同性介质球；本征阻抗与波数表达是为：

根据边界条件建立方程组，在z＝-d的分界平面上，应有

得出反射系数r和透射系数t的表达式分别为：

将(9)(10)带入式(3)(4)得到在z＝-d处的反射波和透射波的电场值

进一步，步骤2中介质球界面的反射波、透射波的电场和磁场的具体计算如下：

各向同性介质球的入射波根据本征矢量和平面波因子乘积的解析展开式，将入射波的电磁场用各向同性球矢量波函数展开的具体表达式为：

各向同性介质球的入射波沿着正z方向入射，其电场的极化方向是正x方向，其电场的幅度为t，其中

经各向同性介质球面散射后，利用均匀平面波在无界均匀媒质中的传播特性，反射波和透射波的电场与磁场分别为：

其中，表示n从1到∞，m从-n到+n的求和，表示反射电场与透射电场的反射系数和透射系数；表示球矢量波函数，l＝1，l＝3分别表示矢量波函数由第一类、第三类球贝塞尔函数构成，表示球坐标系中的一个矢量；通过边界条件，在r＝a时，电磁场在切向边界上连续，即有：

联立方程组得反射电场和透射电场的展开系数

其中jn(kr)为球贝塞尔函数，为第一类球汉克尔函数,p的具体表达式：

将(17a-d)带入(15a-d)得出各向同性介质球界面的反射波与透射波的电场和磁场的具体值

进一步，步骤3中全局坐标系下的透射场与镜像坐标系下的反射场的具体值计算如下：

全局坐标系：球坐标系中的三个坐标变量r,θ,φ；

局部坐标系：球坐标系中的三个坐标变量rimg,θimg,φimg；

在球坐标系中，过空间任意一点的三个相互正交的坐标单位矢量分别是r，θ，增加的方向，且遵循右手螺旋法则，它们与之间的变换关系为：

经各向同性介质球散射，其反射波作为入射波入射到分界平面上，利用介质平面的镜像原理其反射波和透射波的电场和磁场可表示为：

其中反射波的电场和磁场用局部坐标系表示，透射波的电场和磁场用全局坐标系表示，根据边界条件及p点θ＝π，θimg＝0，r＝rimg＝d，可以计算得到展开系数

其中mnm(θ,φ)，nnm(θ,φ)的具体表达式为：

联立(20a-d)(21)得出全局坐标系下的透射场与镜像坐标系下的反射场的具体值。

进一步，步骤4中利用球矢量波函数的叠加原理来求得在全局坐标系下对应反射场然后采用递归t矩阵算法和传统理论方法得出矩阵变换系数，并建立数据库，提高计算效率与精度，最后得出的具体值，计算如下：

利用球矢量波函数的叠加原理将局部坐标系表示的反射电场转换成全局坐标系表示的反射电场全局坐标系下的反射电场的具体表达式为：

三维矢量波函数加法定理下的之间的关系可以得到具体的表达式：

式中表示v从1到∞，u是从-v到+v的累加。是第一类三维矢量波函数；称为坐标系变换系数；

联立(20a)，(23)，(24)，(25)式，得出矩阵形式表示如下：

使用传统的方法计算矩阵变换系数首先由平面波的展开公式：

其中θk,φk表示k的指向方位角，θ,φ表示r的指向方位角，jl表示l阶球bessel，ylm表示归一化的球贝塞尔函数

由r＝r'+r”得

e^ik·r＝e^ik·r'·e^ik·r”(29)

联立(27)(28)(29)式根据ylm(θk,φk)的正交归一性得

其中是wigner3-j符号,

由此得

其中zl(x)等于四类贝塞尔函数，

由以上推导可以看出，此方法计算坐标系变换系数的运算量是相当大的，每个βl'm',lm,αl'm',lm需要对l”进行求和，而且wigner3-j符号依赖与多种因素，使计算量增加。

采用递归t矩阵算法进行求解坐标系变换系数根据经典的mie散射理论可以得出单个散射球的t矩阵表示为

根据矩阵变换系数建立并矢坐标变换矩阵

其中的下角标ji表示坐标系i到坐标系j的变换；的上角标(j),(h)表示坐标变换系数中基函数的选取为球的bessel函数或第一类球hankel函数；递推方法的初值：

故球谐的加法公式可以写成：

根据球面波函数的性质

其中是连带的勒让德函数，联立(38)(39)(40)(41)得

再利用连带的勒让德函数、贝塞尔函数、球矢量波函数得：

其中为系数

根据(43)式发现，当给定一定的n，如果知道nn的一系列的转移矩阵元素值，那么就可以知道一系列的(n+1,n+1)的转移矩阵元素，所以(43)式提供的是一条沿m,n的对角线的递推路径。

根据(44)式发现当知道(n+1,m),(n,m)的波函数的加法定理系数后，通过递推关系可以求得((n+1,m)系数，在这个式子里m是不变的，它的递归路线是向上的。

联立(37)(42)(43)(44)得

其中

在计算矩阵变换系数时，建立阶乘数据库，把递推的数据存放在数据库中，根据需求调用数据库中的数值，这样计算方式比直接计算方法大大节省时间，提高了计算效率。最后联立(23)(26)(43)(48)得出全局坐标系下的反射电场的值。

进一步，步骤1计算过程中自由空间的相对介电常数εr0＝1,相对磁导率μr0＝1，媒质1是土壤媒质，土壤媒质的相对介电常数εr1＝5或εr1＝30,相对磁导率μr1＝1,电导率σ1＝0.002或σ1＝0.1，各向同性介质球的相对介电常数εr2＝10,相对磁导率μr2＝1。

进一步，通过三维电磁仿真软件feko，设置近场点，得出的值即为所求的电场值，然后通过fortran编程语言根据求解的电场的具体值，再通过迭代循环，建立方程组，最后得到在p点处的电场值即程序所得到的电场值。

本发明的有益效果如下：

采用基于递归t矩阵算法与传统理论方法求解地下各向同性介质球的电磁散射，并将递推数据建立数据库保存起来，提高数据的运算速度，进而提高计算效率。通过改变频率，媒质1的电导率、相对介电常数，埋藏深度，将在p点的电场的仿真值以及程序所得到的电场值进行对比来验证该方法的正确性，进而得出参数对电场值的影响。本发明适用于快速、高效地求解地下三维目标的电磁散射特性。

附图说明

图1是本发明给出地下介质球目标电磁散射的几何图。

图2是本发明给出全局坐标系和局部坐标系的示意图。

图3是本发明给出加法定理中r，r′,r″之间关系的示意图。

图4是本发明给出三维矢量波函数加法定理递推公式示意图。

图5是本发明给出εr1＝5，σ1＝0.1，εr2＝10，a＝5m，d＝25m时，不同频率的仿真和程序所得电场值的对比图。

图6是本发明给εr1＝5，σ1＝0.1，εr2＝10，a＝5m，d＝25m和εr1＝5，

σ1＝0.002，εr2＝10，a＝5m，d＝25m，不同电导率下仿真和程序所得电场值的对比图。

图7是本发明给出εr1＝30，σ1＝0.1，εr2＝10，a＝5m，d＝25m和εr1＝5，

σ1＝0.1，εr2＝10，a＝5m，d＝25m时，不同相对介电常数下仿真和程序所得电场值的对比图。

图8是本发明给出εr1＝5，σ1＝0.002，εr2＝10，a＝5m，d＝25m和εr1＝5，

σ1＝0.002，εr2＝10，a＝5m，d＝10m，不同埋藏深度下仿真和程序所得电场值的对比图。

图9是本发明给出εr1＝5，σ1＝0.002，εr2＝10，a＝5m，d＝25m和地下无掩埋物时，仿真和程序所得电场值的对比图。

具体实施方式

下面结合具体实施例来对本发明进行进一步说明，但并不将本发明局限于这些具体实施方式。本领域技术人员应该认识到，本发明涵盖了权利要求书范围内所可能包括的所有备选方案、改进方案和等效方案。

一种基于递归t矩阵算法的地下各向同性介质球电磁散射的计算方法，其具体步骤包括：

步骤1.一均匀平面波从自由空间沿坐标系正z方向垂直入射到距离坐标系原点d处，即z＝-d的无限大的有耗介质分界平面上，利用均匀平面波在无界均匀媒质中的传播特性，根据边界条件建立方程组得出z＝-d处的反射电场和透射电场

步骤2.根据本征矢量和平面波因子乘积的解析展开式，将各向同性介质球的入射波的电磁场用各向同性球矢量波函数展开，然后利用均匀平面波的传播特性及各向同性介质球表面的边界条件，得出介质球的反射波、透射波的电场和磁场

步骤3.入射到无限大的有耗介质分界平面上进行反射和透射，而反射波又一次作为入射波入射到各向同性介质球上，以此循环反复；在此循环过程中根据介质平面的镜像原理得出反射场和透射场的表达式，再利用边界条件联立方程组得出全局坐标系下的透射场与镜像坐标系下的反射场的具体值。

步骤4.利用球矢量波函数的叠加原理来求得在全局坐标系下对应反射场然后采用递归t矩阵算法和传统理论方法得出矩阵变换系数，得出的具体值，最后联立方程组得出某点的电场值。

本实施例步骤1中在z＝-d处的反射电场和透射电场的具体计算如下：

研究目标的几何模型参见图1，在z＜-d空间是自由空间，参数ε0,μ0,σ0，自由空间是无损耗空间其σ0＝0，ε0＝ε0εr0,μ0＝μ0μr0，其中σ0、ε0、εr0、μ0、μr0分别表示自由空间的电导率、介电常数、相对介电常数、磁导率、相对磁导率，μ0＝4π×10^-7h/m，εr0＝1,μr0＝1；在此空间一均匀平面波从自由空间垂直入射到z＝-d的无限大的有耗介质分界平面上，该入射波是均匀平面波，其电场的幅度为1，其沿正z方向垂直入射，其电场的极化方向为正x方向；故在自由空间中沿正z方向传播的平面波的入射电场和磁场表示为：

其中表示自由空间中的本征阻抗、波数；黑色加粗字体表示该值是矢量，表示直角坐标系的单位矢量，i表示虚数单位；在z＝-d的无限大的有耗介质分界平面上，利用均匀平面波在无界均匀媒质中的传播特性，得出反射波和透射波的电场、磁场分别为：

透射波所在的空间即z＞-d，空间介质是媒质1，该空间是有耗介质空间，电导率σ1≠0，ε1＝ε0εr1,μ1＝μ0μr1，其中σ1、ε1、εr1、μ1、μr1分别表示媒质1的电导率、介电常数、相对介电常数、磁导率、相对磁导率，；该空间掩埋有一半径为a的各向同性介质球，埋藏深度为d，且d＞a，各向同性介质球的参数为ε2,μ2,σ2，其中其中σ2、ε2、εr2、μ2、μr2分别表示介质球的电导率、介电常数、相对介电常数、磁导率、相对磁导率；r,t为反射波和透射波的反射和透射系数；η1、k1为媒质1的本征阻抗、波数，η2、k2为介质球的本征阻抗、波数；0表示自由空间，1表示媒质1空间，2表示各向同性介质球；本征阻抗与波数表达是为：

根据边界条件建立方程组，在z＝-d的分界平面上，应有

得出反射系数r和透射系数t的表达式分别为：

将(9)(10)带入式(3)(4)得到在z＝-d处的反射波和透射波的电场值

本实施例步骤2中介质球界面的反射波、透射波的电场和磁场的具体计算如下：

各向同性介质球的入射波根据本征矢量和平面波因子乘积的解析展开式，将入射波的电磁场用各向同性球矢量波函数展开的具体表达式为：

各向同性介质球的入射波沿着正z方向入射，其电场的极化方向是正x方向，其电场的幅度为t，其中

经各向同性介质球面散射后，利用均匀平面波在无界均匀媒质中的传播特性，反射波和透射波的电场与磁场分别为：

其中，表示n从1到∞，m从-n到+n的求和，表示反射电场与透射电场的反射系数和透射系数；表示球矢量波函数，l＝1，l＝3分别表示矢量波函数由第一类、第三类球贝塞尔函数构成，表示球坐标系中的一个矢量；通过边界条件，在r＝a时，电磁场在切向边界上连续，即有：

联立方程组得反射电场和透射电场的展开系数

其中jn(kr)为球贝塞尔函数，为第一类球汉克尔函数,p的具体表达式：

将(17a-d)带入(15a-d)得出各向同性介质球界面的反射波与透射波的电场和磁场的具体值

本实施例步骤3中全局坐标系下的透射场与镜像坐标系下的反射场的具体值计算如下：

全局坐标系：球坐标系中的三个坐标变量r,θ,φ；

局部坐标系：球坐标系中的三个坐标变量rimg,θimg,φimg；

参见图2，在球坐标系中，过空间任意一点的三个相互正交的坐标单位矢量分别是r,θ,增加的方向，且遵循右手螺旋法则，它们与之间的变换关系为：

经各向同性介质球散射，其反射波作为入射波入射到分界平面上，利用介质平面的镜像原理其反射波和透射波的电场和磁场可表示为：

其中反射波的电场和磁场用局部坐标系表示，透射波的电场和磁场用全局坐标系表示，根据边界条件及p点θ＝π，θimg＝0，r＝rimg＝d，可以计算得到展开系数

其中mnm(θ,φ)，nnm(θ,φ)的具体表达式为：

联立(20a-d)(21)得出全局坐标系下的透射场与镜像坐标系下的反射场的具体值。

本实施例步骤4中利用球矢量波函数的叠加原理来求得在全局坐标系下对应反射场然后采用递归t矩阵算法和传统理论方法得出矩阵变换系数，并建立数据库，提高计算效率与精度，最后得出的具体值，计算如下：

利用球矢量波函数的叠加原理将局部坐标系表示的反射电场转换成全局坐标系表示的反射电场全局坐标系下的反射电场的具体表达式为：

参见图3，三维矢量波函数加法定理下的之间的关系可以得到具体的表达式：

式中表示v从1到∞，u是从-v到+v的累加。是第一类三维矢量波函数；称为坐标系变换系数；

联立(20a)，(23)，(24)，(25)式，得出矩阵形式表示如下：

使用传统的方法计算矩阵变换系数首先由平面波的展开公式：

其中θk,φk表示k的指向方位角，θ,φ表示r的指向方位角，jl表示l阶球bessel，ylm表示归一化的球贝塞尔函数

由r＝r'+r”得

e^ik·r＝e^ik·r'·e^ik·r”(29)

联立(27)(28)(29)式根据ylm(θk,φk)的正交归一性得

其中是wigner3-j符号,

由此得

其中zl(x)等于四类贝塞尔函数，

由以上推到可以看出，由此方法计算坐标系变换系数的运算量是相当大的，每个βl'm',lm,αl'm',lm需要对l”进行求和，而且wigner3-j符号依赖与多种因素，使计算量增加。

采用递归t矩阵算法进行求解坐标系变换系数根据经典的mie散射理论可以得出单个散射球的t矩阵表示为

根据矩阵变换系数建立并矢坐标变换矩阵

其中的下角标ji表示坐标系i到坐标系j的变换；的上角标(j),(h)表示坐标变换系数中基函数的选取为球的bessel函数或第一类球hankel函数；递推方法的初值：

故球谐的加法公式可以写成：

根据球面波函数的性质

其中是连带的勒让德函数，联立(38)(39)(40)(41)得

再利用连带的勒让德函数、贝塞尔函数、球矢量波函数得：

其中为系数

参见图4，根据(43)式发现，当给定一定的n，如果知道nn的一系列的转移矩阵元素值，那么就可以知道一系列的(n+1,n+1)的转移矩阵元素，所以(43)式提供的是一条沿m,n的对角线的递推路径。

根据(44)式发现当知道(n+1,m),(n,m)的波函数的加法定理系数后，通过递推关系可以求得((n+1,m)系数，在这个式子里m是不变的，它的递归路线是向上的。

联立(37)(42)(43)(44)得

其中

在计算矩阵变换系数时，建立阶乘数据库，把递推的数据存放在数据库中，根据需求调用数据库中的数值，这样计算方式比直接计算方法大大节省时间，提高了计算效率。最后联立(23)(26)(43)(48)的出全局坐标系下的反射电场的值。

通过改变频率，媒质1的电导率、相对介电常数，埋藏深度，将在p点的电场的仿真值以及程序所得到的电场值进行对比来验证该方法的正确性，进而得出参数对电场值的影响。

通过三维电磁仿真软件feko，设置近场点，得出的值即为所求的电场值，然后通过fortran编程语言根据求解的电场的具体值，再通过迭代循环，建立方程组，最后得到在p点处的电场值即程序所得到的电场值。

由图5可以看出频率在70khz-700khz的情况下仿真所得到的电场值和根据公式推导所得到的电场值匹配度良好，随着频率的增加电场值逐渐增加。

由图6可以得出结论：当土壤介质的电导率的减小时，电场值变大。

图7对比可以得出结论：当土壤介质的相对介电常数增加时，电场值基本上不变。

由图8可以得出结论：当埋藏深度增加时，电场值变小。

由图9可以得出结论：地下有掩埋物时，会使电场值变小。

以上对本发明的优选实施例及原理进行了详细说明，对本领域的普通技术人员而言，依据本发明提供的思想，在具体实施方式上会有改变之处，而这些改变也应视为本发明的保护范围。