CN113673070B - 一种任意姿态电性天线辐射电磁场的快速计算方法 - Google Patents

Abstract

本发明涉及一种任意姿态电性天线辐射电磁场的快速计算方法，包括：获取包括天线长度、天线电流和天线姿态参数在内的天线参数，并设置高斯‑切比雪夫积分参数；基于天线姿态参数对天线参数进行投影变换，得到天线长度在x、y和z三个方向的天线长度值和天线电偶极矩三分量；根据设置的高斯‑切比雪夫积分参数和天线长度值计算出各积分节点坐标和对应的积分系数，并利用各积分节点坐标和天线电偶极距三分量得到各积分节点位置处电偶极子天线三分量的频域电磁响应；利用各积分节点对应的积分系数和各积分节点位置处电偶极子天线三分量的频域电磁响应得到任意姿态电性天线的辐射电磁场。本发明实现快速计算辐射电磁场的同时提升了计算的精度和效率。

Description

一种任意姿态电性天线辐射电磁场的快速计算方法

技术领域

本发明涉及电磁探测技术领域，尤其涉及一种任意姿态电性天线辐射电磁场的快速计算方法。

背景技术

目前，海洋目标电磁探测技术领域中如海底石油、水合物等探测采用海洋可控源电磁方法和海水中目标探测方法。其中，海洋可控源电磁方法(MCSEM)通常采用几百米长的水平电性天线在海水中(距海底几十米的位置)辐射峰值电流几百安培至千安培、基频n×10^-1Hz至n×10Hz范围内的矩形波电流，通过布置在海底或者拖曳在距水平电性天线固定偏移距的电场或磁场传感器观测电场/磁场响应信号，然后通过采用相适应的数据处理手段对电磁信号进行处理，运用预先建立的层状海洋模型正反演算法获得对实测电磁信号的定量反演解释，从而得到海洋中目标电阻率信息。

MCSEM中现有的正演仿真计算方法主要由两种，一种是将有限长电性天线看作电偶极子天线，然后利用层状海洋模型下电偶极子天线辐射的电磁场模拟解释海洋目标及海底电性参数；另一种方法是将有限长度电性天线进行均匀分割，然后把分割后的每段等效为电偶极子，然后对每个等效电偶极子天线的电磁响应求和得到有限长度电性天线的辐射电磁场。第一种方法将有限长度电性天线看作电偶极子天线，MCSEM层状海洋模型的正演计算即是计算海水中电偶极子天线辐射的电磁场。这种方法认为，当海洋目标(海底或海水目标)和观测传感器距离电性天线很远(相比于电性天线长度)时电性天线可以等效为电偶极子天线。然而，实际MCSEM中的电性天线长几百米(一般100米至300米)，距离海底几十米距离(一般50米左右)，而且航行作业时，电性天线从靠近海底的传感器阵列到远离传感器阵列的过程中，无法始终保证电性天线距离海底和传感器的距离远大于电性天线的长度(如5倍电性天线长度)，从而导致该方法在计算小偏移距的电磁场时产生较大的计算误差。第二种方法均匀分割有限长度电性天线，然后把分割后的每段等效为电偶极子，并对每个等效电偶极子天线的电磁响应求和得到有限长度电性天线的辐射电磁场。应用中该方法可以通过迭代实现也可以直接设定分割数求解辐射电磁场，但是无论是直接设定分割数进行计算还是通过迭代法进行辐射电磁场的计算，总之均匀分割法在保证计算精度的同时会极大的降低计算效率，不利于工程应用。

发明内容

鉴于上述的分析，本发明旨在提供一种任意姿态电性天线辐射电磁场的快速计算方法，以解决现有电磁探测技术中采用均匀分割法计算有限长度电性天线辐射电磁场时在保证计算精度的同时会极大的降低计算效率，致使不利于工程应用的问题。

本发明的目的主要是通过以下技术方案实现的：

本发明提供了一种任意姿态电性天线辐射电磁场的快速计算方法，包括以下步骤：S1，获取包括天线长度、天线电流和天线姿态参数在内的天线参数，并设置高斯-切比雪夫积分参数；S2，基于所述天线姿态参数对所述天线参数进行投影变换，得到天线长度在x、y和z三个方向的天线长度值和天线电流在x、y和z三个方向的天线电偶极矩三分量；其中，x、y和z分别为天线的X轴、Y轴和Z轴；S3，根据设置的所述高斯-切比雪夫积分参数和所述天线长度值计算出各积分节点坐标和对应的积分系数，并利用所述各积分节点坐标和所述天线电偶极距三分量得到各积分节点位置处电偶极子天线三分量的频域电磁响应；S4，利用所述各积分节点对应的积分系数和所述各积分节点位置处电偶极子天线三分量的频域电磁响应进行高斯-切比雪夫积分计算得到任意姿态电性天线的辐射电磁场。

