CN110083904B - 基于gpu加速的量子雷达散射截面计算方法 - Google Patents

Abstract

本发明公开了一种基于GPU加速的量子雷达散射截面计算方法，具体为：对电磁场进行量子化，根据二能级自发辐射原理，求得光子状态、光子波函数，并得到量子雷达散射截面计算公式；对目标建模并进行表面三角形面元离散，得到目标剖分文件以及每个三角形面元的中心点以及外法向量，并根据目标剖分文件进行八叉树分组；对目标进行遮挡以及亮暗面判断并记录亮面三角形面元对应的顶点坐标，计算目标在入射方向上的正交投影面积；利用获得的数据，计算量子雷达散射截面。本发明使用GPU加速，能够大幅度缩短计算时间，并且能够分析电大尺寸复杂模型的量子雷达目标特性。

Description

基于GPU加速的量子雷达散射截面计算方法

技术领域

本发明属于量子雷达目标探测技术领域，具体为一种基于GPU加速的量子雷达散射截面计算方法。

背景技术

按照经典理论进行雷达设计，则回波的信噪比将受到标准量子极限的限制而无法进一步提高。理论上证明，通过有效的量子操作提高雷达探测的性能，能够突破标准量子极限乃至达到海森堡极限，极大地提高雷达测距、测角分辨率及成像分辨率，具有广阔的应用前景。而量子雷达的发射信号表现为一小束乃至单个光子，其雷达-目标的相互作用可以描述成光子-原子的散射过程，因此传统电磁理论将无法解释这一过程，并且现阶段只能分析一些典型目标的量子雷达散射截面积，模型较为简单，对于分析一些更贴近实际军事应用的目标仍然迫在眉睫。

发明内容

本发明的目的在于提供一种基于GPU加速的量子雷达散射截面计算方法。

实现本发明目的的技术解决方案为：一种基于GPU加速的量子雷达散射截面计算方法，步骤如下：

步骤1、对电磁场进行量子化，根据二能级自发辐射原理，求得光子状态、光子波函数，并得到量子雷达散射截面计算公式；

步骤2、对目标建模并进行表面三角形面元离散，得到目标剖分文件以及每个三角形面元的中心点以及外法向量，并根据目标剖分文件进行八叉树分组；

步骤3、对目标进行遮挡以及亮暗面判断并记录亮面三角形面元对应的顶点坐标，计算目标在入射方向上的正交投影面积；

步骤4、利用步骤3所获得的数据，计算量子雷达散射截面。

本发明与现有技术相比，其显著优点为：(1)本发明将计算任意复杂三维目标入射方向正交投影面积具体化，只需给出入射方向矢量以及亮面三角形面元的法向量，就可以得到其正交投影面积，计算简单有效；(2)本发明引入GPU加速技术，对于一些大型复杂目标，能够快速获得其量子雷达散射特性，提升计算效率。

下面结合附图对本发明作进一步详细描述。

附图说明

图1是一种基于GPU加速的雷达散射截面计算方法的流程图。

图2是本发明实施例1中立方体模型示意图。

图3是本发明实施例1中立方体模型单站QRCS与RCS曲线对比图。

图4是本发明实施例2中B2飞机模型示意图。

图5是本发明实施例2中B2飞机模型单站QRCS与RCS曲线对比图。

具体实施方式

一种基于GPU(图形处理器)加速的量子雷达散射截面计算方法，首先对电磁场进行量子化，根据二能级自发辐射原理，求得光子状态，从而求得光子波函数，类比经典雷达理论，推导得到量子雷达散射截面积计算公式；然后目标建模并进行表面三角形面元离散，并根据目标剖分文件进行八叉树分组；接下来对目标进行遮挡以及亮暗面判断，计算目标在入射方向上的正交投影面积；最后根据上一步获得的信息计算量子雷达散射截面积。

如图1所示，一种基于GPU加速的量子雷达散射截面计算方法，具体步骤为：

步骤1、对电磁场进行量子化，根据二能级自发辐射原理，求得光子状态，从而求得光子波函数，并得到量子雷达散射截面计算公式。具体步骤如下：

步骤1.1、电磁场的正则量子化。为了揭示电磁场的波粒二象性，充分反映电磁场与物质相互作用的量子特性，必须将电磁场进行量子化。在经典电磁场中，在体积为V的空腔内，矢量势可以表示为：

相应的电场强度和磁场强度为：

其中，i表示虚数单位，r表示空间位置，t表示时间，k为波矢量，e_k为电磁场的极化单位矢量，表征电磁场的偏振方向，A_k(t)表示电磁波的振幅，

表示对A_k(t)取共轭。

应用正则量子化方法，引入产生算符

和湮灭算符

