CN110765687B - 基于时域有限差分方法的涡旋波束源合成方法 - Google Patents

Abstract

基于时域有限差分方法的涡旋波束源合成方法，本发明涉及涡旋波束源合成方法。本发明目的是为了解决现有方法研究粒子散射带来硬件上、频段上的局限，以及计算效率低、计算时间长、计算结果准确率低的问题。基于时域有限差分方法的涡旋波束源合成方法具体过程为：步骤一、获取频域的简谐平面波；步骤二、根据频域的简谐平面波，得出时域高斯脉冲；步骤三、基于时域高斯脉冲函数，得到涡旋电磁场在时域所需要的脉冲。本发明属于微波工程技术领域、计算电磁学算法领域。

Description

基于时域有限差分方法的涡旋波束源合成方法

技术领域

本发明涉及涡旋波束源合成方法，属于微波工程技术领域、计算电磁学算法领域。

背景技术

随着科技的发展，现代研究表明电磁波不仅能够携带能量，也能携带动量，而除了能携带基于坡印廷矢量的线性动量，也能携带角动量，包括自旋角动量(Spin AngularMomentum，简称SAM)和轨道角动量(Orbital Anglar Momentum,简称OAM)。目前对于携带OAM的电磁波合成方式主要有天线、阵列天线、超表面等，但是当利用携带OAM的电磁波研究粒子散射特性时，这几种方法表显得力不从心。首先，这几种方式均为频域求解方法，当想观察诸多个频点处粒子的散射特性时会大大地增加计算量，因此受限于其较低的计算效率；其次，无法捕捉到问题空间中丰富的时域信息，进而在微观角度上无法对一些现象做出解释；最后，若所需频段提高，则所需硬件尺寸缩小，现有的硬件加工技术无法达到所需要的精度，因此也受限于目前的硬件加工技术的精度。

对于数值方法，目前有基于积分方程的矩量法(Method of Moments，MoM)、基于变分原理的有限元法(Finite Element Method)，但由于人们解决问题逐渐复杂，复杂程度或体现在材料上，例如德拜模型；或体现在结构上，诸如孔，缝，槽；或体现在体积上，诸如飞机模型等；或体现在不同领域上，如人体皮肤的生物学领域上等，致使数值频域求解技术也受限于其较低的计算效率，从而不满足科研需求。对于这些问题，频域求解往往受限于其较低的计算效率，而显得力不从心。虽然数值计算方法一直作为各类工程试验的理论依据，但还是由于计算效率低以及物理现象的直观度不高，从而极大限制了数值技术的应用前景。计算机硬件的发展推动了解决问题由频域到时域的步伐，时域有限差分方法(FiniteDifference Time Domain Method，FDTD)的出现，使人们逐渐具备了在时域对具有宽频带特性的瞬变电磁场问题的计算分析条件，从而可以对物理量和物理现象更直观的理解，时域有限差分法技术的一个突出优点是可以给出问题空间丰富的时域信息，只需要经过简单的时频变换，便可得到宽频带范围的频域信息，相对于频域方法显著节约了计算量。

发明内容

本发明目的是为了解决现有方法研究粒子散射带来硬件上、频段上的局限，以及计算效率低、计算时间长、计算结果准确率低的问题，而提出基于时域有限差分方法的涡旋波束源的合成方法。

基于时域有限差分方法的涡旋波束源合成方法具体过程为：

步骤一、获取频域的简谐平面波；

步骤二、根据频域的简谐平面波，得出时域高斯脉冲；

步骤三、基于时域高斯脉冲函数，得到涡旋电磁场在时域所需要的脉冲。

发明效果

由于在粒子散射领域没有人利用时域有限差分方法合成OAM源，因此本发明填补了这一空白，为研究粒子对携带OAM电磁场的散射特性提供了解析解，在散射领域增添重要的一笔。

对于粒子散射，前人们利用平面波这种远场源对目标进行研究。在频域，平面波具有等幅等相的性质，在FDTD仿真空间中，平面波的时域波形往往采用高斯脉冲，由于对其进行频谱变换后，能够得到宽频带的频域信息，因此被大部分研究者所采用的，解决了现有方法研究粒子散射带来硬件上、频段上的局限。