进一步地，所述S4中利用所述各积分节点对应的积分系数和所述各积分节点位置处电偶极子天线三分量的频域电磁响应进行高斯-切比雪夫积分计算的构造公式为：

其中，L为天线长度，A_k为各积分节点对应的积分系数，p_k为天线各积分节点位置处电偶极矩，θ和均为天线姿态参数，为天线各积分节点坐标的三分量和与天线各积分节点的电偶极矩三分量和一一对应，为观测点相对于天线各积分节点的位置矢量，x_k为高斯-切比雪夫积分节点，k＝0,1,2,…,n。

进一步地，所述S3具体包括以下步骤：

S31，根据n+1次切比雪夫多项式T_n+1(x)求解n+1个零点，得到n+1个积分节点x_k，其中，k＝0,1,2,…,n，n为设置的所述高斯-切比雪夫积分参数；

S32，根据各积分节点x_k对应的拉格朗日插值基函数L_k(x)计算各积分节点对应的积分系数A_k；

S33，根据所述天线长度值对应的坐标计算出天线各积分节点坐标的三分量和

S34，采用所述天线参数对天线各积分节点位置处电偶极矩进行投影变换，得到与所述天线各积分节点坐标的三分量和一一对应的天线电偶极拒的三分量和

S35，利用所述天线各积分节点坐标的三分量和及一一对应的所述天线电偶极矩的三分量和并通过快速Hankel数值滤波得到各积分节点位置处电偶极子天线三分量的频域电磁响应。

进一步地，所述S31中T_n+1(x)求解n+1个零点得到积分节点x_k的计算公式如下：

其中，T_n+1(x)的n+1个零点为x_k，k＝0,1,2,…,n。

进一步地，所述高斯-切比雪夫积分节点x_k计算公式如下：

其中，n为设置的所述高斯-切比雪夫积分参数。

进一步地，所述各积分节点x_k对应的积分系数A_k计算公式如下：

其中，x_k对应的拉格朗日插值基函数L_k(x)公式如下：

进一步地，所述S33中根据电性天线长度值对应的坐标计算天线各积分节点坐标的三分量和包括：

当计算天线的X轴分量时，b＝L_x/2，a＝-L_x/2，计算得到天线X轴分量各分点坐标

当计算天线的Y轴分量时，b＝L_y/2，a＝-L_y/2，计算得到天线Y轴分量各分点坐标

当计算天线的Z轴分量时，b＝L_z/2，a＝-L_z/2，计算得到天线Z轴分量各分点坐标

其中，L_x、L_y和L_z为天线长度在x、y和z三个方向的天线长度值，a和b分别为所述天线长度值对应的坐标。

进一步地，所述天线长度在x、y和z三个方向的天线长度值L_x、L_y和L_z的公式如下：

其中，L为天线长度，θ和均为天线姿态参数，θ为有限长电性天线在XOY平面内的投影与x轴夹角，为有限长电性天线与XOY平面的夹角。

进一步地，所述S34中采用所述天线参数对天线各积分节点位置处电偶极矩进行投影变换，得到与所述天线各积分节点坐标的三分量和一一对应的天线电偶极拒的三分量和的公式如下：

其中，I为天线电流，p_k＝Idl_k为天线各积分节点位置处电偶极矩，θ为有限长电性天线在XOY平面内的投影与x轴夹角，为有限长电性天线与XOY平面的夹角。

进一步地，所述任意姿态电性天线的辐射电磁场为频域电磁场，通过GS变换转换为任意姿态电性天线的辐射时域电磁场；

所述GS变换的计算公式如下：

其中，函数F(*)表示频域函数，函数f(t)表示时域函数，n_f为变换时采用的频点数目，f为相应的频点值，K_m是GS变换系数。

本发明技术方案的有益效果如下：本发明公开了一种任意姿态电性天线辐射电磁场的快速计算方法，考虑到电性天线的长度影响，为提高短偏移距辐射电磁场的计算精度，将任意姿态有限长度电性天线进行少点数剖分，以保证各分点处电偶极子天线进行辐射电磁场的精确仿真计算，并解决了现有电性天线均匀分割技术在保证计算精度时需要精细剖分而严重影响计算效率的问题。

本发明的其他特征和优点将在随后的说明书中阐述，并且，部分的从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在所写的说明书、权利要求书、以及附图中所特别指出的结构来实现和获得。

