瞬态矢量辐射传输理论箔条云散射处理方法与流程

本发明属于箔条云时域散射计算技术领域，尤其涉及一种瞬态矢量辐射传输理论箔条云散射处理方法。

背景技术：

目前，最接近的现有技术：在现代电磁战争中，箔条因其历史悠久、制作简单、作用明显、成本低廉，被广泛应用。在电子对抗中作为主要侦察装备的雷达，在各国防御中发挥着至关重要的作用。因而现代战争特别是在未来战争中，箔条干扰技术通过在空间中释放数量庞大的箔条，致使敌方雷达在获得的回波与目标真实的回波之间的差异较大，从而敌方雷达不能准确和快速识别目标，这种干扰方式不需要预先了解敌方雷达的特性，通用性较强，其对雷达的干扰主要有两种方式，一是通过干扰物自身产生回波、掩盖目标的真实回波，二是改变干扰物空间位置和结构，使它的回波与真实的目标回波相似，来诱骗敌方雷达对目标的侦查判断。箔条弹由发射至空中炸开到形成箔条云分为非成熟期和成熟期。非成熟期是指箔条从箔条弹中炸开但还没有完全扩散开的这一段时期。在这一时期，箔条密度急剧减小，雷达后向散射截面急剧增大，其空间分布是不均匀和时变的，不具有干扰能力，成熟期是指箔条充分扩散后无论密度还是雷达后向散射截面都达到相对稳定的时期，成熟期的箔条云才具有防护目标的能力，且先进的箔条弹从爆炸到完全扩散开只需要30～50，而成熟的箔条云可以在空中滞留几分钟，所以箔条云大部分时间处于成熟阶段。对于成熟阶段内的箔条云的时域散射特性计算的研究对箔条云干扰和抗箔条云干扰的机理研究至关重要，已有的文献对于箔条云时域散射特性的研究大多数基于实验的实测数据进行，而通过计算机仿真获得数据是一种经济、可行的途径。

目前，对箔条云时域散射的计算主要使用半实物仿真方法、数值算法、独立散射方法、基于单箔条的散射截面的统计叠加方法。半实物仿真方法采用部分实际数据作为仿真的一部分，能够获得更加贴近实际情况的模拟结果，是一种比较可靠的仿真技术，从提出开始，就被广泛应用在系统研制和模拟实验中；数值算法(矩量法等)是求解maxwell方程组的积分或微分形式实现箔条云团散射的计算；独立散射方法是基于单根箔条的全空间散射特性，通过旋转变换和非相干统计叠加实现箔条云团散射计算，当然，也有很多研究是基于实测结果展开进一步分析，通过大量实验实现箔条云团散射的预估和特性研究。

综上所述，现有技术存在的问题是：

(1)基于数值算法的箔条云时域散射计算虽然精确，但由于计算资源的限制和需要求解的未知数数量巨大，致使数值算法无法计算大规模箔条云的时域散射。

(2)基于实测结果的箔条云团时域散射计算具有一定的现实意义，却无法穷尽除实验具体条件外的其他可能情况，实测结果有一定局限性，很难进行推广。

(3)基于独立散射的箔条云团时域散射计算虽然预估速度很快，但由于箔条之间存在电磁耦合，而独立散射方法却忽略箔条间的多次耦合，使之计算结果的可信度较低。

解决上述技术问题的难度：为了解决上述技术问题，主要有以下技术难点：如何在考虑箔条间耦合的情况下实现快速时域散射计算是一个主要的难题。

解决上述技术问题的意义：实现时域箔条云团散射的快速计算，可通过计算结果分析箔条云的时域信息，同时为大规模电小散射体电磁计算提供一种基于时域的高效计算方法。

技术实现要素：

针对现有技术存在的问题，本发明提供了一种瞬态矢量辐射传输理论箔条云散射处理方法。

本发明是这样实现的，一种瞬态矢量辐射传输理论箔条云散射处理方法，所述瞬态矢量辐射传输理论箔条云散射处理方法包括以下步骤：

第一步，对入射波——线性调频信号进行分段傅立叶变换，获取频率信息；

第二步，基于箔条云空间分布模型使用k-means算法进行分区处理，并使用线矩量法计算瞬态矢量辐射传输理论所需的分区衰减系数信息；

读取箔条云团空间分布数据，并基于k-means算法对箔条云团进行区域划分包括：

(1)调用箔条云的空间分布及取向信息，数据的读入格式为：

其中，(x,y,z)是箔条中心点的坐标，单位为m，是箔条的取向，单位为弧度；

(2)箔条相互之间的距离对箔条云进行分类，间接反映箔条云不同区域的分布数密度，箔条云分区包括：1)在箔条中心点模型中，随机选择k个中心点；2)遍历所有箔条云空间中心点数据，将每个箔条中心点归属为k个拟定中心点所对应的类别中；3)计算k个聚类中每一个子类的重心，并作为新的k个分类的中心点；4)重复2)、3)，直至k类中心点不再变化或小于指定阈值；将箔条云进行分区，分区后数据存储格式为：

(x′,y′,z′,n)；

其中，(x′,y′,z′)是分区后的箔条中心点的坐标，n为分区编号；

第三步，基于分区数据和获得的区域内能量衰减系数，使用蒙特卡洛方法对瞬态矢量辐射传输方程求解，使用光子发射、输运、碰撞、接收过程追踪的方法实现瞬态矢量辐射传输方程的并行快速求解；