本发明对于脉冲形式依然基于高斯脉冲而提出的一种OAM脉冲，在时域依然是一个可以观察到的高斯脉冲，但是在频域由此脉冲产生的电磁场的等相位面是螺旋形状，产生的电磁场幅度为中心带有奇点的圆形分布的场，与平面波的等相位面为平面大大不同。

首先利用平面波验证本算法的正确性，其次将源换成OAM源，激励源波形依然是高斯脉冲波形，在不放入介质球的情况下，观察涡旋电磁场在时域的能量分布，此波形可以在自由空间中在微波段合成频率为10GHz的涡旋波束，模式数为3；其次在仿真空间中放入介质球，研究介质球对OAM波束的散射特性，受FDTD数值稳定度影响，时域信号在频域最大频率可达到15GHz，而采样频率为60GHz，符合奈奎斯特采样定律，研究介质球对其的散射特性，并给出RCS的一些结果，解决了现有方法研究粒子散射带来硬件上、频段上的局限。本发明所提出的脉冲形式可以广泛用于微波、毫米波、太赫兹波段的雷达，本发明中对粒子的散射特性研究成果可以为微波、毫米波、太赫兹波段雷达成像提供有力的参考，解决现有方法研究粒子散射带来计算效率低、计算时间长、计算结果准确率低的问题。

附图说明

图1为在FDTD建立的仿真空间中的介质球轮廓图；

图2为介质球差分近似的细节图；

图3为太赫兹频段入射高斯脉冲采样图，Incident electric field为入射电场强度，Volt为电压单位(伏)，meter为距离单位(米)，Volt/meter为电场强度单位(伏每米)，E_θ,inc为入射电场强度沿θ方向的分量，E_φ,inc为入射电场强度沿

方向的分量，time为时间，ns为纳秒；

图4为太赫兹频段在原点处电场采样图，sampled electric field为采样电场强度；

图5为介质球在xy平面收发分离RCS极坐标下曲线图，xy plane为xy平面，θ为球坐标系下自变量(球坐标系下空间任意一点到原点的连线与z轴形成的角度)，φ为球坐标系下自变量(直角坐标系下空间任意一点到xoy平面的垂线所在的包含z轴的平面与xoz平面形成的角度在球坐标系下的表示)，RCS_θ为雷达散射截面积沿θ方向的分量，RCS_φ为雷达散射截面积沿

方向的分量，f为频率；

图6为介质球在xz平面收发分离RCS极坐标下曲线图，xz plane为xz平面；

图7为介质球在yz平面收发分离RCS极坐标下曲线图，yz plane为yz平面；

图8为介质球收发分离RCS在yz平面中的

分量与解析结果的对比图，DielectricSphere为介质球，freq代表频率，THZ为太赫兹，r为介质球的半径，ε_r为介质球相对介电常数，μ_r为介质球的相对磁导率，Ana为解析解，FDTD为本仿真的结果，θ⁰为本仿真中的球坐标系下自变量，φ⁰为本仿真中的球坐标系下自变量，

为本仿真中的雷达散射截面积

方向上的分量；dB为雷达散射截面积的单位；

图9为太赫兹频段在时域，介质球对平面波的散射图，时间步为800；

图10为太赫兹频段在时域，介质球对平面波的散射图，时间步为1130；

图11为太赫兹频段在时域，介质球对平面波的散射图，时间步为1880；

图12在不放介质球的情况下，时域OAM脉冲能量分布图；

图13a为不放介质球时，加入OAM源，在频域不同频点处的介质球散射能量在0GHz～7.5GHz处xy平面的截面幅度图；

图13b为不放介质球时，加入OAM源，在频域不同频点处的介质球散射能量在8GHz～15.5GHz处xy平面的截面幅度图；

图14a为不放介质球时，加入OAM源，在频域不同频点处的介质球散射能量在0GHz～7.5GHz处xy平面的截面相位图；

图14b为不放介质球时，加入OAM源，在频域不同频点处的介质球散射能量在8GHz～15.5GHz处xy平面的截面相位图；

图15为在频点10GHz处xy平面的截面幅度图；

图16为在频点10GHz处xy平面的截面相位图；

图17为微波波段时域入射OAM脉冲采样图；

图18为微波波段仿真空间内放入介质球时原点处电场采样值曲线图；

图19a为电场采样值的频谱图，magnitude为电厂采样值的幅度，frequency为本仿真的频率范围；

图19b为采样电场相位的频谱图，phase为采样电场的相位；

图20为介质球在xy平面收发分离RCS极坐标下曲线图；

