基于频率角度二维混合插值的降雨粒子散射运算加速方法与流程

本发明属于雷达技术领域，尤其涉及一种基于频率角度二维混合插值的降雨粒子散射运算加速方法。

背景技术：

雷达的概念起源于20世纪初，因具有全天候、全天时远距离遥感测能力而成为军事领域不可或缺的电子设备。随着科技水平的飞速发展，针对不同应用背景的雷达体制也广泛出现在了国民经济的各个领域。降雨雷达作为天气雷达的一个重要分支，在对降水、冰雹等天气现象进行监测和预报方面起着重要作用，得到了各国政府的高度重视。从降雨的电磁极化特性角度考虑，降雨粒子对雷达波束的散射是降雨探测的基础。但是，降水粒子的散射特性很复杂，不仅与降水粒子的尺寸、形状及相态有关，还在很大程度上依赖于粒子在空间的方位、取向等因素。因此，从理论上分析降雨粒子的散射特性具有重要的研究价值。

近年来，在有关环境的遥感应用中，对具有复杂的降雨粒子群电磁波相互作用的表征研究是非常重要的。常见的解决复杂的降雨粒子群散射的方法是依靠数字技术，例如，矩量法(MOM)。矩量法(MOM)可以提供复杂介质散射的精确解决方案，但是在求解积分方程时运算非常复杂并且需要很大的存储空间。为了减轻这些瓶颈，人们提出了一些快速、有效的方法，如：共轭梯度快速傅里叶变换(CG-FFT)法、快速多极子法(FMM)等来降低计算机内存和计算量的需求，提高矩量法的效率。曾亮和陈志豪教授等在分析随机粗糙面散射问题时，将阻抗矩阵分为强相互作用矩阵和弱相互作用矩阵与待求向量的乘积，大大降低了计算复杂度和存储复杂度。近年来，一个快速数值方法——稀疏矩阵/规则网格(SM/CG)法被提出，应用于三维降雨粒子群的散射分析中。

雷达散射截面(RCS)同时包含频率和角度信息。在许多实际应用中，理想中预测的单站RCS目标同时存在频率域和空间域。虽然可以通过SM/CG方法来降低计算复杂度和内存的要求，但我们还是需要重复的在感兴趣频段的各频率或角度计算，以便能够在很宽的角频域波段获得RCS。为了解决这一困难，提出了许多内插方法提高效率，如渐近波形估计(AWE)法和基于模型的参数估计(MBPE)法。但是这些方法应用起来也很困难，AWE方法使用时，有时难以得到的阻抗矩阵和感应电流矢量的导数。MBPE方法的困难之处在于需要矩阵求逆得到替代模型的系数，导致计算频率和角度扫描时浪费了大量时间。

技术实现要素：

针对现有技术存在的问题，本发明提供一种基于频率角度二维混合插值的降雨粒子散射运算加速方法，在计算频率和角度扫描时节省了大量时间并且有效的提高了降雨粒子群带宽散射的运算速度。

本发明针对上述技术问题主要是通过下述技术方案得以解决的：

基于频率角度二维混合插值的降雨粒子散射运算加速方法，包括如下步骤：

(1)在角度域应用三次样条插值法求得除了两个边界点以外的所有采样点的一阶导数；

(2)在宽频带域应用Stoer-Bulirsch算法，通过计算采样点得到有理插值函数；

(3)采用二维混合自适应插值方法对目标波段的角度域和频率域进行二维自适应采样，并对采样结果进行由粗到细的结构划分。

进一步的，所述步骤(1)在角度域应用三次样条插值法求得除了两个边界点以外的所有采样点的一阶导数具体包括：

假设在角度域中采样节点将整个角度范围划分为多个区间；在每个区间中是一个在处展开的三阶多项式，通过Hermite插值理论，的插值公式为：

其中表示采样点，为了计算角度范围中的值，需要得到和和的一阶导数可以在零自然边界条件下求出，然后通过式(2)得到其他采样点的一阶导数；

其中h_i表示和之间的间隔，

进一步的，所述步骤(2)在宽频带域应用Stoer-Bulirsch算法，通过计算采样点得到有理插值函数具体为：

2a)用已知的采样数据为表示有理函数，并且已知的采样数据为(f_i,P(f_i))i＝1,...,k；其中，A(f)，B(f)分别是次数不超过p和q的多项式，P(f_i)为采样点f_i对应的采样值；