附图说明

附图仅用于示出具体实施例的目的，而并不认为是对本发明的限制，在整个附图中，相同的参考符号表示相同的部件。

图1为本发明实施例的一种任意姿态电性天线辐射电磁场的快速计算方法流程图；

图2为本发明实施例的典型层状海洋模型下任意姿态电性天线及高斯-切比雪夫分点示意图；

图3为本发明实施例的典型层状海洋模型下任意姿态电性天线辐射电磁场的非均匀稀疏分割点的高斯-切比雪夫积分计算方法流程图；

图4为本发明实施例的任意姿态天线的辐射时域电磁场的计算方法流程图；

图5为本发明实施例的海底沿x轴在0-1km范围内的磁感应强度总场与各磁场分量的幅度分布曲线图；

图6为本发明实施例的海底沿x轴在0-1km范围内的电场强度总场与各电场分量的幅度分布曲线图；

图7为本发明实施例的海底沿x轴在0-1km范围内的磁感应强度总场与各磁场分量的相位分布曲线图；

图8为本发明实施例的海底沿x轴在0-1km范围内的电场强度总场与各电场分量的相位分布曲线图；

图9为本发明实施例的海底沿y轴在0-1km范围内的磁感应强度总场与各磁场分量的幅度分布曲线图；

图10为本发明实施例的海底沿y轴在0-1km范围内的电场强度总场与各电场分量的幅度分布曲线图；

图11为本发明实施例的海底沿y轴在0-1km范围内的磁感应强度总场与各磁场分量的相位分布曲线图；

图12为本发明实施例的海底沿y轴在0-1km范围内的电场强度总场与各电场分量的相位分布曲线图；

图13为本发明实施例的海底(0,500m)观测位置的磁感应总场强度及各分量的负阶跃响应曲线图；

图14为本发明实施例的海底(0,500m)观测位置的电场总场强度及各分量的负阶跃响应曲线图。

具体实施方式

下面结合附图来具体描述本发明的优选实施例，其中，附图构成本申请一部分，并与本发明的实施例一起用于阐释本发明的原理，并非用于限定本发明的范围。

本发明的技术构思：针对现有方案无法兼顾计算精度和计算效率的不足，提出了一种典型层状海洋模型中任意姿态有限长度电性天线辐射电磁场的精确快速计算方法，采用切比雪夫多项式的零点作为分点，对原有的均匀分割方案进行了以切比雪夫多项式零点为分点的不均匀分割方案的改进，从而实现了有限长电性天线的精确快速计算，兼顾了计算的精度和效率。现有主要技术方案及其缺点或不足：

现有技术方案一：将有限长度电性天线看作电偶极子天线，MCSEM层状海洋模型的正演计算即是计算海水中电偶极子天线辐射的电磁场。这种方法认为，当海洋目标(海底或海水目标)和观测传感器距离电性天线很远(相比于电性天线长度)时电性天线可以等效为电偶极子天线。然而，实际MCSEM中的电性天线长几百米(一般100米至300米)，距离海底几十米距离(一般50米左右)，而且航行作业时，电性天线从靠近海底的传感器阵列到远离传感器阵列的过程中，无法始终保证电性天线距离海底和传感器的距离远大于电性天线的长度(如5倍电性天线长度)，从而导致该方法在计算小偏移距的电磁场时产生较大的计算误差。

现有技术方案二：均匀分割有限长电性天线，然后把分割后的每段等效为电偶极子，并对每个等效电偶极子天线的电磁响应求和得到有限长度电性天线的辐射电磁场。应用中该方法可以通过迭代实现也可以直接设定分割数求解辐射电磁场。

迭代法的操作方法为，首先给定初始分割数N₀和阈值Δ，然后按照N₀和2N₀对电性天线长度进行均匀剖分，并计算各分段电偶极子天线的辐射电磁场E_xi、E_yi、E_zi、B_xi、B_yi、B_zi和E_xj、E_yj、E_zj、B_xj、B_yj、B_zj，其中i＝1…N₀和j＝1…2N₀。然后分别求N₀和2N₀分段数对应各分段电偶极子天线的辐射电磁场分量并计算相对误差，若满足收敛条件，则选定均匀分割点数为N₀，否则加倍更新均匀分割数后再次比较两种分割数下辐射电磁场计算结果的误差，直到满足收敛条件并输出当前分割点数N₀和辐射电磁场结果。无论是直接设定分割数进行计算还是通过迭代法进行辐射电磁场的计算，总之均匀分割法在保证计算精度的同时会极大的降低计算效率，不利于工程应用中数据资料的快速解释。