第四步，基于瞬态矢量辐射传输方程的求解结果，经过数据处理得到stokes矢量四个分量随时间变化结果。

进一步，所述第一步包括：对入射波——线性调频信号进行分段处理，并做傅立叶变换，获取其频率信息；

(1)线性调频信号是指瞬态频率随时间呈线性变化的信号，chrip信号的时域表达式写为：

其中，t是时间变量，单位为秒；t为脉冲持续时间；k是线性调频斜率，单位是hz/s；

角度表达式：

对时间取微分后的瞬态频率为：

信号的带宽是chrip信号的斜率和时间的乘积：

bw＝|k|t；

带宽决定了能够达到的分辨率；

(2)基于线性调频信号的瞬态频率随时间呈线性变化的情况，对线性调频信号进行分段处理，并利用傅立叶变换fft将信号分解成频率谱进而求得各分段信号的中心频率，傅立叶变换公式如下：

(3)利用定积分法计算入射线性调频信号的概率分布与累计概率分布，把直角坐标系上的函数的图像用平行于y轴的直线把其分割成无数个矩形，把某个区间[a,b]上的矩形累加起来，所得到的就是函数的图像在区间[a,b]的面积，定积分法的表达式如下：

对于所有实数x，累积分布函数定义如下：

fx(x)＝p(x≤x)；

将累积分布函数应用到分段入射波信号的概率分布的求解中，概率的获得是利用定积分求得的分段信号的面积与整个时段内函数的面积比值求得。

进一步，所述第三步使用随机抽样算法，在每个分区使用离散事件抽样算法抽取相邻少量箔条利用线矩量法计算电磁波的极化传输损耗，建立数据库，包括以下步骤：

(1)在某一分区中，使用均匀分布随机数选中该区域中某一箔条：

n(i)＝unrand(1,mi)；

其中，函数unrand(1,mi)依均匀分布随机产生大于等于1且小于等于mi的正整数n(i)，标记i是区域的索引号，mi表示第i个区域箔条总数量；

(2)在某一分区中，同样使用k-means算法寻找与第n(i)号箔条最为临近的q根箔条，并使用线矩量法计算被选择的q根箔条电磁散射结果，每根箔条被分成t段，在每一段上定义矢量基函数：

其中，xm为第m段的中点，δxm为第m段的长度，所有箔条电流表示为：

其中，αn是待定系数，定义在第n段之上，在其他段为零，当箔条是理想导体时，在箔条表面建立如下电场积分方程：

其中，和分别为场点和源点，为入射波电场强度，将函数作为未知函数f代入泛函方程：

l(f)＝g；

得到：

其中，l是积分方程线性算子，使用伽略金法，将权重函数wn作为检验函数，与方程两边进行内积运算，得到：

建立矩阵方程：

[lmn][αn]＝[gm]；

其中，

求解上述矩阵方程，得出αn的解，得到q根箔条表面的电流分布函数的具体表达式，使用：

求得散射远场，并确定q根箔条的远场雷达散射截面：

(3)在获取q根箔条的雷达散射截面的基础上计算损耗系数κ：

κ＝-nσ；

其中，n是单位体积箔条的数密度；

(4)对每个分区进行delaunay三角剖分和外包络建立，并针对不同区域的外包络每个三角形顶点、法向量、中心点索引建立数据库。

进一步，针对不同区域的外包络每个三角形顶点、法向量、中心点索引建立数据库包括以下步骤：

(1)采用bowyer逐点插入法实现delaunay三角剖分；将点集内所有的点包含在一个超级三角形中进行处理，实现箔条云每个分区的delaunay三角剖分；

1)在箔条空间分布点中取相邻三个点a、b、c，在点a、b、c外面建立一个超级三角形pqr；

2)分析a点，因为a点在超级三角形pqr内部，把a点和p、q、r三点分别连线；

3)再考虑b点，看b点在哪个子三角形的内部，结果发现只在aqr三角形内部，于是b点分别和aqr这个三角形的三个顶点连线；

4)此时一共有五个三角形，分别画出这五个三角形的外接圆，对c点检查是否在这五个三角形的外接圆内；

5)判定结果可知c点在三角形apr和三角形abr的外接圆内，删除这两个三角形的公共边ar，然后从c点分别向这两个三角形组合起来的四边形的四个顶点连线；

6)最后把含有超级三角形pqr顶点的三角形全部删除，就能得到确定的delaunay三角形；

(2)剖分寻找外包络点；

1)依据给定的箔条云空间数据，寻找三个维度的最大值和最小值：

xmax，xmin，ymax，ymin，zmax，zmin；

并在空间使用长方体进行剖分，每个长方体的三个方向上的边长为(δx,δy,δz)，于是，每个长方体区域的中心点为：

且

其中，ceil()函数表示向上取整数；

2)使用收缩算法，自外向内搜索每一层盒子，遇到盒子有点存在时停下；

(3)确定外包络三角剖分外法线方向，使用连续临边外法向方向计算准则，两相邻三角形的顶点分别为a，b，c和b，c，d；m，n，p是三角形abc三个顶点序号，w，g，h为三角形bcd的三个顶点，按照右手螺旋定则和从小到大的顺序计算三角形abc法向量为：

三角形bcd的三个顶点的标号自小到大依次为bdc、cbd、dcb，继续下去，编号所有的三角形顶点，则在循环每一个三角形时，按照右手定则确定外法线方向，对于外包络中心点坐标使用三顶点坐标的代数均值计算。