2b)令R_i,1＝P(f_i)i＝1,...,k，R(f)＝R_1,k，应用Stoer-Bulirsch算法，根据式(3)-(5)递归求解出任介有理函数在任意位置f插值R(f)：

R_j,k＝((f-f_j)R_j+1,k-1+(f_j+k-f)R_j,k-1)/(f_j+k-f_j) (3)

R_j,k＝(f_j+k-f_j)/((f-f_j)/R_j+1,k-1+(f_j+k-f)/R_j,k-1) (4)

R_j,k＝R_j+1,k-2+(f_j+k-f_j)/((f-f_j)/(R_j+1,k-1-R_j+1,k-2)+(f_j+k-f)/(R_j,k-1-R_j+1,k-2)) (5)。

5、进一步的，所述步骤(3)采用二维混合自适应插值方法对目标波段的角度域和频率域进行二维自适应采样，并对采样结果进行由粗到细的结构划分具体为：

3a)取四个采样点(f₁,φ₁),(f₁,φ_n),(f_m,φ₁)和(f_m,φ_n)做线性插值，然后得到I(f₁,φ₁),I(f₁,φ_n),I(f_m,φ₁)和I(f_m,φ_n)；

3b)假设在一个平面上有m×n个采样点，行为m，列为n，通过直线插值得到近似雷达横截面rcs_a1(f,φ)，通过三次样条插值近似得到雷达横截面rcs_a2(f,φ)；

3c)在每个区间选取这两条误差曲线中平均误差最大的点的位置也就是φ_tj,得到采样线；

3d)：当收敛误差CE(φ)，大于步骤3c)中的最大误差ME，此时采样点最好，否则将插值函数的曲线添加到采样区间中；

3e)在每个区间[f_si,f_si+1](i＝1,2,…,m),当ME出现时通过Stoer-Bulirsch算法得到采样线；

3f)当出现CE(f)大于步骤3e)中的ME时，保留此时对应的采样点，否则将插值函数的曲线添加到采样区间中；

第3g)如果每个区间满足步骤3d)中的条件，并且在[f_si,f_si+1]区间满足步骤3f)中出现的条件，则采样结束。

本发明的有益效果是：

1.应用Stoer-Bulirsch算法计算采样点得到有理插值函数，并且不需要求解矩阵的逆，所以它能通过许多采样点得到有理插值函数，同时避免了奇异性问题，大幅节省了计算机的运算时间。

2.三次样条插值法通过另一种方式来获得每个采样节点的一阶导数而不是反复的求解线性方程组，在角度域扫描时大大降低了扫描成本。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是基于频率角度二维混合插值的降雨粒子散射运算加速方法的结构示意图；

图2是雨滴粒子的形状近似图；

图3是不同大小降雨粒子的横截面；

图4(a)是模型I的示意图；

图4(b)是模型I在垂直极化下的RCS曲线；

图5(a)是模型II的示意图；

图5(b)是模型II在垂直极化下的RCS曲线；

图6(a)是模型III的示意图；

图6(b)是模型III在垂直极化下的RCS曲线。

具体实施方式

下面结合附图对本发明的优选实施例进行详细阐述，以使本发明的优点和特征能更易于被本领域技术人员理解，从而对本发明的保护范围做出更为清楚明确的界定。

本发明的主要思路为：在感兴趣波段(即目标波段)的感应电流角度域应用三次样条算法；在宽频带降雨粒子群电磁特性分析中应用Stoer-Bulirsch算法；通过二维混合自适应插值技术，采用由粗略到精细的分层结构作为一个反复的过程，以生成一组不均匀的采样节点模型，来提高运算效率。

基于频率角度二维混合插值的降雨粒子散射运算加速方法，包括：

1、在角度域应用三次样条插值法，不需要反复求解线性方程组来获得每个采样点的导数，可直接求得除了两个边界点以外的所有采样点的一阶导数；具体为：假设在角度域中采样节点将整个角度范围划分为多个区间；在每个区间中是一个在处展开的三阶多项式，通过Hermite插值理论，的插值公式为：

其中表示采样点，为了计算角度范围中的值，需要得到和和的一阶导数可以在零自然边界条件下求出，然后通过式(2)得到其他采样点的一阶导数；