图21为介质球在xz平面收发分离RCS极坐标下曲线图；

图22为介质球在yz平面收发分离RCS极坐标下曲线图；

图23为微波波段在时域介质球对OAM波的散射图，时间步为250；

图24为微波波段在时域介质球对OAM波的散射图，时间步为530；

图25为微波波段在时域介质球对OAM波的散射图，时间步为1090；

图26a为在频域不同频点处的介质球散射能量分布图0GHz～7.5GHzxy平面的截面幅度图；

图26b为在频域不同频点处的介质球散射能量分布图8GHz～15.5GHzxy平面的截面幅度图；

图27a在频域不同频点处的介质球散射相位分布图0GHz～7.5GHzxy平面的截面相位图；

图27b在频域不同频点处的介质球散射相位分布图8GHz～15.5GHzxy平面的截面相位图；

图28a为在频点10GHz频点处xy截面的幅度图；

图28b为在频点10GHz频点处xy截面的相位图；

图29a为在频点10GHz频点处yz截面的幅度图；

图29b为在频点10GHz频点处yz截面的相位图。

具体实施方式

具体实施方式一：本实施方式基于时域有限差分方法的涡旋波束源合成方法具体过程为：

步骤一、获取频域的简谐平面波；

步骤二、根据频域的简谐平面波，得出时域高斯脉冲；

步骤三、本发明首先利用这个脉冲研究了介质球对平面波的散射，通过与介质球散射的解析解进行对比，验证算法的正确性(仿真验证)；

基于时域高斯脉冲函数，得到涡旋电磁场在时域所需要的脉冲。

具体实施方式二：本实施方式与具体实施方式一不同的是，所述步骤一中获取频域的简谐平面波；具体过程为：

频域的简谐平面波的表达式为：

式中，

为简谐平面波的指数项；E_m1为频域的简谐平面波幅度值；ω为角频率；t为时间；

为频域的简谐平面波的波矢量；

为空间任意点的矢径；exp为自然对数的底数；j是虚数单位，j²＝-1。

其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是，所述步骤二中根据频域的简谐平面波，得出时域高斯脉冲；具体过程为：

由于对于时域仿真，频域的简谐平面波形式可以直接由一个宽带的时域入射脉冲来代替，而时域和频域是通过傅里叶变换来联系的；由于高斯脉冲具有丰富的频域信息，因此利用高斯脉冲作为产生频域平面波的入射脉冲，替换后形式为：

式中，E(t)为时域高斯脉冲函数，t₀为高斯脉冲最大值对应的横坐标，量纲是时间，单位是秒；ω为角频率；E_m2为时域高斯脉冲的幅度值；

利用含有能够体现高斯脉冲性质参数的表示式(2)，将式(2)转化为：

式中，τ为决定高斯脉冲在时域和频域的宽度的参数，量纲是时间，单位是秒。

式(3)为本发明方法中平面波源所应用的最终脉冲形式。

其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是，所述步骤三中本发明首先利用这个脉冲研究了介质球对平面波的散射，通过与介质球散射的解析解进行对比，验证算法的正确性(仿真验证)；

基于时域高斯脉冲函数，得到涡旋电磁场在时域所需要的脉冲；具体过程为：

平面波的表达式及其时域脉冲形式由是(1)及(3)给出，对于携带轨道角动量(OAM)的电磁波，常称为涡旋波束，其表达形式上与平面波表达式相比多了一个相位因子

因此涡旋电磁场电场表达式为：

式中，

表示涡旋电磁场为入射场，in表示入射，l为拓扑电荷数，也就是涡旋电磁场的模式数，也可以理解为涡旋电磁场的相位在一个平面内变化2π弧度的次数，l＝3表示涡旋电磁场的相位在一个平面内变化三个2π弧度；

为方位角，此方位角为涡旋电磁场的相位在与传播方向垂直平面上变化的弧度数，范围为[0,2π]；E_m3为频域涡旋电磁场的电场幅度值。

接下来介绍如何将频域的涡旋电磁场形式转化为时域脉冲形式。利用傅里叶变换的时移性质如式(5)：

式中，h(t)为时域函数；h(t-t′)为平移t′秒的时域函数；F(f)为h(t)对应的频域函数；t′为时移因子，在高斯脉冲里面相当于t₀，对于式(5)，t′泛指原时域函数中在自变量时间t基础上平移的量，只不过在高斯脉冲里面是t₀，因此t₀只是众多t′中的一种情况。f为频域的频率；