进一步，所述使用收缩算法，自外向内搜索每一层盒子，遇到盒子有点存在时停下包括：

1)对于盒子序号为n(n,:,:)区域1：如果(n-1,m,p)盒子中已经存在标志位为1的点，则后续(:,m,p)盒子不再进行判断；否则，循环(m,p)，如果盒子n(n,m,p)不为0，则该盒子内最外侧点标记为1；

2)对于盒子序号为n(:,m,:)的区域2：如果盒子(n,m-1,p)中已经存在标志位为1的点，则后续(n,:,p)盒子不再进行判断；否则，循环(n,p)，如果盒子n(n,m,p)不为0，则该盒子内最外侧点标记为1；

3)对于盒子序号为n(:,:,p)的区域3：如果盒子(n,m,p-1)中已经存在标志位为1的点，则后续(n,m,:)盒子不再进行判断；否则，循环(n,m)，如果盒子n(n,m,p)不为0，则该盒子内最外侧点标记为1。

进一步，所述第三步包括：建立瞬态矢量辐射传输方程，并使用蒙特卡洛方法模拟光子在箔条云团中的碰撞过程，求解时域方程包括：

(1)考虑天线方向图信息的源平面产生；

(2)源发射光子与箔条云外三角剖分面相交判，抽取源点：

p(m)＝(xsm,ysm,zsm)；

并已知入射波的传播方向为：

则写出粒子运动的直线方程为：

x＝xsm+l·nx，y＝ysm+l·ny，z＝zsm+l·nz；

其中，l为参数；

判断某外包络三角形与直线是否相交的方法为，需要判交的三角形的三个顶点可记为tr1，tr2，tr3，则三角形内部的点表示为：

将已经推导的直线表达式代入得到：

λ1(tr1-tr2)+λ2(tr1-tr3)+tnxnynz＝tr1-pm)；

简写为：

λ1a+λ2b+lc＝d；

写成矩阵形式为：

其解为：

如果λ1,λ2∈[0,1]，则射线和箔条云外包络有交点，当有多个交点时，需要使用最小l法则判断第一交点；

根据计算出的第一交点坐标(x1,y1,z1)，求出单个光子发射源点——(xsm,ysm,zsm)与第一交点之间的距离lb，表达式为：

lb＝[(xsm-x1)²+(ysm-y1)²+(zsm-z1)²]^1/2；

通过lb的求解得出光子从源点发射到进入箔条云团之前所经历的时间tb：

c是真空中的光速，其值为3×10⁸m/s；

(3)光子在箔条云中依概率传输，使用蒙特卡洛方法求解瞬态矢量辐射传输方程，n根箔条随机分布在长为ds，体积为dv的柱形体单元中，一辐射强度为i的电磁波通过该单元，由能量守恒定律得到电磁波辐射强度的变化量为：

其中，是入射强度和散射强度之间的过渡系数矩阵，第一项指的是由于箔条吸收κas、散射κs、背景吸收κab导致的能量耗散，第二项表示的是空间其他辐射源辐射的能量，第三项表示的是多散射体耦合散射在方向上的散射强度总和，其中；

如果在空间中，箔条的占空比为fs，则背景的吸收系数表示如下：κab＝2k″·(1-fs)，k″是波数的虚部箔条吸收率是箔条损耗场能和入射场能的比，使用损耗场表示为：

引入stokes矢量的四个分量如下：

(4)光子出射与能量收集，产生光子、光子与散射体碰撞并追踪光子的轨迹与光子经历这些轨迹所用的时间，当光子离开散射体后将对散射能量与时间进行统计；有n个能量为i0i的光子被发射，有m个光子经连续碰撞逃出云团，散射方向分别为自身携带能量分别为ii，光子经历的时间为t′i，对全空间的对进行网格划分，θ角的步长为δθ，角的步长为同时有θ＝mδθ、判断出射光子是否进入第(m,n)单元的方法是：

于是，在方向接收到的总光子能量为：

而第i个光子遍历整个过程所经历的时间则由t′i表示：

t′i＝t0+tb+tl；

其中t0表示光子在源平面出射时的初始时间，tb表示光子从出射到进入箔条云之前所经历的时间，tl表示光子进入箔条云，在其中折射和反射后经历自由程，最终出射云边界所用的时间。

进一步，所述考虑天线方向图信息的源平面产生包括以下步骤：

1)按照外部给定的天线，读入天线的远场方向图信息，数据存储格式为：

其中，θ、是天线远场方向角，etheta和ephi是远场指定位置电场的方向和方向分量；

2)在给定的电磁波入射角情况下，建立源平面的方程如下：

其中，源平面的中心点为：

获得给定箔条云团的外包络信息及其剖分数据，将外包络点在平面π上进行投影，并再次使用delaunay三角剖分给出平面外包络；

3)使用矩形网格，及得到的平面外包络产生源平面上的采样点，在源平面上定义坐标系如下：

并单位化：

记外包络的各顶点的位置矢量为(xi,outline,yi,outline,zi,outline)，在定义的矢量上投影，在和维度平面，使用矩形网格策略产生网格如下：

其中δvn，δwm为划分的网格大小，同时需要判定，如果产生的点(vn,wn)不在最大外包络内，则去除；