其中h_i表示和之间的间隔，

2、在宽频带域应用Stoer-Bulirsch算法，通过计算采样点得到有理插值函数；具体为：

首先，用已知的采样数据为表示有理函数，并且已知的采样数据为(f_i,P(f_i))i＝1,...,k；其中，A(f)，B(f)分别是次数不超过p和q的多项式，P(f_i)为采样点f_i对应的采样值；

然后，令R_i,1＝P(f_i)i＝1,...,k，R(f)＝R_1,k，应用Stoer-Bulirsch算法，根据式(3)-(5)递归求解出任介有理函数在任意位置f插值R(f)；该方法不需要进行矩阵求逆运算，从而在根本上避免了奇异性问题。

R_j,k＝((f-f_j)R_j+1,k-1+(f_j+k-f)R_j,k-1)/(f_j+k-f_j) (3)

R_j,k＝(f_j+k-f_j)/((f-f_j)/R_j+1,k-1+(f_j+k-f)/R_j,k-1) (4)

R_j,k＝R_j+1,k-2+(f_j+k-f_j)/((f-f_j)/(R_j+1,k-1-R_j+1,k-2)+(f_j+k-f)/(R_j,k-1-R_j+1,k-2)) (5)。

其中，应用Stoer-Bulirsch算法的具体计算流程如下：

第1步:计算宽频带两端f1、f2的真实值EM(f₁)、EM(f₂)，并将f₁、EM(f₁)、f₂、EM(f₂)加入到采样组中

第2步:利用式(3)和式(4)作为递归公式，分别计算宽频带内各感兴趣频率点f_t的值P₁(f_t)和P₂(f_t)，并找出最大差异点f₃及差异值error12。

$error12=|P_1\left(f_3\right)-P_2\left(f_3\right)|=\underset{i}{max}|P_1\left(f_{ti}\right)-P_2\left(f_{ti}\right)|---\left(6\right)$

第3步:如果最大值error12比设定的误差门限eps小，则结束计算，转到第6步。如果error12比eps大，则计算f₃处的真实值EM(f₃)，并将f₃和EM(f₃)加入到采样组中。

第4步:利用式(3)和式(5)作为递归公式，分别计算宽频带内各感兴趣频率点f_t的值P₁(f_t)和P₃(f_t)，并找出最大差异点f₄及差异值error13。

$error13=|P_1\left(f_4\right)-P_3\left(f_4\right)|=\underset{i}{max}|P_1\left(f_{ti}\right)-P_3\left(f_{ti}\right)|---\left(7\right)$

第5步:如果最大值error13比设定的误差门限eps小，则结束计算，转到第6步。如果error13比eps大，则计算f₄处的真实值EM(f₄)，并将f₄和EM(f₄)加入到采样组中，然后转入第2步进行下一个采样点的计算。

第6步:根据采样中已有的数据用式(3)建模，对各感兴趣频率点进行插值计算。

3、采用二维混合自适应插值方法对目标波段的角度域和频率域进行二维自适应采样，并对采样结果进行由粗到细的结构划分；具体为：

第1步：取四个采样点(f₁,φ₁),(f₁,φ_n),(f_m,φ₁)和(f_m,φ_n)做线性插值。然后得到I(f₁,φ₁),I(f₁,φ_n),I(f_m,φ₁)和I(f_m,φ_n)。

第2步：假设在一个平面上有m×n个采样点，行为m，列为n，通过直线插值得到近似雷达横截面rcs_a1(f,φ)，通过三次样条插值近似得到雷达横截面rcs_a2(f,φ)。

第3步：在每个区间选取这两条误差曲线中平均误差最大的点的位置也就是φ_tj,得到采样线，比如

第4步：当收敛误差CE(φ)，大于步骤3中的最大误差ME，此时采样点最好，否则将插值函数的曲线添加到采样区间中。

第5步：在每个区间[f_si,f_si+1](i＝1,2,…,m),当ME出现时通过Stoer-Bulirsch算法得到采样线。

第6步：当出现CE(f)大于步骤5中的ME时，此时采样点最好，否则将插值函数的曲线添加到采样区间中。

第7步：如果每个区间满足步骤4中的条件，并且在[f_si,f_si+1]区间满足步骤6中出现的条件，此时自适应采样过程就结束了。