则A_k(t)、

分别表示为：

得到量子化之后的矢量势：

其中，

为约化普朗克常数，ε₀为自由空间介电常数，ω_k为角频率。

在自由场情况下，湮灭算符满足运动方程：

由此得到，

将公式(4)、(7)代入式(2)和(3)得到量子化之后的电磁场：

通常情况下，场算符可以写成正、负频两部分：

可见正频部分只包含湮灭算符，负频部分只包含产生算符。

步骤1.2、利用二能级自发辐射原理，得到光子状态以及光子波函数。

二能级原子自发辐射是量子电动力学中最基本的作用过程之一。二能级原子只有两个非简并能级E₊和E_-，其本征跃迁频率为

根据二能级自发辐射原理得到光子状态为：

其中r₀表示二能级原子位置，下标k代表辐射光子波矢量，ω_k为辐射光子角频率，Γ为衰变速率，|1_k>表示产生一个光子，g_k表示耦合系数，具体为：

根据光子波函数定义以及式(8)和(11)，得到光子波函数的表达式为：

其中，c为光速，t表示时间，ΔR₀＝r-r₀表示光子探测位置和二能级原子之间的位置矢量。

步骤1.3、类比于传统雷达散射截面(RCS)的定义，量子雷达散射截面表示为：

其中

和

分别表示散射场和入射场强度，下标i和s分别代表入射和散射，r_d和r_s分别为量子雷达的接收机位置和发射机位置，R为雷达作用距离。

入射到目标处的光子将以一定的概率沿各个方向散射出去，则雷达接收机处检测到的散射场强度为：

其中，η为自由空间波阻抗，N为目标表面原子总个数，ΔR_m为总的干涉距离(即雷达发射机到目标以及目标到雷达接收机的距离)。

当忽略吸收效应时，体系能量守恒原则要求目标截获的所有能量将全部向半球空间中散射出去，则得到入射场强度与散射场强度对应关系为：

其中，θ_i,φ_i分别表示入射光子的俯仰角和方位角，A_⊥为目标对于入射波矢量的正交投影面积，其大小依赖于目标实际面积以及入射光子的俯仰角和方位角。

则将量子雷达散射截面表达式写为：

其中，k_i和k_s分别表示入射光子波矢量和散射光子波矢量，x_m表示目标上原子的位置。

步骤2、目标建模并进行表面三角形离散，计算得到每个三角形面元的中心点、外法向量以及目标模型剖分文件，每个三角形面元的外法向量具体为：

其中，

分别表示沿着三角形面元两个边a,b的方向矢量。

根据目标剖分文件进行八叉树分组，具体步骤如下：

步骤2.1、首先用一个立方体包围目标模型，将该立方体定义为第0层，将立方体等分为八个小立方体，将小立方体定义为第1层，每个小立方体再继续等分位八个小立方体，直到第L层，将第l层的每个立方体定义为一组，第l层所含有的组数为8^l，使每个组中的平均离散边的个数不超过50个；

步骤2.2、将每一层中的所有组按照组中心的位置依次编号为1到8^l，其中1≤l≤L；

步骤2.3、将第l-1层定义为第l层的父层，第l层定义为第l-1层的子层，则子层中组i对应的父层中的组编号定义为i^p，由子层中的组i索引它的父层中的组i^p的方法为首先把编号i转化成二进制序列，把该二进制序列去掉右边三位并且转化成十进制即为编号i^p。

步骤3、对目标进行遮挡以及亮暗面判断，并记录亮面三角形面元对应的顶点坐标，计算目标在入射方向上的投影面积，具体步骤如下：

步骤3.1、判别三角形面元的亮暗面并记录亮面三角形面元对应的顶点坐标，判别方法为：

面元位于照亮区

面元位于阴影区

为入射光子波矢量，

为目标三角形面元的单位外法向分量；

步骤3-2、对目标模型进行遮挡判断。从三角形面元S的中心点出发，沿着入射光子波矢的反方向发射一条测试射线，当测试射线与目标表面三角形面元没有交点时，则表明该三角形面元没有被遮挡，否则该面元被遮挡处于阴影区。其中利用测试射线与三角形面元所在平面进行求交，判断测试射线与目标三角形有无交点，这一过程较为耗时，所以采用了八叉树算法来加速，即先计算测试射线与第一层中的每个组的组中心的距离，如果距离小于该组的对角线长度的一半，则认为测试射线有可能和该组内的三角形面元有交点，通过八叉树索引到该组的下一层，如此往复。最后，将最细层空间组内的面元与测试射线方程联立求出交点。整个测试过程，剔除了许多不相交的空间组的面元，大大减少了求交的面元数量，从而减少了计算量，提高计算速度。