4)依据方向函数和离散事件抽样方法进行源发射点抽样，在指定的外包络内使用矩形网格产生的离散点为p(i)，每个点对应的归一化方向图值为q(i)，则产生源抽样；

计算累积概率：

按照均匀分布产生一个随机数：

randnumber＝rand()；

其中，函数可产生[0,1]区间内的随机数。若m(m)≤randnumber≤m(m+1)，则被抽取的源点序号为m，其对应的点坐标为p(m)；

5)利用[0,1]区间内的随机数得到每一分段信号内光子的初始时间t0，具体表达式如下：

t0＝ti+rand()·δt；

其中ti为第i段信号的起始时间，rand()为[0,1]区间内的随机数，δt为分段信号的时间长度；并根据时间计算出该光子的初始能量i0

i0＝|s(t)|·sin·iapm；

其中，s(t)为chrip信号表达式，sin为入射的stokes矢量，而iapm则表示源发射点位置处的方向图表征的能量值。

进一步，所述基于蒙特卡洛方法求解瞬态矢量辐射传输方程包括以下步骤：

1)自由程抽样，得到光子第一次进入箔条云的具体位置，并按照下面的自由程公式进行抽样，得到下次碰撞光子前所走路程：

其中，ξ是[0,1]之间的均匀分布随机数，ke是消光系数；对于时域问题，自由行程的距离意味时间的消耗，流逝的时间取决于自由行程的距离以及光在介质中的光速；

传播l距离所耗费的时间为：

其中n为介质的折射率，c0为真空中的光速，则光子在传播l距离后时间变为：

t′＝t0+tb+tl；

2)箔条取向抽样，箔条的取向由空间角确定，服从的概率密度函数为p(θ)和使用舍选法抽取散射方向；

首先，在θ∈[0,π]，之间按照均匀分布抽取箔条取向角

将箔条取向角带入概率密度函数

其次，在区间上使用上述的ξ再次产生一个均匀分布随机数p0，按照判断准则进行发射角选择：

如果判断准则给出的是拒绝选择则需要返回在θ∈[0,π]，之间按照均匀分布抽取箔条取向角继续进行，直到产生满足要求的数对为止；

3)散射方向抽样，基于能量守恒和舍选法进行散射方向抽样，在某一确定的入射角下，归一化的散射能量表示为：

且：

则粒子沿着方向出射的概率密度为：

其中

其中，方向散射系数定义如下：

散射方向由空间角确定，的选取采用舍选法进行；

4)散射传输能量获取，计算每次碰撞的衰减因子在给定的入射波条件下，有：

单次碰撞前的能量wi-1和碰撞后的能量wi之间关系为：

直到光子射出箔条云，遍历设定的最大光子数进行数据处理，判断光子是否在云内外的基本判定方法是，从该点向任意方向发射射线，判断射线与几何体的交点个数，奇数为云内，偶数则为云外；单个光子在箔条云中经历的时间是通过自由程的叠加得到的，光子每一次碰撞，就会经历一个自由程l，直至最终光子从箔条云中逃逸或光子携带能量在碰撞过程损耗殆尽，单个光子经历的所有自由程叠加的结果与光在箔条云中的速度的比值就得到了光子在箔条云中经历的时间tl。

进一步，所述第三步还包括：对求解数据进行运算与处理，得到时域stokes矢量结果，stokes矢量的定义及入射光子能量i0i的表示形式，在方向的雷达散射截面表示为：

按照stokes矢量四个分量形式表示，方向的时域散射结果已通过运算得到，而方向收集的光子携带的时间信息通过t′i反映，通过光子的索引信息i将t′i与σi一一对应，求得σi随t′i变化的波形即为所求的箔条云的时域stokes矢量结果。

本发明的另一目的在于提供一种应用所述瞬态矢量辐射传输理论箔条云散射处理方法的智能终端。

综上所述，本发明的优点及积极效果为：本发明基于并行瞬态矢量辐射传输理论(transientvectorradiativetransfertheory，tvrt)的大体量箔条云电磁散射计算方法，主要使用蒙特卡洛方法(montecarlomethod，mc)，解决电子对抗领域大数量的箔条云时域斯托克斯矢量(stokes矢量)快速高效计算问题。本发明具体涉及线性调频信号的处理方法、箔条云团的分区方法、时域散射计算方法，可用于电子对抗领域箔条云团时域stokes矢量结果计算及使用方法探索。

本发明基于瞬态矢量辐射传输理论，可以实现箔条云团电磁散射的快速准确计算，在实际应用中可实现快速实时计算。本发明使用线性调频信号作为入射波信号，基于该信号瞬态频率随时间变化，在仿真实现中提取了chrip信号的不同频率分量，针对不同频率对箔条云团的散射进行了计算，最后得到了针对入射波信号的时域stokes矢量结果，提高了仿真的准确性。

本发明基于delaunay算法对箔条云团进行分区处理，提升了模型的处理效率，并将射线判交方法引入算法之内，提高了散射算法的适应性和效率。本发明使用蒙特卡洛方法对瞬态矢量辐射传输方程进行求解，原理简单可操作性强，便于后续并行处理，在指定的仿真数目内可实现全时域stokes矢量结果获取。

附图说明

图1是本发明实施例提供的瞬态矢量辐射传输理论箔条云散射处理方法流程图。

图2是本发明实施例提供的瞬态矢量辐射传输理论箔条云散射处理方法的实现流程图。