本发明的一个具体实施例，如图1所示，公开了一种任意姿态电性天线辐射电磁场的快速计算方法，包括以下步骤：

S1，获取包括天线长度、天线电流和天线姿态参数在内的天线参数，并设置高斯-切比雪夫积分参数；

S2，基于所述天线姿态参数对所述天线参数进行投影变换，得到天线长度在x、y和z三个方向的天线长度值和天线电流在x、y和z三个方向的天线电偶极矩三分量；其中，x、y和z分别为天线的X轴、Y轴和Z轴；

S3，根据设置的所述高斯-切比雪夫积分参数和所述天线长度值计算出各积分节点坐标和对应的积分系数，并利用所述各积分节点坐标和所述天线电偶极距三分量得到各积分节点位置处电偶极子天线三分量的频域电磁响应；

S4，利用所述各积分节点对应的积分系数和所述各积分节点位置处电偶极子天线三分量的频域电磁响应进行高斯-切比雪夫积分计算得到任意姿态电性天线的辐射电磁场。

与现有均匀分割法计算任意姿态电性天线辐射电磁场相比，本发明实现任意姿态电性天线辐射电磁场快速计算的同时提升了计算精度和计算效率。

具体地，典型层状海洋模型如图2所示，分为三层，即空气层、海水层和海底层。空气-海水-海底的层状模型参数：海水及海底各层厚度d_i,i＝1,2,...,n,n≥2，空气(i＝0)、海水(i＝1)及海底(i＞1)各层介质电导率σ_i、介电参数ε_i和磁导率μ_i＝μ_rμ₀，其中i＝0,1,2,...,n,n≥2。一般将海水及海底相对磁导率μ_r设置为1。天线参数包括天线长度L，天线峰值电流I，天线电矩P＝IL，天线姿态参数θ和及天线中线距离海水表面距离(天线深度)H。观测参数：观测点坐标(x,y,z)，相对于天线中心坐标(0,0,H)的位置矢量及相对于天线各积分节点的位置矢量其中，空气层设定为半无限大均匀空间，海底层为半无限大均匀空间或层状空间，海水层深度为d₁，设定分解面互相平行，水平方向上无限延伸，空气、海水和海底(可以为多层海底)的介电参数分别为σ_i、ε_i(i＝0,1,2,...,n，n≥2)，其中i＝0表示空气层，i＝1表示海水层。真空磁导率为μ₀，相对磁导率μ_r设置为1。AB表示电性天线首尾，天线长度为L，天线电极矩为P＝IL，天线中点位于原点O正下方，目标所在深度为H，距离海底深度为h，各层厚度为d_i。笛卡尔坐标系正z方向垂直水平面向下，xyz满足右手螺旋定则。任意姿态电性天线在XOY平面内的投影与x轴夹角为θ，天线与XOY平面的夹角为

如图3所示，本发明技术方案的具体实施步骤如下：

第一步，设置海洋模型参数、天线参数和观测参数。

海洋模型参数包括空气-海水-海底层状模型的层厚度、电导率、介电常数和磁导率，海洋模型参数用于计算电偶极子天线辐射电磁场的积分核函数。天线参数包括天线位置、天线长度、天线姿态、天线电流和天线电矩。观测参数包括观测点坐标、频率和时间等信息。

第二步，设置频-时变换参数和高斯-切比雪夫积分参数。

频-时变换参数用于实现频域解到时域解的变换计算，本发明中涉及的频-时变换采用现有频时变换成熟的GS变换技术，通过数值滤波的方式实现，实现相同的目标也可以采用逆傅里叶变换技术。高斯-切比雪夫积分参数包括分割点数(亦可称为积分阶数)和误差阈值。一般采用的分割点数在7-13之间选取，误差阈值可以在10^-8-10^-12之间选择，得到的积分结果精度相当。

第三步，采用天线姿态参数进行投影变化，得到天线长度的x、y和z三个方向的天线长度值和天线电流矢量的三分量。天线长度在x、y和z三个方向的天线长度值L_x、L_y和L_z将被用于高斯-切比雪夫积分节点坐标的计算。天线电流矢量的三分量计算源自于对天线各积分节点位置处电偶极矩进行投影变换的公式。由于dl_k是天线l_k位置处的单位长度的电偶极子天线，因此该投影变化的公式变为：

第四步，根据设置的高斯-切比雪夫积分参数、天线长度在三个方向的长度值L_x、L_y和L_z计算积分节点坐标和对应的积分系数。