步骤3.3、设每个三角形面元的面积为A_i，则可以表示为：

则正交投影面积

其中，T表示总的亮面的个数，

表示入射波矢，

表示每个三角形面元的外法向量。

步骤4，将获得的亮面上的点的坐标代入相应公式中，计算量子雷达散射截面。具体步骤如下：

步骤4.1、根据步骤3，将量子雷达散射截面积计算公式重新写为：

步骤4.2、对于公式(18)中的积分项，先将被积函数在积分区间内离散，因为是二重积分，可以认为是填写一个二维矩阵M，分别用p,q来表示矩阵的行和列，矩阵的每一个元素记为M_pq，矩阵的每一个元素M_pq通过GPU线程计算，计算公式为：

其中，

K＝2π/λ，λ为波长；

GPU中每个线程计算矩阵M中元素M_pq对应的索引为：

p＝blockdim％x*(blockidx％x-1)+threadidx％x

q＝blockdim％y*(blockidx％y-1)+threadidx％y

其中，blockdim表示一个线程块中包含的线程的个数，blockidx代表线程块在线程网格中的标号，threadidx表示一个线程在其线程块中的标号。

步骤4.3、利用梯形积分法对矩阵M进行离散求和积分，得到量子雷达散射截面。

实施例1

为了验证本发明的正确性与有效性，下面给出数值算例。所有算例均在主频2.83GHz、内存8GB的个人计算机上实现。GPU型号为NVIDIA GeForce GTX 1080 Ti。

考察一个立方体模型，结构如图2所示：立方体尺寸为4λ×4λ×4λ，λ＝0.25m，入射波频率为1.2GHz，扫描角度为

θ^inc＝0°～90°。图3为立方体模型单站QRCS(该发明计算结果)与RCS曲线对比图，从图中可以看出，本发明具有很高精度的同时具有极高的效率，并且能够高效快速分析大型复杂目标的量子雷达散射特性。

如表1所示为立方体模型为4λ×4λ×4λ其中λ＝0.25m，其中当θ＝0°时，未知量个数为9380658，θ＝45°时，未知量个数为18761127。从表1中可以看出，利用GPU加速之后，能够大幅节省计算时间，使得分析更大更为复杂的目标的量子雷达散射特性成为可能。

表1

求解用时	θ＝0°	θ＝45°
			CPU并行计算(16个核)	14669.26s	31133.01s
GPU并行计算	95.41s	190.63s

实施例2

图4考察一B2飞机模型，该模型翼展约52.4m，机身长约20.5m，机身高约3m。入射波频率为300MHz,扫描角度为

θ^inc＝0°～180°。图5为B2飞机模型单站QRCS(本发明方法计算结果)分别与物理光学法(PO)以及矩量法(MoM)雷达散射截面积曲线对比图，所得QRCS曲线与传统方法计算结果趋势保持一致。该目标较为复杂且体型较大，且未知量达到一亿之多，更进一步表明这种基于GPU加速的方法能够快速分析电大复杂目标的量子雷达散射特性。

Claims (8)

1.一种基于GPU加速的量子雷达散射截面计算方法，其特征在于，步骤如下：

步骤1、对电磁场进行量子化，根据二能级自发辐射原理，求得光子状态、光子波函数，并得到量子雷达散射截面计算公式，具体步骤为：

步骤1-1、对电磁场进行正则量子化得到量子化后的矢量势以及量子化之后的电磁场，并将量子化后的电场强度表示为正负频两部分，其中，量子化后的矢量势为：

式中，

为约化普朗克常数，ε₀为自由空间介电常数，ω_k为角频率，

为产生算符，

为湮灭算符，e_k为电磁场的极化单位矢量，i表示虚数单位，r表示空间位置，t表示时间，k为波矢量；

量子化之后的电磁场为：

式中，V为空腔体积，

为量子化后的电场强度，

为量子化后的磁场强度；

将量子化后的电场强度表示为正负频两部分为：

步骤1-2、根据二能级自发辐射原理，得到光子状态以及光子波函数，其中，所述光子状态为：

式中，ω₀为本征跃迁频率，具体为：

E₊和E_-为二能级原子的两个非简并能级，r₀表示二能级原子位置，下标k代表辐射光子波矢量，ω_k为辐射光子角频率，Γ为衰变速率，|1_k>表示产生一个光子，g_k表示耦合系数；

所述光子波函数的表达式为：

式中，c为光速，t表示时间，ΔR₀＝r-r₀表示光子探测位置和二能级原子之间的位置矢量；

步骤1-3、根据光子状态以及光子波函数，确定量子雷达散射截面表达式：