图3是本发明实施例提供的线性调频信号预处理模块子流程图。

图4是本发明实施例提供的箔条云团分区模块子流程图。

图5是本发明实施例提供的分区损耗系数计算模块子流程图。

图6是本发明实施例提供的活动特征提取模块子流程图。

图7是本发明实施例提供的分区三角剖分和外包络建立模块子流程图。

图8是本发明实施例提供的基于瞬态矢量辐射传输理论和蒙特卡洛方法的单个光子散射计算模块子流程图。

图9是本发明实施例提供的散射计算坐标系示意图。

图10是本发明实施例提供的喇叭天线方向图。

图11是本发明实施例提供的喇叭天线与箔条云关系图。

图12是本发明实施例提供的时域stokes矢量i分量示意图。

图13是本发明实施例提供的时域stokes矢量q分量示意图。

图14是本发明实施例提供的时域stokes矢量u分量示意图。

图15是本发明实施例提供的时域stokes矢量v分量示意图。

图16是本发明实施例提供的时域stokes矢量i分量示意图。

图17是本发明实施例提供的时域stokes矢量q分量示意图。

图18是本发明实施例提供的时域stokes矢量u分量示意图。

图19是本发明实施例提供的时域stokes矢量v分量示意图。

图20是本发明实施例提供的时域stokes矢量i分量示意图。

图21是本发明实施例提供的时域stokes矢量q分量示意图。

图22是本发明实施例提供的时域stokes矢量u分量示意图。

图23是本发明实施例提供的时域stokes矢量v分量示意图。

图24是本发明实施例提供的时域stokes矢量i分量示意图。

图25是本发明实施例提供的时域stokes矢量q分量示意图。

图26是本发明实施例提供的时域stokes矢量u分量示意图。

图27是本发明实施例提供的时域stokes矢量v分量示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

针对现有技术存在的问题，本发明提供了一种瞬态矢量辐射传输理论箔条云散射处理方法，下面结合附图对本发明作详细的描述。

如图1所示，本发明实施例提供的瞬态矢量辐射传输理论箔条云散射处理方法包括以下步骤：

s101：对入射波——线性调频信号进行分段傅立叶变换，获取其频率信息；

s102：基于箔条云空间分布模型使用k-means算法对其进行分区处理，并使用线矩量法计算瞬态矢量辐射传输理论所需的分区衰减系数信息；

s103：基于分区数据和获得的区域内能量衰减系数，使用蒙特卡洛方法对瞬态矢量辐射传输方程进行求解，即，使用光子发射、输运、碰撞、接收过程追踪的方法实现瞬态矢量辐射传输方程的并行快速求解；

s104：基于瞬态矢量辐射传输方程的求解结果，经过数据处理得到stokes矢量四个分量随时间变化结果。

下面结合附图对本发明的技术方案作进一步的描述。

本发明已解决大体量箔条云团时域电磁散射计算问题；在对入射波信号——线性调频信号进行处理的基础上，调用箔条云团几何分布数据，利用k-means算法、delaunay算法对云团按照指定的分类标准进行分类与剖分，之后基于瞬态矢量辐射传输理论和蒙特卡洛方法模拟入射光子与箔条云的作用过程，最终经过运算得到箔条云时域stokes矢量结果。

如图2所示，本发明实施例提供的瞬态矢量辐射传输理论箔条云散射处理方法包括以下步骤：

(1)预处理，如图3所示，对入射波——线性调频信号进行分段处理，并对其做傅立叶变换，获取其频率信息。针对时域问题，入射波信号选取线性调频信号，仿真频率范围为4～7ghz。chrip信号可在增大射频脉冲宽度、提高平均发射功率、加大通信距离同时又保持足够的信号频谱宽度，不降低雷达的距离分辨率。

(1.1)线性调频信号是指瞬态频率随时间呈线性变化的信号，chrip信号的时域表达式可以写为：

其中，t是时间变量，单位为秒(s)；t为脉冲持续时间(周期)；k是线性调频斜率，单位是hz/s。

角度表达式：

对时间取微分后的瞬态频率为：

信号的带宽是chrip信号的斜率和时间的乘积：

bw＝|k|t；

带宽决定了能够达到的分辨率。

(1.2)基于线性调频信号的瞬态频率随时间呈线性变化的情况，对该信号进行分段处理，并利用傅立叶变换(fft)将信号分解成频率谱进而求得各分段信号的中心频率。傅立叶变换公式如下：

(1.3)利用定积分法计算入射线性调频信号的概率分布与累计概率分布。该方法把直角坐标系上的函数的图像用平行于y轴的直线把其分割成无数个矩形，然后把某个区间[a,b]上的矩形累加起来，所得到的就是这个函数的图像在区间[a,b]的面积。定积分法的具体表达式如下：

对于所有实数x，累积分布函数定义如下：

fx(x)＝p(x≤x)；

将累积分布函数应用到分段入射波信号的概率分布的求解中，概率的获得是利用定积分求得的分段信号的面积与整个时段内函数的面积比值求得。

(2)如图4所示，读取箔条云团空间分布数据，并基于k-means算法对箔条云团进行区域划分，以备后续进行电参数获取，包含以下步骤：

(2.1)调用箔条云的空间分布及取向信息，数据的读入格式为：

其中，(x,y,z)是箔条中心点的坐标，单位为m，是箔条的取向，单位为弧度。