因此要想在频域产生具有涡旋相位的电磁波，必须将高斯脉冲t₀的形式改变，所以根据式(5)，带有入射方位角

模式为l的涡旋电磁场的时移因子表示为式(6)：

t′_v表示涡旋电磁场的时移因子，可以看出t′_v随着频率f变化，因此要是想获得目标频点的涡旋电磁场只需要设置频率f即可；利用式(6)中t′_v替换式(3)中的t₀得到涡旋电磁场的时域脉冲式(7)：

式中，E_inv(t)为涡旋电磁场对应的时域脉冲；E_m4为涡旋电磁场对应的时域脉冲的幅度值。

此式为在FDTD仿真空间中涡旋电磁场时域脉冲的最终形式。

其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是，所述方位角

由与平面波传播方向垂直的平面上的网格形成，表达式为：

式中，x和y分别为传播方向垂直的二维平面上，网格交点的横坐标和纵坐标。

其它步骤及参数与具体实施方式一至四之一相同。

采用以下实施例验证本发明的有益效果：

实施例一：

本实施例具体是按照以下步骤制备的：

参照图1至图11具体说明本实施方式，对于本算法验证部分，设计为太赫兹波段介质球对平面波散射特性，并将本算法结果即在5THz频点处的xz平面以及yz平面的RCS与经典Mie散射理论进行对比，验证本算法的正确性，并以此为基础将源采用OAM源进行研究。入射脉冲传播方向为z轴正方向，仿真区域边界采用CPML边界，层数为8层，用于近远场变换的空气盒子层数为10层，仿真空间网格大小为0.0015mm，介质球半径是0.021mm，目标区域网格数是64×65×65Δx³，时间步为2000，步长为2.6×10^-6ns，高斯脉冲中τ＝5×10^-5ns，首先基于平面波源对粒子的散射进行与经典的Mie散射理论进行验证。

图1中最外层边框显示整个体FDTD仿真空间，内层虚线边框到外层边框之间为CPML边界，层数是8，内部虚线边框构成的六个平面每一平面到介质球的距离即为空气盒子厚度，为10层，介质球材料ε_r＝3，μ_r＝2，图2所展示为在FDTD仿真空间中对于介质球的差分近似细节图，图3为激励源波形为入射高斯脉冲，仿真时，脉冲直接进入介质球，这样能够节省程序运行时间，图4为在原点处对电场进行采样之后的电压值波形，可以看到2×10^-3ns之后场值趋近于0V/m，这说明CPML边界吸收效果明显，证明算法符合数值稳定度条件并收敛，图5、6、7分别为xy平面、xz平面、yz平面收发分离RCS极坐标下曲线,；图8为介质球收发分离RCS在yz平面中的

分量与解析结果的对比图，图8中蓝色实线为解析解，红色虚线问本算法的仿真结果，可以看出二者几乎吻合，说明本发明算法的正确性；图9、10、11分别为时域分别在时间步800、1130、1880下的散射能量分布图。

参照图12至图16具体说明本实施方式，在不放介质球的情况下，将平面波源换为OAM源，频段换为微波波段，依据想在频域得到10GHz的涡旋电磁场，重新设计高斯脉冲的参数，因此时间步为1.67×10^-3ns，高斯脉冲参数τ＝3.34×10^-2ns，时间步数为1200便可达到收敛，由于频率降低，为能够显示整个脉冲幅度，使空间网格尺寸增加为0.001m，将仿真空间拉伸为40×40×120Δx³，图13a、13b和图14a、14b在频域各个频点上的幅度和相位，我们想要的在10GHz频点上的模式为3的OAM源，在图15和图16分别中对应10GHz的带有轨道角动量(OAM)的电磁场的幅度图和相位图，由图15和图16我们可以很直观地观察到10GHz处幅度和相位的具体细节。图15为幅度图，其中心为奇点，几乎没有能量，大体能量分布为一个圆环，而图16为相位图，能够很清晰地观察到在xy平面内相位变化了3个2π，也就是涡旋相位，依此说明OAM源的实现。