式中，k_i和k_s分别表示入射光子波矢量和散射光子波矢量，x_m表示目标上原子的位置，A_⊥为目标对于入射波矢量的正交投影面积，N为目标表面原子总个数；

步骤2、对目标建模并进行表面三角形面元离散，得到目标剖分文件以及每个三角形面元的中心点以及外法向量，并根据目标剖分文件进行八叉树分组；

步骤3、对目标进行遮挡以及亮暗面判断并记录亮面三角形面元对应的顶点坐标，计算目标在入射方向上的正交投影面积；

步骤4、利用步骤3所获得的数据，计算量子雷达散射截面。

2.根据权利要求1所述的基于GPU加速的量子雷达散射截面计算方法，其特征在于，步骤2中三角形面元的外法向量为：

其中，

分别表示沿着三角形面元两个边a,b的方向矢量。

3.根据权利要求1所述的基于GPU加速的量子雷达散射截面计算方法，其特征在于，步骤2中根据目标模型剖分文件，对目标进行八叉树分组，具体操作如下：

步骤2.1、首先用一个立方体包围目标模型，将该立方体定义为第0层，将立方体等分为八个小立方体，将小立方体定义为第1层，每个小立方体再继续等分位八个小立方体，直到第L层，将第l层的每一个立方体定义为一个组，第l层所含有的组数为8^l，并使每个组中的平均离散边的个数不超过50个；

步骤2.2、将每一层中的所有组按照组中心的位置依次编号为1到8^l，其中1≤l≤L；

步骤2.3、将第l-1层定义为第l层的父层，第l层定义为第l-1层的子层，则子层中组i对应的父层中的组编号定义为ip，由子层中的组i索引它的父层中的组ip。

4.根据权利要求3所述的基于GPU加速的量子雷达散射截面计算方法，其特征在于，索引方法为：把编号i转化成二进制序列，把该二进制序列去掉右边三位并且转化成十进制即为编号ip。

5.根据权利要求1所述的基于GPU加速的量子雷达散射截面计算方法，其特征在于，步骤3中对目标模型进行遮挡以及亮暗面判断，计算目标在入射方向的正交投影面积的具体方法为：

步骤3.1、判别三角形面元的亮暗面并记录亮面三角形面元对应的顶点坐标，判别方法为：

面元位于照亮区

面元位于阴影区

为入射光子波矢量，

为目标三角形面元的单位外法向分量；

步骤3.2、从三角形面元的中心点出发，沿着入射光子波矢的反方向发射一条测试射线，当测试射线与目标表面三角形面元没有交点时，则表明该三角形面元没有被遮挡，否则该面元被遮挡处于阴影区；

步骤3.3、计算目标在入射方向上的正交投影面积，计算公式为：

式中，T表示总的亮面的个数，

表示入射波矢，

表示每个三角形面元的外法向量，A_i为每个三角形面元的面积，具体为：

6.根据权利要求5所述的基于GPU加速的量子雷达散射截面计算方法，其特征在于，判断测试射线与目标三角形有无交点的具体方法为：

步骤3.2.1、计算测试射线与第一层中的每个组的组中心的距离，如果距离小于该组的对角线长度的一半，则认为测试射线和该组内的三角形面元有交点；

步骤3.2.2、通过八叉树索引到该组的下一层，对下一层的每个组重复步骤3.2.1，直到第L层；

步骤3.2.3、将第L层立方体组内的所有面元与测试射线方程联立求交点，并判断交点在不在三角形面元内部，若在，则判断测试射线与目标表面三角形面元有交点。

7.根据权利要求1所述的基于GPU加速的量子雷达散射截面计算方法，其特征在于，步骤4利用步骤3所获得的数据，计算量子雷达散射截面的具体方法为：

步骤4.1、根据目标在入射方向上的正交投影面积，将量子雷达散射截面积计算公式重写为：

步骤4.2、对于重写的量子雷达散射截面积计算公式的积分项，将被积函数在积分区间内离散，填写为一个二维矩阵M，分别用p,q来表示矩阵的行和列，矩阵的每一个元素为M_pq；

步骤4.3、利用梯形积分法对矩阵M进行离散求和积分，得到量子雷达散射截面。

8.根据权利要求7所述的基于GPU加速的量子雷达散射截面计算方法，其特征在于，矩阵的每一个元素M_pq通过GPU线程计算，计算公式为：

其中，

K＝2π/λ，λ为波长；

GPU中每个线程计算矩阵M中元素M_pq对应的索引为：

p＝blockdim％x*(blockidx％x-1)+threadidx％x

q＝blockdim％y*(blockidx％y-1)+threadidx％y

其中，blockdim表示一个线程块中包含的线程的个数，blockidx代表线程块在线程网格中的标号，threadidx表示一个线程在其线程块中的标号。

Images