(2.2)k-means算法中k值一般根据工程实际选择，它决定聚类结果中类别的数量。依据是箔条相互之间的距离对箔条云进行分类，其间接反映箔条云不同区域的分布数密度，箔条云分区包含下面的几个步骤：

1)在箔条中心点模型中，随机选择k个中心点；

2)遍历所有箔条云空间中心点数据，将每个箔条中心点归属为k个拟定中心点所对应的类别中；

3)计算k个聚类中每一个子类的重心，并作为新的k个分类的中心点；

4)重复步骤2)、3)，直至k类中心点不再变化或小于指定阈值；这样便可按照要求将箔条云进行分区。分区后数据存储格式为：

(x′,y′,z′,n)；

其中，(x′,y′,z′)是分区后的箔条中心点的坐标，n为分区编号。

(3)如图5所示，使用随机抽样算法，在每个分区使用离散事件抽样算法抽取相邻少量箔条利用线矩量法计算电磁波的极化传输损耗(衰减系数)，针对其建立数据库，包含以下步骤：

(3.1)在某一分区中，使用均匀分布随机数选中该区域中某一箔条，即：

n(i)＝unrand(1,mi)；

其中，函数unrand(1,mi)可以依均匀分布随机产生大于等于1且小于等于mi的正整数n(i)，标记i是区域的索引号，mi表示第i个区域箔条总数量；

(3.2)在某一分区中，同样使用k-means算法寻找与第n(i)号箔条最为临近的q根箔条，并使用线矩量法计算被选择的q根箔条电磁散射结果，设每根箔条被分成t段，在每一段上定义矢量基函数；

其中，xm为第m段的中点，δxm为第m段的长度。于是所有箔条电流可以表示为：

其中，αn是待定系数，定义在第n段之上，在其他段为零，当箔条是理想导体时，在箔条表面可建立如下电场积分方程：

其中，和分别为场点和源点，为入射波电场强度，将函数作为未知函数f代入泛函方程：

l(f)＝g；

可以得到：

其中，l是积分方程线性算子。使用伽略金法，将权重函数wn作为检验函数，与上述方程两边进行内积运算，可以得到：

在此基础上可以建立矩阵方程：

[lmn][αn]＝[gm]；

其中，

求解上述矩阵方程，得出αn的解，进而可以得到q根箔条表面的电流分布函数的具体表达式，在此基础上使用：

求得散射远场，并确定q根箔条的远场雷达散射截面：

(3.3)在获取q根箔条的雷达散射截面的基础上可计算损耗系数κ如下：

κ＝-nσ；

其中，n是单位体积箔条的数密度。其余分区均使用上述方法进行计算即可。

(4)如图6所示，对上述每个分区进行delaunay三角剖分和外包络建立，并针对不同区域的外包络每个三角形顶点、法向量、中心点索引建立数据库，具体包含以下步骤：

(4.1)delaunay三角剖分算法基本实现步骤：

delaunay三角剖分算法的优秀特性是空圆特性和最大化最小角特性，这两个特性避免了狭长三角形的产生。空圆特性其实就是对于两个共边的三角形，任意一个三角形的外接圆中都不能包含有另一个三角形的顶点，这种形式的剖分产生的最小角最大。采用bowyer逐点插入法实现，基本步骤如下(这里采用三个点的情况为例，多点情况类似)：

1)在箔条空间分布点中取相邻三个点a、b、c，在他们外面建立一个足够大的超级三角形pqr；

2)分析a点，因为a点在超级三角形pqr内部，把a点和p、q、r三点分别连线；

3)再考虑b点，看b点在哪个子三角形的内部，结果发现只在aqr三角形内部，于是b点分别和aqr这个三角形的三个顶点连线；

4)此时一共有五个三角形，分别画出这五个三角形的外接圆，对c点检查是否在这五个三角形的外接圆内；

5)判定结果可知c点在三角形apr和三角形abr的外接圆内，删除这两个三角形的公共边ar，然后从c点分别向这两个三角形组合起来的四边形的四个顶点连线；

6)最后把含有超级三角形pqr顶点的三角形全部删除，就能得到确定的delaunay三角形。

具体实施方法是将点集内所有的点包含在一个超级三角形中进行上述步骤处理，就可实现箔条云每个分区的delaunay三角剖分。

(4.2)进行上述剖分后，为了寻找外包络点(外包络三角形)，可以采用下面的算法步骤：

①依据给定的箔条云空间数据，寻找三个维度的最大值和最小值：

xmax，xmin，ymax，ymin，zmax，zmin；

并在空间使用长方体进行剖分，每个长方体的三个方向上的边长为(δx,δy,δz)，于是，每个长方体区域的中心点为：

且

其中，ceil()函数表示向上取整数。

②由于外包络点只可能存在长方体外侧的网格中，因此可以使用收缩算法，自外向内搜索每一层盒子，遇到盒子有点存在时停下。以最外层长方体盒子判断为例，有：

1)对于盒子序号为n(n,:,:)区域1：如果(n-1,m,p)盒子中已经存在标志位为1的点，则后续(:,m,p)盒子不再进行判断；否则，循环(m,p)，如果盒子n(n,m,p)不为0，则该盒子内最外侧点标记为1。

2)对于盒子序号为n(:,m,:)的区域2：如果盒子(n,m-1,p)中已经存在标志位为1的点，则后续(n,:,p)盒子不再进行判断；否则，循环(n,p)，如果盒子n(n,m,p)不为0，则该盒子内最外侧点标记为1。