需要说明的是，后续流程中将获得的积分节点坐标和对应的积分系数代入高斯-切比雪夫积分计算公式，得到各积分节点位置处电偶极子天线的辐射电磁场，并与对应积分系数相乘求和得到任意姿态有限长电性天线的辐射电磁场。计算中高斯-切比雪夫积分计算公式中p_k变为I。

第五步，分别计算天线x分量、y分量和z分量的辐射电磁场。

利用得到的积分节点位置、电流矢量的三分量I_x、I_y和I_z分别采用沿不同方向(x轴、y轴和z轴)电偶极子电磁场的计算式计算各积分节点位置处电偶极子天线x,y,z三分量产生的辐射电磁场。这一步输出各积分节点位置、三分量的电偶极子天线在不同频率时的辐射电磁场值。

电偶极子天线辐射电磁场的计算采用了快速Hankel数值滤波的方法，该方法是现在计算层状海洋模型水平或垂直电偶极子天线辐射场的常用方法，故在此不做赘述。

第六步，若需要计算任意姿态有限长电性天线的辐射时域电磁场，则利用频-时转换技术将输出的辐射电磁场的频域结果转化为时域结果，否则跳转至第七步。

第七步，将计算得到的各积分节点位置处电偶极子天线的辐射电磁场的结果进行高斯-切比雪夫积分计算，获得最终任意姿态有限长电性天线的辐射电磁场。

本发明的一个具体实施例，所述S4中利用所述各积分节点对应的积分系数和所述各积分节点位置处电偶极子天线三分量的频域电磁响应进行高斯-切比雪夫积分计算的构造公式为：

其中，L为天线长度，A_k为各积分节点对应的积分系数，p_k为天线各积分节点位置处电偶极矩，θ和均为天线姿态参数，为天线各积分节点坐标的三分量和与天线各积分节点的电偶极矩三分量和一一对应，为观测点相对于天线各积分节点的位置矢量，x_k为高斯-切比雪夫积分节点，k＝0,1,2,…,n。

具体地，设电性天线l_i位置的电偶极子天线在观测位置处产生的响应(电场或磁场)表示为其中f(…)可以表示辐射电磁场E_x、E_y、E_z、B_x、B_y、B_z的计算式，θ和为天线姿态的已知参数，p_i＝Idl_i，p_i表示天线l_i位置处的电偶极矩，I为天线电流强度，dl_i表示l_i位置的电偶极子单位长度，和分别为电偶极天线和观测点的相对于坐标原点的位置矢量，观测点坐标(x,y,z)相对于天线中心坐标(0,0,H)的位置矢量因此，任意姿态天线的辐射电磁场计算式如下，

现有技术方案在求解层状海洋中任意姿态电性天线的辐射电磁场时，采用两种方法处理计算上式，即一种忽略天线的长度影响，则上式沿天线长度的积分将消失，如下式所示，

其中，表示任意姿态电性天线中心位置的位置矢量。

另一种方法对天线长度的均匀精细分割则是将上面积分式离散化为众多分点对应的电偶极子天线响应之和的形式，如下式所示，

其中，N为均匀分割的分点数，一般为保证计算精度，对于长度100米以上的天线，N取值不小100。

而本发明为解决现有方案的不足，采用高斯-切比雪夫积分方法能够实现非均匀稀疏分割和代数精度为2n+1的计算方法，兼顾计算精度和计算效率，该方法利用具备带权正交性的切比雪夫多项式的零点作为高斯积分点，采用带权值的拉格朗日插值多项式计算积分系数，并通过对积分核函数的构造实现快速精确的计算方案。观测点坐标(x,y,z)相对于天线各积分节点的位置矢量本发明的技术方案对应的计算式如下：

其中，节点数n一般取值7-13，A_k表示积分节点l_k对应的积分系数，该积分系数不依赖积分区间。

本发明的一个具体实施例，所述S3具体包括以下步骤：

S31，根据n+1次切比雪夫多项式T_n+1(x)求解n+1个零点，得到n+1个积分节点x_k，其中，k＝0,1,2,…,n，n为设置的所述高斯-切比雪夫积分参数；

S32，根据各积分节点x_k对应的拉格朗日插值基函数L_k(x)计算各积分节点对应的积分系数A_k；

S33，根据所述天线长度值对应的坐标计算出天线各积分节点坐标的三分量和

S34，采用所述天线参数对天线各积分节点位置处电偶极矩进行投影变换，得到与所述天线各积分节点坐标的三分量和一一对应的天线电偶极拒的三分量和

S35，利用所述天线各积分节点坐标的三分量和及一一对应的所述天线电偶极矩的三分量和并通过快速Hankel数值滤波得到各积分节点位置处电偶极子天线三分量的频域电磁响应。