下面给出了几个数值算例对本方法加以说明。

为了验证S-B(Stoer-Bulirsch)自适应采样方法在计算宽频带降雨粒子群上的有效性，本节给出了几个数值算例加以说明。求解矩阵方程的迭代算法采用的是GPBi-CG，迭代的初始值均设为零，最大迭代步数设置为10000，迭代精度为10^-5。本节中所有算例计算的工作环境均为Intel Core(TM)II Quad处理器，3.50GB内存，2.83GHz主频，单精度。

模型I.1000个位置随机大小任意的球形粒子群

在0.3*0.3*0.3m³空间内随机分布着1000个半径0-2mm的球形降雨粒子，未知量个数为35308。雨滴粒子的介电常数随着频率变化，介于4+0.04i和3.5245+0.08755i之间。

模型II.500个位置随机大小任意的椭球形粒子群

在实际的降雨中，由于受到重力和空气阻力的双重影响，大多数的降雨粒子形状并非是标准球形，而是呈现为类似椭球的结构，如图2所示。其中，椭球的短半轴半轴长为a₃，长半轴半轴长分别为a₁,a₂(a₁＝a₂)。于是，降雨粒子的形状就可以通过椭球的轴比(a₁/a₃)来描述。

本算例设计了在0.3*0.3*0.3m³空间内随机分布着500个半短轴0-2mm的椭球形降雨粒子，椭球轴比在1-2之间随机变化，未知量个数为34604。雨滴粒子的介电常数随着频率变化，介于4+0.04i和3.4+0.1i之间。在1-13GHz的频率范围内每0.01GHz取一个采样点，用SM/CG法逐点扫描需计算1201个频率点。

模型III.500个位置随机大小满足Keenan分布的降雨粒子群

上一个例子描述了降雨粒子为椭球体时的情形，而在真实降雨中，降雨粒子的形状与大小有着密切的关系，图3给出了不同大小的雨滴与对应等效球体的横截面示意图。由图可以看出，随着雨滴体积的增大，真实的雨滴形状呈现出越来越接近扁椭球的特点。

本算例中，我们假设在0.5*0.5*0.5m³空间内随机分布着500个半短轴0-2mm的椭球形降雨粒子粒子，它们的轴比关系满足Keenan et al.模型，未知量个数为27498。雨滴粒子的介电常数随着频率变化，介于4+0.04i和3.3495+0.105i之间。在10-23GHz的频率范围内每0.01GHz取一个采样点，用SM/CG法逐点扫描需计算1301个频率点。

图4、5和6分别给出了上述三个例子直接采用SM/CG方法和自适应S-B插值方法的计算结果。从图中可以看出，三种情况下的计算结果均吻合良好。算例1中，如图4所示，在1-20GHz的频率范围内以20MHz为采样间隔，当采用SM/CG法逐点扫描时需计算951个频率点，而当结合自适应S-B插值方法时，直接计算的采样点个数降为224。算例2、3中，如图5、6所示，频率采样间隔为10MHz，当采用SM/CG法逐点扫描时分别需计算1201和1301个频率点，而当结合自适应S-B插值方法时，直接计算的采样点个数分别降为387和310。

表1给出了以上三个算例分别直接采用SM/CG方法和自适应S-B插值方法在计算时间上的比较。从中可以看出，在使用自适应S-B插值方法之后，三种情况下，计算总频带RCS值所需的时间分别只有原来的23.66％，26.62％和28.20％。值得注意的是，该采样算法具有自适应特点，因此可以根据实际情况，自动选择较为适合的采样频率点位置和采样点个数，即RCS曲线较复杂，波动较剧烈时采样点数较多，而RCS曲线较简单，变化较平缓时采样点数较少。以上算例均表明，自适应S-B插值方法是一种计算降雨粒子宽频带电磁特性的有效方法。

以上所述，一种基于三次样条结合Stoer-Bulirsch模型二维自适应混合插值技术应用在降水粒子宽频带电磁特性的快速计算。相比传统的AWE技术，这种新方法无须矩阵求逆，避免一阶导数的计算，计算结果表明，自适应S-B插值方法是一种计算降雨粒子宽频带电磁特性的有效方法。

表1直接采用SM/CG方法计算与使用自适应S-B采样方法计算的对比

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何不经过创造性劳动想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书所限定的保护范围为准。