3)对于盒子序号为n(:,:,p)的区域3：如果盒子(n,m,p-1)中已经存在标志位为1的点，则后续(n,m,:)盒子不再进行判断；否则，循环(n,m)，如果盒子n(n,m,p)不为0，则该盒子内最外侧点标记为1。

(4.3)确定外包络三角剖分外法线方向。使用连续临边外法向方向计算准则，设两相邻三角形的顶点分别为a，b，c和b，c，d。并设m，n，p是三角形abc三个顶点序号，w，g，h为三角形bcd的三个顶点，按照右手螺旋定则和从小到大的顺序计算三角形abc法向量为：

若三角形bcd的三个顶点的标号自小到大依次为bdc、cbd、dcb之一即可，如此继续下去，编号所有的三角形顶点，则在循环每一个三角形时，只要按照右手定则即可确定外法线方向，对于外包络中心点坐标可以使用三顶点坐标的代数均值计算。

(5)建立瞬态矢量辐射传输方程，并使用蒙特卡洛方法模拟光子在箔条云团中的碰撞过程，单个光子散射计算流程图如图7，光子与箔条云作用流程概念图如图8，求解时域方程，包含以下步骤：

(5.1)考虑天线方向图信息的源平面产生，该步骤是后续瞬态矢量辐射传输理论的前提和基础，主要包含以下步骤：

1)按照外部给定的天线，读入天线的远场方向图信息，数据存储格式为：

其中，θ、是天线远场方向角，etheta和ephi是远场指定位置电场的方向和方向分量。

2)在给定的电磁波入射角情况下，建立源平面的方程如下

其中，源平面的中心点可以选为:

上述的步骤(4)已经获得给定箔条云团的外包络信息及其剖分数据，将外包络点在平面π上进行投影，并再次使用delaunay三角剖分给出平面外包络。

3)使用矩形网格，及上一步得到的平面外包络产生源平面上的采样点，在源平面上定义坐标系如下：

并将其单位化：

记外包络的各顶点的位置矢量为(xi,outline,yi,outline,zi,outline)，将其在上述定义的矢量上投影，在和维度平面，使用矩形网格策略可以产生网格如下

m,n取0，±1，±2，±3，…

其中δvn，δwm为划分的网格大小，同时需要判定，如果产生的点(vn,wn)不在最大外包络内，则将其去除。

4)依据方向函数和离散事件抽样方法进行源发射点抽样，设在指定的外包络内使用矩形网格产生的离散点为p(i)，每个点对应的归一化方向图值为q(i)，则可以按照下面的方式产生源抽样。

计算累积概率：

按照均匀分布产生一个随机数：

randnumber＝rand()；

其中，函数可产生[0,1]区间内的随机数。若m(m)≤randnumber≤m(m+1)，则被抽取的源点序号为m，其对应的点坐标为p(m)。

5)利用[0,1]区间内的随机数得到每一分段信号内光子的初始时间t0，具体表达式如下：

t0＝ti+rand()·δt；

其中ti为第i段信号的起始时间，rand()为[0,1]区间内的随机数，δt为分段信号的时间长度；并根据时间计算出该光子的初始能量i0

i0＝|s(t)|·sin·iapm；

其中，s(t)为chrip信号表达式，sin为入射的stokes矢量，而iapm则表示源发射点位置处的方向图表征的能量值。

(5.2)源发射光子与箔条云外三角剖分面相交判定。上述步骤已经抽取源点：

p(m)＝(xsm,ysm,zsm)；

并已知入射波的传播方向为：

则可以写出粒子运动的直线方程为：

x＝xsm+l·nx，y＝ysm+l·ny，z＝zsm+l·nz；

其中，l为参数。

判断某外包络三角形与直线是否相交的方法为，设需要判交的三角形的三个顶点可记为tr1，tr2，tr3，则三角形内部的点可以表示为：

将已经推导的直线表达式代入可以得到：

λ1(tr1-tr2)+λ2(tr1-tr3)+tnxnynz＝tr1-pm)；

简写为：

λ1a+λ2b+lc＝d；

将其写成矩阵形式为：

其解为：

如果λ1,λ2∈[0,1]，则射线和箔条云外包络有交点，当有多个交点时，需要使用最小l法则判断第一交点。

根据计算出的第一交点坐标(x1,y1,z1)，可以求出单个光子发射源点——(xsm,ysm,zsm)与第一交点之间的距离lb，其表达式为：

lb＝[(xsm-x1)²+(ysm-y1)²+(zsm-z1)²]^1/2；

通过lb的求解可以得出光子从源点发射到进入箔条云团之前所经历的时间tb：

c是真空中的光速，其值为3×10⁸m/s。

(5.3)光子在箔条云中依概率传输，使用蒙特卡洛方法求解瞬态矢量辐射传输方程，设n根箔条随机分布在长为ds，体积为dv的柱形体单元中，一辐射强度为i的电磁波通过该单元，由能量守恒定律可以得到该电磁波辐射强度的变化量为：

其中，是入射强度和散射强度之间的过渡系数矩阵，上式的第一项指的是由于箔条吸收κas、散射κs、背景吸收κab导致的能量耗散，第二项表示的是空间其他辐射源辐射的能量，第三项表示的是多散射体耦合散射在方向上的散射强度总和，其中，