本发明的一个具体实施例，所述S31中T_n+1(x)求解n+1个零点得到积分节点x_k的计算公式如下：

其中，T_n+1(x)的n+1个零点为x_k，k＝0,1,2,…,n。

本发明的一个具体实施例，所述高斯-切比雪夫积分节点x_k计算公式如下：

其中，n为设置的所述高斯-切比雪夫积分参数。

本发明的一个具体实施例，所述各积分节点x_k对应的积分系数A_k计算公式如下：

其中，x_k对应的拉格朗日插值基函数L_k(x)公式如下：

本发明的一个具体实施例，所述S33中根据电性天线长度值对应的坐标计算天线各积分节点坐标的三分量和包括：

当计算天线的X轴分量时，b＝L_x/2，a＝-L_x/2，计算得到天线X轴分量各分点坐标

当计算天线的Y轴分量时，b＝L_y/2，a＝-L_y/2，计算得到天线Y轴分量各分点坐标

当计算天线的Z轴分量时，b＝L_z/2，a＝-L_z/2，计算得到天线Z轴分量各分点坐标

其中，L_x、L_y和L_z为天线长度在x、y和z三个方向的天线长度值，a和b分别为所述天线长度值对应的坐标。

本发明的一个具体实施例，所述天线长度在x、y和z三个方向的天线长度值L_x、L_y和L_z的公式如下：

其中，L为天线长度，θ和均为天线姿态参数，θ为有限长电性天线在XOY平面内的投影与x轴夹角，为有限长电性天线与XOY平面的夹角。

本发明的一个具体实施例，所述S34中采用所述天线参数对天线各积分节点位置处电偶极矩进行投影变换，得到与所述天线各积分节点坐标的三分量和一一对应的天线电偶极拒的三分量和的公式如下：

其中，I为天线电流，p_k＝Idl_k为天线各积分节点位置处电偶极矩，θ为有限长电性天线在XOY平面内的投影与x轴夹角，为有限长电性天线与XOY平面的夹角。

具体地，本发明实施例中各积分节点及对应积分系数的求解，以及利用所述各积分节点坐标和所述天线电偶极矩三分量得到各积分节点位置处电偶极子天线三分量的频域电磁响应的计算，具体实施步骤如下：

第一步，求n+1次切比雪夫多项式的n+1个零点，用于后续产生积分节点。

由于切比雪夫多项式在区间[-1,1]的带权值正交性，n+1次切比雪夫多项式的n+1个零点也是高斯点，对应的高斯-切比雪夫积分方法的代数精度也是2n+1次的。切比雪夫多项式满足下面的递推关系式：

其对应的n+1个零点的计算式如下：

第二步，计算各节点对应的拉格朗日插值基函数L_k(x)，用于计算各积分节点对应的积分系数。

x_k对应的拉格朗日插值基函数L_k(x)的计算式如下，

第三步，求各积分节点对应的积分系数A_k，计算式如下，

第四步，积分区间变换，调整积分区间与有限长电性天线长度值匹配。

其中，a和b分别为电性天线长度值对应的坐标。

当计算天线的X轴分量时，b＝L_x/2，a＝-L_x/2，计算得到天线X轴分量各分点坐标

当计算天线的Y轴分量时，b＝L_y/2，a＝-L_y/2，计算得到天线Y轴分量各分点坐标

当计算天线的Z轴分量时，b＝L_z/2，a＝-L_z/2，计算得到天线Z轴分量各分点坐标

求解过程中的x_k保持不变，L_x、L_y和L_z通过投影变化得到：

第五步，通过如下高斯-切比雪夫积分计算公式得到任意姿态有限长电性天线的辐射电磁场。

在计算过程中，l_k和p_k应对应，对应关系为，与对应，与对应，与对应。其中和的计算方法如下，分别对应各积分节点位置处天线电偶极矩的三分量：

需要说明的是，高斯-切比雪夫积分计算实现中一般取dl_k＝1，因此p_k变为I。

本发明的一个具体实施例，所述任意姿态电性天线的辐射电磁场为频域电磁场，通过GS变换转换为任意姿态电性天线的辐射时域电磁场；

所述GS变换的计算公式如下：

其中，函数F(*)表示频域函数，函数f(t)表示时域函数，n_f为变换时采用的频点数目，f为相应的频点值，K_m是GS变换系数。