参照图17至图28具体说明本实施方式，远场源为OAM源时，在FDTD仿真空间中放入介质球，介质球材料ε_r＝3，μ_r＝2，介质球半径为20×10^-3m，依据想在频域得到10GHz的涡旋电磁场，高斯脉冲的参数与实施方式二中相同，因此时间步为1.67×10^-3ns，高斯脉冲参数τ＝3.34×10^-2ns，时间步数为1200便可达到收敛，空间网格尺寸为0.001m，图17展示为激励源时域OAM脉冲波形曲线，图18为原点处电场强度采样值曲线，图19a为电场采样值的频谱曲线，图19b为采样电场的相位的频谱曲线，图20、图21、图22分别为10GHz频点下介质球在xy平面、xz平面、yz平面收发分离RCS极坐标下曲线，图23、图24、图25分别为在时域时间步250、530、1090下由介质球引起的散射场在xz平面上的截面场，图26a、26b、图27a、27b展示了在频段0GHz～15.5GHz，每隔0.5GHz的各频点处xy截面的幅度图和相位图，将10GHz频点处传播方向上幅度图和相位图的xy截面图细节展示于图28a和图28b，由图28a可以观察到在球内存在奇点，但是受到介质球的影响，幅度环外侧与图15相比，变得不均匀，在球边界处还出现了向外辐射趋势，说明介质球作为二次源向外进行辐射能量；图28b与图16相对比，可以观察到介质球内部依然为模式为3的涡旋相位，但是球的外侧相位发生变化，这是介质球内部粒子相互作用导致的。介质球散射的yz截面的幅度图和相位图细节展示与图29a和图29b，由图29a可以观察到在yz截面介质球的散射的能量分布，可观察到在y＝0m的平面依然存在能量奇点，但是由于介质球散射作用使得介质球内部粒子相互作用，使奇点半径缩小，但是奇点不会消失，并且可以观察到介质球的前向散射场和后向散射场；对于图29b可以观察到的是yz截面的相位，可以观察到介质球的前向散射场和后向散射场以及介质球内部的相位分布，这是之前从未观察到过的现象，可以应用于毫米波成像，图像识别等领域，并为粒子散射领域填补了空白。

本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (1)

1.基于时域有限差分方法的涡旋波束源合成方法，其特征在于：所述方法具体过程为：

步骤一、获取频域的简谐平面波；

步骤二、根据频域的简谐平面波，得出时域高斯脉冲；

步骤三、基于时域高斯脉冲函数，得到涡旋电磁场在时域所需要的脉冲；

所述步骤一中获取频域的简谐平面波；具体过程为：

频域的简谐平面波的表达式为：

式中，

为简谐平面波的指数项；E_m1为频域的简谐平面波幅度值；ω为角频率；t为时间；

为频域的简谐平面波的波矢量；

为空间任意点的矢径；exp为自然对数的底数；j是虚数单位，j²＝-1；

所述步骤二中根据频域的简谐平面波，得出时域高斯脉冲；具体过程为：

式中，E(t)为时域高斯脉冲函数，t₀为高斯脉冲最大值对应的横坐标，量纲是时间，单位是秒；ω为角频率；E_m2为时域高斯脉冲的幅度值；

将式(2)转化为：

式中，τ为决定高斯脉冲在时域和频域的宽度的参数，量纲是时间，单位是秒；

所述步骤三中基于时域高斯脉冲函数，得到涡旋电磁场在时域所需要的脉冲；具体过程为：

涡旋电磁场电场表达式为：

式中，

表示涡旋电磁场为入射场，in表示入射，l为涡旋电磁场的相位在一个平面内变化2π弧度的次数；

为方位角，范围为[0,2π]；E_m3为频域涡旋电磁场的电场幅度值；

利用傅里叶变换的时移性质如式(5)：

h(t-t′)→e^-j2πft′·F(f) (5)

式中，h(t)为时域函数；h(t-t′)为平移t′秒的时域函数；F(f)为h(t)对应的频域函数；t′为时移因子，f为频域的频率；

带有入射方位角

模式为l的涡旋电磁场的时移因子表示为式(6)：

式中，t′_v表示涡旋电磁场的时移因子；

利用式(6)中t′_v替换式(3)中的t₀得到涡旋电磁场的时域脉冲式(7)：

式中，E_inv(t)为涡旋电磁场对应的时域脉冲；E_m4为涡旋电磁场对应的时域脉冲的幅度值；

所述方位角

表达式为：

式中，x和y分别为传播方向垂直的二维平面上，网格交点的横坐标和纵坐标。

Images