如果在空间中，箔条的占空比为fs，则背景的吸收系数表示如下：κab＝2k″·(1-fs)，k″是波数的虚部箔条吸收率是箔条损耗场能和入射场能的比，使用损耗场表示为：

这里引入stokes矢量的四个分量如下：

基于蒙特卡洛方法求解瞬态矢量辐射传输方程包含以下步骤：

1)自由程抽样。通过前面的步骤，可以得到光子第一次进入箔条云的具体位置，并按照下面的自由程公式进行抽样，得到下次碰撞光子前所走路程：

其中，ξ是[0,1]之间的均匀分布随机数，ke是消光系数。对于时域问题，自由行程的距离意味时间的消耗，流逝的时间取决于自由行程的距离以及光在介质中的光速。

传播l距离所耗费的时间为：

其中n为介质的折射率，c0为真空中的光速。则光子在传播l距离后时间变为：

t′＝t0+tb+tl；

2)箔条取向抽样。箔条的取向由空间角确定，设它们服从的概率密度函数为p(θ)和使用舍选法抽取散射方向；

首先，在θ∈[0,π]，之间按照均匀分布抽取箔条取向角即有如下式子：

ξ是[0,1]之间的均匀随机数；

将箔条取向角带入概率密度函数

其次，在区间上使用上述的ξ再次产生一个均匀分布随机数p0，按照下面的判断准则进行发射角选择：

如果上述准则给出的是“拒绝选择”则需要返回第一步继续进行，直到产生满足要求的数对为止。

3)散射方向抽样。基于能量守恒和舍选法进行散射方向抽样，设在某一确定的入射角下，归一化的散射能量可以表示为：

且：

则粒子沿着方向出射的概率密度为

其中

其中，方向散射系数定义如下：

散射方向由空间角确定，的选取采用舍选法进行，和上文中选择入射方向方法一致，这里不再赘述。

4)散射传输能量获取。即计算每次碰撞的衰减因子在给定的入射波条件下，有：

这样单次碰撞前的能量wi-1和碰撞后的能量wi之间关系为：

重复上述步骤，直到光子射出箔条云，遍历设定的最大光子数便可以进行后续的数据处理，判断光子是否在云内外的基本判定方法是，从该点向任意方向发射射线，判断射线与几何体的交点个数，奇数为云内，偶数则为云外。

单个光子在箔条云中经历的时间是通过自由程的叠加得到的，光子每一次碰撞，就会经历一个自由程l，直至最终光子从箔条云中逃逸或光子携带能量在碰撞过程损耗殆尽。单个光子经历的所有自由程叠加的结果与光在箔条云中的速度的比值就得到了该光子在箔条云中经历的时间tl。

(5.4)光子出射与能量收集。连续地按照上述的过程产生光子、光子与散射体碰撞(能量损失)并追踪光子的轨迹与光子经历这些轨迹所用的时间，当光子离开散射体后本发明将对其散射能量与时间进行统计。

设有n个能量为i0i的光子被发射，有m个光子经连续碰撞逃出云团，散射方向分别为自身携带能量分别为ii，此时光子经历的时间为t′i，对全空间的对进行网格划分，θ角的步长为δθ，角的步长为同时有θ＝mδθ、判断出射光子是否进入第(m,n)单元的方法是：

于是，在方向接收到的总光子能量为：

而第i个光子遍历整个过程所经历的时间则由t′i表示：

t′i＝t0+tb+tl；

其中t0表示光子在源平面出射时的初始时间，tb表示光子从出射到进入箔条云之前所经历的时间，tl表示光子进入箔条云，在其中折射和反射后经历自由程，最终出射云边界所用的时间。

(6)对求解数据进行运算与处理，得到所需要的时域stokes矢量结果，见于上述stokes矢量的定义及入射光子能量i0i的表示形式，在方向的雷达散射截面可以表示为：

上述结果按照stokes矢量四个分量形式表示。鉴于此时方向的时域散射结果已通过运算得到，而方向收集的光子携带的时间信息也通过t′i反映出来，通过光子的索引信息i可将t′i与σi一一对应，求得σi随t′i变化的波形即为所求的箔条云的时域stokes矢量结果。

下面结合仿真对本发明的技术效果作详细的描述。

一、仿真条件：仿真采用不同光子数入射箔条云，参数如表1所示，散射计算坐标系如图9所示，一种可能的天线方向图形式如图10所示，天线与箔条云的位置关系如图11所示。

表1实验条件一览表

二、仿真内容与结果

仿真1，采用表1中1号实验参数，结果如图12～图15所示，这里展示stokes矢量四个分量i、q、u、v随时间变化波形。

仿真2，采用表1中2号实验参数，结果如图16～图19所示，这里展示stokes矢量四个分量i、q、u、v随时间变化波形。

仿真3，采用表1中3号实验参数，结果如图20～图23所示，这里展示stokes矢量四个分量i、q、u、v随时间变化波形。

仿真4，采用表1中4号实验参数，结果如图24～图27所示，这里展示stokes矢量四个分量i、q、u、v随时间变化波形。

下面结合实验对本发明的应用效果详细的描述。

通过一系列的实验仿真对系统性能进行了评估，对大体量箔条云散射进行计算，在相同的条件下，鲜见公开报道的对比结果，与实测结果对比，上述四种仿真的误差均在15％以内，由于行业公开箔条云散射计算结果，因此本发明的适用性和准确度都较高，对于简单类型的大体量箔条云本发明效果更好。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。