需要说明的是，本发明可用于计算频域电磁场也可用于计算时域电磁场，并利用频-时转换技术将输出的辐射电磁场的频域结果转化为时域结果，本发明中的频-时转换技术采用的现有成熟的GS变化技术。当计算时域电磁场时，可以直接将有限长电性天线的辐射频域电磁场通过GS变换得到有限长电性天线的辐射时域电磁场，或者先将各积分节点位置处电偶极子天线三分量的频域电磁响应通过GS变换得到各积分节点位置处电偶极子天线三分量的时域电磁响应，再通过高斯-切比雪夫积分计算得到有限长电性天线的辐射时域电磁场，如图4所示。

本发明实施例通过非均匀稀疏分割积分技术、投影变换技术、数值滤波和时频转换技术实现的，其中非均匀稀疏分割积分技术通过采用高斯分点作为非均匀稀疏分割点，并采用分点计算时对应的高斯-切比雪夫多项式积分进行快速精确积分实现非均匀稀疏分割积分技术，保证计算精度的同时极大提升计算效率；通过投影变换技术实现任意姿态电性天线向XY、XZ和YZ平面的正交投影，转换为计算水平及垂直电性天线的辐射电磁场，实现对任意姿态电性天线电磁场的计算，通过采用数值滤波和时频转换技术实现天线辐射电磁场的时频域快速计算。也就是说，本发明实施例采用高斯-切比雪夫积分方法进行典型层状海洋模型任意姿态有限长电性天线电磁场的计算方法，是以n+1次切比雪夫多项式的零点作为积分节点构建代数精度为2n+1次的积分求解方法，实现了非均匀稀疏分割积分方法，并在保证计算精度的同时能够极大地提升计算效率。

与现有技术方案一相比，本发明考虑了天线的有限长度的影响，设计的具备2n+1代数精度的高斯-切比雪夫积分方法，即使在非均匀稀疏分割的情况下仍能保证计算精度，特别是保证了MCSEM短收发偏移距时的辐射电磁场的计算精度，提升了数据资料解释准确性。

与现有技术方案二相比，本发明设计的技术方案无需进行均匀精细分割即可保证有限长度电性天线辐射电磁场的计算精度，通过以n+1次切比雪夫多项式的零点作为高斯积分点实现非均匀稀疏分割，同时由于积分点和对应系数与积分区间无关，而仅与阶数有关，因此，应用中可以先将常用阶数的积分节点和积分系数存储，在后续MCSEM原理仿真与数据解释时直接调用即可，进一步提升了计算效率。

为说明本发明技术方案相比现有技术方案的有益效果，给出下列计算案例对比，采用典型的层状海洋模型，详细的仿真参数见表1。

表1高斯-切比雪夫方法仿真参数

图5-图12给出了海底沿x方向和y方向0-1km范围内磁感应总场、各磁场分量、电场总场和各电场分量的幅度和相位分布曲线，图13和图14给出了海底(0,500m)位置的磁感应总场、各磁场分量、电场总场和各电场分量的负阶跃响应曲线。其中，实线为电偶极子天线近似结果，划线为均匀精细分割计算结果，点线为本发明的高斯-切比雪夫积分方法计算结果。

计算效率：针对单频点，200个观测点六个分量的仿真计算耗用时间见表2。所用仿真计算机配置为Win7系统，六核Intel i5-8400，主频2.80GHz，8GB RAM，仿真软件为MatlabR2018a-64位版本。计算耗时表明本发明技术方案的效率提升。

表2计算耗时统计，单位s

注:仿真计算平台配置计算机配置为Win7系统，六核Intel i5-8400，主频2.80GHz，8GB RAM，仿真软件为MatlabR2018a-64位版本

计算精度：以200点均匀精细分割计算结果为参考，频域计算结果见图5-图12和时域计算结果见图13和图14，说明了本发明技术方案的计算精度提升。

综上所述，本发明公开了一种任意姿态电性天线辐射电磁场的快速计算方法，包括以下步骤：S1，获取包括天线长度、天线电流和天线姿态参数在内的天线参数，并设置高斯-切比雪夫积分参数；S2，基于所述天线姿态参数对所述天线参数进行投影变换，得到天线长度在x、y和z三个方向的天线长度值和天线电流在x、y和z三个方向的天线电偶极矩三分量；其中，x、y和z分别为天线的X轴、Y轴和Z轴；S3，根据设置的所述高斯-切比雪夫积分参数和所述天线长度值计算出各积分节点坐标和对应的积分系数，并利用所述各积分节点坐标和所述天线电偶极距三分量得到各积分节点位置处电偶极子天线三分量的频域电磁响应；S4，利用所述各积分节点对应的积分系数和所述各积分节点位置处电偶极子天线三分量的频域电磁响应进行高斯-切比雪夫积分计算得到任意姿态电性天线的辐射电磁场。本发明采用切比雪夫多项式的零点作为分点，对原有的均匀分割方案进行了以切比雪夫多项式零点为分点的不均匀分割方案的改进，从而实现了任意姿态电性天线的精确快速计算，兼顾了计算的精度和效率。

本领域技术人员可以理解，实现上述实施例中方法的全部或部分流程，可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于计算机可读存储介质中。其中，所述计算机可读存储介质为磁盘、光盘、只读存储记忆体或随机存储记忆体等。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。

Claims (2)

1.一种任意姿态电性天线辐射电磁场的快速计算方法，其特征在于，包括以下步骤：

S1，获取包括天线长度、天线电流和天线姿态参数在内的天线参数，并设置高斯-切比雪夫积分参数；

S2，基于所述天线姿态参数对所述天线参数进行投影变换，得到天线长度在x、y和z三个方向的天线长度值和天线电流在x、y和z三个方向的天线电偶极矩三分量；其中，x、y和z分别为天线的X轴、Y轴和Z轴；

S3，根据设置的所述高斯-切比雪夫积分参数和所述天线长度值计算出各积分节点坐标和对应的积分系数，并利用所述各积分节点坐标和所述天线电偶极距三分量得到各积分节点位置处电偶极子天线三分量的频域电磁响应，包括：

S31，根据n+1次切比雪夫多项式T_n+1(x)求解n+1个零点，得到n+1个积分节点x_k，其计算公式为：

其中，T_n+1(x)的n+1个零点为x_k，k＝0,1,2,…,n，n为设置的所述高斯-切比雪夫积分参数；

S32，根据各积分节点x_k对应的拉格朗日插值基函数L_k(x)计算各积分节点对应的积分系数A_k，其计算公式如下：

S33，根据所述天线长度值对应的坐标计算出天线各积分节点坐标的三分量和包括：

当计算天线的X轴分量时，b＝L_x/2，a＝-L_x/2，计算得到天线X轴分量各分点坐标当计算天线的Y轴分量时，b＝L_y/2，a＝-L_y/2，计算得到天线Y轴分量各分点坐标当计算天线的Z轴分量时，b＝L_z/2，a＝-L_z/2，计算得到天线Z轴分量各分点坐标

其中，L_x、L_y和L_z为天线长度在x、y和z三个方向的天线长度值，a和b分别为所述天线长度值对应的坐标；

其中，所述天线长度在x、y和z三个方向的天线长度值L_x、L_y和L_z的公式如下：

其中，L为天线长度，θ和均为天线姿态参数，θ为有限长电性天线在XOY平面内的投影与x轴夹角，为有限长电性天线与XOY平面的夹角；

S34，采用所述天线参数对天线各积分节点位置处电偶极矩进行投影变换，得到与所述天线各积分节点坐标的三分量和一一对应的天线电偶极拒的三分量和其计算公式如下：

其中，I为天线电流，p_k＝Idl_k为天线各积分节点位置处电偶极矩，θ为有限长电性天线在XOY平面内的投影与x轴夹角，为有限长电性天线与XOY平面的夹角；

S35，利用所述天线各积分节点坐标的三分量和及一一对应的所述天线电偶极矩的三分量和并通过快速Hankel数值滤波得到各积分节点位置处电偶极子天线三分量的频域电磁响应；

S4，利用所述各积分节点对应的积分系数和所述各积分节点位置处电偶极子天线三分量的频域电磁响应进行高斯-切比雪夫积分计算得到任意姿态电性天线的辐射电磁场；

所述利用所述各积分节点对应的积分系数和所述各积分节点位置处电偶极子天线三分量的频域电磁响应进行高斯-切比雪夫积分计算的构造公式为：

其中，L为天线长度，A_k为各积分节点对应的积分系数，p_k为天线各积分节点位置处电偶极矩，θ和均为天线姿态参数，为天线各积分节点坐标的三分量和与天线各积分节点的电偶极矩三分量和一一对应，为观测点相对于天线各积分节点的位置矢量，x_k为高斯-切比雪夫积分节点，k＝0,1,2,…,n。

2.根据权利要求1所述的方法，其特征还在于，所述任意姿态电性天线的辐射电磁场为频域电磁场，通过GS变换转换为任意姿态电性天线的辐射时域电磁场；

所述GS变换的计算公式如下：

其中，函数F(*)表示频域函数，函数f(t)表示时域函数，n_f为变换时采用的频点数目，f为相应的频点值，K_m是GS变换系数。