基于T矩阵法的机载极化气象雷达降水粒子回波仿真方法

基于t矩阵法的机载极化气象雷达降水粒子回波仿真方法
技术领域
1.本发明属于气象雷达信号处理技术领域，具体涉及基于t矩阵法的机载极化气象雷达降水粒子回波仿真方法。


背景技术：

2.球变暖问题的不断加剧导致世界范围内暴雨、暴雪、冰雹等强对流天气频发，使得飞机遭遇危险天气的概率大大增加，对民航飞机安全稳定的运行造成严重影响。现代航空利用机载气象雷达设备实时监测飞行航路上的危险气象并告警，辅助机组做出正确的飞行决策，确保飞行顺利稳定。现有机载气象雷达采用单极化天线进行气象探测，在遭遇对流云团等复杂气象天气情况下，无法提供更精细化的危险区域降水粒子检测与分类信息。作为即时监测飞机航线上的气象环境的重要航空电子设备，机载气象雷达亟需更进一步提升其探测的准确性，为安全飞行提供更可靠的气象监测信息，因此持续开展极化气象雷达相关技术研究对保障民航安全非常必要。
3.2012年，作为民用飞机机载气象雷达行业非常重要的设备制造商，美国柯林斯公司(rockwell collins)申请了机载气象雷达双极化天线专利，表明了机载气象雷达将借助双极化技术来推动机载气象雷达进一步发展。但在机载双极化气象雷达相关技术起步研究阶段，获取实测数据成本代价大、困难程度高，仿真带有极化信息的雷达回波数据有助于相关双极化气象雷达检测与算法技术的深入研究。


技术实现要素：

4.针对现有技术的存在的问题，本发明旨在提供基于t矩阵法的机载极化气象雷达降水粒子回波仿真方法。
5.为了实现上述目的，本发明采用如下技术方案，
6.基于t矩阵法的机载极化气象雷达降水粒子回波仿真方法，所述方法包括以下步骤：
7.s1利用数值预报模式对气象数据进行准确建模仿真；
8.s2利用降水粒子微物理特性和t矩阵法求解不同降水粒子反射率因子；
9.s3根据雷达基本参数和雷达气象方程实现机载极化气象雷达降水粒子回波仿真。
10.需要说明的是，在步骤s1中，利用数值预报模式对气象数据进行准确建模仿真方法是：将全球地形数据、初始驱动场和模拟区域的方案参数输入至wrf模式，利用地图插值和数值积分方法获取六种降水粒子的混合比、粒子数浓度、风速和空间坐标的参数。
11.需要说明的是，在步骤s2中，利用降水粒子微物理特性和t矩阵法求解不同降水粒子反射率因子方法是：利用t矩阵法结合不同降水粒子形状、尺寸、复折射率的参数计算单一粒子散射矩阵，再通过空间分辨网格内降水粒子微物理特性求解反射率因子。
12.需要说明的是，在步骤s3中，根据雷达基本参数和雷达气象方程实现机载极化气象雷达降水粒子回波仿真方法是：运用雷达气象方程，设定雷达参数，结合降水粒子空间坐
标、反射率因子等参数，实现机载双极化雷达回波i/q回波数据仿真。
13.本发明的有益效果在于，可准确反映气象特征，与实测数据的对比分析证实了所提方法的有效性、可靠性。
附图说明
14.图1为本发明提供的基于wrf的机载极化气象雷达回波仿真原理框图；
15.图2为wrf算法流程图；
16.图3为本发明提供的wrf嵌套结构示意图；
17.图4为本发明提供的基于wrf模式的气象场景仿真结果；
18.图5为散射坐标系示意图；
19.图6为本发明提供的仿真场景示意图；
20.图7为本发明提供的回波仿真结果与地基雷达反射率因子对比图；
21.图8为本发明提供的六种降水粒子反射率因子与差分反射率因子。
具体实施方式
22.以下将对本发明作进一步的描述，需要说明的是，以下实施例以本技术方案为前提，给出了详细的实施方式和具体的操作过程，但本发明的保护范围并不限于本实施例。
23.本发明为基于t矩阵法的机载极化气象雷达降水粒子回波仿真方法，所述方法包括以下步骤：
24.s1利用数值预报模式对气象数据进行准确建模仿真；
25.s2利用降水粒子微物理特性和t矩阵法求解不同降水粒子反射率因子；
26.s3根据雷达基本参数和雷达气象方程实现机载极化气象雷达降水粒子回波仿真。
27.需要说明的是，在步骤s1中，利用数值预报模式对气象数据进行准确建模仿真方法是：将全球地形数据、初始驱动场和模拟区域的方案参数输入至wrf模式，利用地图插值和数值积分方法获取六种降水粒子的混合比、粒子数浓度、风速和空间坐标的参数。
28.需要说明的是，在步骤s2中，利用降水粒子微物理特性和t矩阵法求解不同降水粒子反射率因子方法是：利用t矩阵法结合不同降水粒子形状、尺寸、复折射率的参数计算单一粒子散射矩阵，再通过空间分辨网格内降水粒子微物理特性求解反射率因子。
29.需要说明的是，在步骤s3中，根据雷达基本参数和雷达气象方程实现机载极化气象雷达降水粒子回波仿真方法是：运用雷达气象方程，设定雷达参数，结合降水粒子空间坐标、反射率因子等参数，实现机载双极化雷达回波i/q回波数据仿真。
30.实施例
31.如图1所示，本发明提供的基于t矩阵的机载极化气象雷达降水粒子回波仿真方法包括按顺序进行的下列步骤：
32.利用数值预报模式对气象数据进行准确建模仿真
33.wrf模式是美国国家大气研究中心、国家环境预报中心和美国环境模拟中心以及众多高等科研院校、气象勘测单位、气象业务应用部门等于1997年联合研究开发的数值预报模式业务应用产品，是全球中尺度预报模式中应用最广泛的天气预报模式，wrf主要考虑10公里以下的网格分辨率，利用数值方法和资料同化技术实现气象数据模拟，是从云尺度
到降水尺度等不同尺度的重要天气预报工具。本文选用arw4.1版本，模式的物理过程包含了大气的水平涡动扩散、垂直涡动扩散，积云对流参数化方案、云物理方案，太阳短波辐射和大气长波辐射方案等。wrf模式的基本配置方案，如下表所示。
34.运算方式计算方案动力框架完全可压缩非静力平衡动力框架控制方程组通量形式水平方向arakawa c网格垂直方向地形追随质量坐标系时间积分四阶的runge-kutta算法
35.基于wrf的气象场景仿真流程图如下图2所示，将dem地形数据与温度、风速、压强等初始场数据作为模式运行输入数据，设置投影方式与嵌套层数将输入数据插值到模拟区域网格中，结合积分步长、模式物理参数方案、积云对流参数等方案设置，迭代模拟输出仿真区域空间坐标、不同降水粒子的粒子数浓度、混合比等数据，实现气象场景仿真。
36.本文以2015年6月17日热带风暴“比尔”为气象场景仿真典型实例。该降水场景于17日12时登陆美国得州马特高尔达岛，该热带风暴最高风速约20米每秒，携带大量水汽，在得州造成局部强降雨，降雨量最高达250毫米，满足对气象场景仿真的气象要求。gmt时间2015年6月17日15:00左右，该风暴经过得州东部，被wsr-88d雷达(编号：kfws，北纬32.34
°
，西经97.18
°
，海拔高度231米)监测。为高保真实现降水粒子场景建模仿真，仿真时长范围至少在24小时，因此模式仿真时间段从2015年6月17日12:00至18日12:00。
37.仿真场景中心经纬度为北纬32.34
°
，西经97.18
°
，地图投影选用兰伯特投影方式。每一层选用合适的微物理方案，保证能输出六种降水粒子的气象参数，准确仿真气象场景。由于不同嵌套网格气象场景精度的不同，本文对于仿真结果验证均采用子网格的气象模拟数据作为研究对象。模式运行嵌套示意图如图3所示。图中外边框表示母网格模拟区域，带标记符号的白色方框表示第二层网格模拟区域，红色矩形框则为最内层模拟区域。
38.图4表示基于wrf模式的降水粒子仿真可视化结果。图4(a)表示雨滴、云滴、冰雹、霰、冰晶和雪晶等六种降水粒子水汽混合比，图4(b)表示雨滴、云滴、冰雹、霰、冰晶和雪晶等六种降水粒子数浓度。如图4(a)所示，wrf输出的水汽混合比，用于描述仿真场景中空气的含水量，单位为g/kg。从图中可以看出，在该时刻热带风暴在中心处产生涡旋，在涡旋中心及边缘区域降水含量较大。如图4(b)所示，该图表示降雨粒子浓度分布，在涡旋附近，降雨粒子浓度较大。对比图4(a)和图4(b)相同区域，区域红线部分处于热带风暴漩涡边缘，水汽混合比相对中心区域略有增加，并粒子数浓度也有所增大，表明风暴运动带来了大量降水，导致降水粒子数据增大，符合风暴气象特征。通过对比图4(a)和图4(b)可以初步判断wrf仿真热带风暴基本符合热带风暴特性，可用于后续极化雷达回波仿真计算与分析。
39.2)利用降水粒子微物理特性和t矩阵法求解不同降水粒子反射率因子
40.为了给降水粒子回波仿真提供不同极化通道的反射率因子数据，分析雨滴、云滴、冰雹、霰、冰晶和雪晶等六种不同降水粒子类型极化散射特性，计算对应散射矩阵求解反射率因子是回波仿真的关键步骤。
41.首先对六种不同类型球形降水粒子的物理特性做统计归纳，用t矩阵法求解散射矩阵时需要考虑六种降水粒子的尺寸、半径、离子积常数、粒子数密度等物理性质。由于降
水粒子的尺寸、形状、复折射系数等都会影响降水粒子的散射能力，进而导致其回波信号的强弱。
42.对六种降水粒子的物理特征进行统计描述，如下表所示。
[0043][0044]
其次应用t矩阵法计算不同降水粒子类型单个粒子散射矩阵，以雨滴粒子散射矩阵计算为例云滴、雪晶、霰、冰雹和霰粒子散射矩阵计算类似，构建散射坐标系l，将坐标原点位于粒子内部，如图5所示，其方向在空间固定，其原点位于粒子内部。平面入射电磁波的传播方向由单位矢量或等效地由一对指定，其中是从z轴正方向看向坐标原点，入射波与坐标轴的夹角，θ是方位角，即当沿z轴正方向看时，从x轴正方向沿顺时针方向旋转角度。电场的分量和φ分量分别表示为和分量位于穿过电磁波和z轴的平面中，而分量垂直于该平面，其中和是相应的单位向量，如下关系：
[0045][0046]
和通常也分别表示为e
ν
和eh，分别称为垂直和水平电场矢量分量。
[0047]
考虑具有电场矢量的单极化平面入射电磁波表达如下：
[0048][0049]
上式中入射波以方向入射到粒子上，这里的k＝2π/λ是自由空间波数，λ是自由空间波长，r是连接标准坐标系原点和观测点的位置矢量，下标l表示标准坐标系中计算的向量分量。由于麦克斯韦方程组和边界条件的线性关系，始终可以用入射电场线性表示散射电场。在远场区域(kr＞＞1,r＝|r|)中，散射波可以用矢量球面波函数表示，由下式给出：
[0050][0051][0052]
其中s
l
是2
×
2的散射矩阵在这里求解散射矩阵选用t矩阵法来计算散射矩阵元素值。
[0053]
t矩阵方法优势在于可以不考虑入射场与散射场，而只从微粒本身几何结构(形状、尺寸)和自然属性(复折射率)、粒子在所选坐标系中的位置出发，以此计算粒子的散射
矩阵在计算时把散射场的矢量球谐函数的展开系数用入射场的矢量球谐函数的展开系数表示，转换矩阵即t矩阵，当降水粒子的大小、形状及复折射系数确定，则粒子的t矩阵也就确定，而与入射场、散射场无关。文献[19]-[21]指出，当粒子形状趋于极端形状以及大尺度情况下，t矩阵法计算误差增大，但在适用范围内计算精度高、速度快。入射场和散射场用矢量球面波函数m
mn
和n
mn
描述，表示如下：
[0054][0055][0056]
式(6)中r0表示粒子外接球半径，在粒子坐标系下，坐标系原点位于粒子的内部。式(5)和式(6)中的展开系数可以表达如下：
[0057][0058][0059][0060][0061]
由于入射场和散射场满足麦克斯韦方程，则粒子的入射场、散射场展开系数是线性相关的，采用矩阵形式来描述，上述方程可写为：
[0062][0063]
其中t＝-q
(1)
(q
(3)
)-1
,q的上角标(1),(3)分别表示第一类和第三类bessel函数。
[0064]
电磁波的入射场展开系数可通过近似值解析表达式获得，当散射粒子的t矩阵
[0065]
被确定时，相应的散射场就可以利用式(4)-(6)求得，具体表达式如下：
[0066][0067]
[0068][0069]
上式中依次为散射电磁波和入射电磁波在坐标系中的和m和n表示散射求和点数，不同粒子的等体积球体直径d、形状因子μ、波长λ和复折射率m用于确定单个粒子散射矩阵求解时的求和上限nmax，计算公式如下式(16)，其他系数α
mnn'
、和计算如下：
[0070]
nmax＝ψ(μ
,
λ,d,m)
ꢀꢀꢀꢀ
(16)
[0071][0072][0073][0074][0075]
再结合代入公式(12-15)，计算出不同降水类型单个粒子的散射矩阵的各个元素值。
[0076]
本专利以气象场景建模仿真结果为研究对象，考虑降水粒子形状为球形、入射波长3.2cm条件下，利用t矩阵法计算了雨滴(等体积球体半径＝10μm，复折射率＝66+33j)、冰雹(等体积球体半径＝5μm，复折射率＝3.17+5.4
×
10-3j)、霰(等体积球体半径＝10μm，复折射率＝1.97+0.271j)、雪晶(等体积球体半径＝5μm，复折射率＝6.05+1.57j)、冰晶(等体积球体半径＝4μm，复折射率＝3.17+5.4
×
10-3j)、云滴(等体积球体半径＝2μm，复折射率＝66+33j)等六种降水粒子物理特征参数下对应的散射矩阵。
[0077]
在进行雨滴、云滴、冰雹、霰、冰雹和雪晶降水粒子的雷达回波仿真计算时，需要引入变量表征空间不同类型群粒子目标对电磁波的反射能力，即反射率因子。雨滴、云滴、冰雹、霰、冰雹和雪晶反射率因子反映了一个空间分辨单元内所有同类降水粒子散射能力的叠加结果，因此在分析时需要结合其微物理特性，如倾角分布、雨滴谱分布和轴比分布，当空间分辨单元尺度远远小于气象场景尺度时，可以认为空间分辨单元内粒子是均匀分布的，因此可以用粒子尺寸d、方位角θ、俯仰角三维积分来计算该空间分辨单元的反射率因子。以雨滴为例，雨滴的反射率因子z
ij
的具体计算公式如下：
[0078][0079]
式中，h和v分别表示雷达电磁波的极化方向，h表示水平方向，v表示垂直方向，λ表示入射电磁波波长，d
rain
表示雨滴的等体积球体直径，具体数值参考表2，|kw|2为雨滴的离子积常数值取0.93。σ
ij
(d
rain
)为雨滴粒子的后向散射截面积，可利用式(12-15)求得的散射
矩阵求得。p(θ),分别表示方位角θ和俯仰角的概率密度函数，通常取高斯型。n(d
rain
)为雨滴粒子雨滴谱函数，通常选用gamma分布，用于表示粒子直径与粒子数之间的关系，利用wrf仿真的不同降水粒子的气象场景结果，如降水粒子混合比和粒子数浓度数据。如上即分析了雨滴粒子的反射率因子计算过程，其他粒子(霰、冰雹、冰晶、雪晶和云滴)的反射率因子计算公式同上所示，离子积常数取值参考表2如：霰、冰雹、冰雹和雪晶取0.2；云滴取0.93。
[0080]
3)根据雷达基本参数和雷达气象方程实现机载极化气象雷达降水粒子回波仿真
[0081]
当气象雷达进行雨滴、云滴、冰雹、霰、冰雹和雪晶降水粒子回波仿真时，仿真场景坐标系示意图如图6所示。建立标准坐标系(x,y,z)，其原点(x0,y0,z0)平行于地球表面。雷达坐标系(x',y',z')以雷达位置(xa,ya,za)为中心，其中雷达天线相对于地球表面的仰角为θ。第p个降水粒子位置的坐标为
[0082]
根据降水粒子雷达气象方程实现空间群粒子极化雷达回波信号的仿真，即利用反射率因子来计算不同极化方向上接收降水粒子的平均电磁波功率第p个降水粒子的表示如下：
[0083][0084]
式中，i和j分别表示雷达电磁波的极化方向，h表示水平方向，v表示垂直方向，r表示雷达与目标的径向距离，z
ij
表示雨滴、云滴、冰雹、霰、冰雹和雪晶不同类型降水粒子不同极化方向下的反射率因子。根据以上公式即可计算出每个极化通道距离r处从雷达分辨单元接收的单个脉冲的平均功率，最后结合雨滴、云滴、冰雹、霰、冰雹和雪晶粒子的相位信息，则距离r处雷达分辨单元所仿真的单个脉冲雷达回波信号可以表示为：
[0085][0086][0087]
式(23)中wd表示时间角频率；f0是载波频率；r表示雷达与目标径向距离；wr表示时延相位；z
ij
表示反射率因子。c为常数，具体计算如式(24)，p
t
表示雷达发射功率，g0表示天线增益，τ表示脉冲宽度，c表示光速。
[0088]
为验证本方法效果，本发明人选取美国俄克拉荷马州北纬32.34
°
，西经97.18
°
为中心经纬度，2015年6月17日15:00的气象场景进行仿真。雷达参数如下表所示。
[0089]
参数参数值参数参数值飞机高度(m)6000波束宽度(
°
)3.5工作频率(ghz)9.375脉冲重复频率(hz)800飞机速度(m/s)200采样脉冲数64距离分辨率(m)150脉冲宽度(μs)1扫描方式ppi天线扫描速度(
°
/s)45
[0090]
图7是混合回波仿真结果和地基雷达反射率因子对比图。因为极化参数是目前气象工作者通过极化雷达回波对气象场景反演的主要参数。因此在仿真结果对比分析时利用反射率因子对其仿真结果进行展示。图7(a)是雨滴、云滴、雪晶、霰、冰雹和霰粒子混合回波仿真结果，图7(b)是kfws地基雷达反射率因子图，从图7(a)(b)两图可以看出根据雷达回波反演的水平反射率因子(左图)与地基雷达kfws实测反射率因子(右图)分布特征基本一致，符合风暴气象变化特征。对比图7中红色圆圈标记区域，仿真结果的反射率因子数值可达到40dbz左右，kfws地基雷达实测数据约为41dbz，说明降水粒子雷达回波仿真结果符合实际情况。
[0091]
图8(a)-(l)是六种降水粒子各自的水平反射率因子和差分反射率因子图，图(a)和图(b)为云滴的回波仿真结果，云滴的仿真能为降水的形成提供水汽条件，其相态变化可以转变为雨滴或冰晶等；从图(c)和图(d)可以看出雨区集中在风暴旋涡中心处，且差分反射率因子数据表明雨滴粒子形状各异，存在球形与扁椭球形；图(e)和图(g)显示了冰雹和霰的仿真结果，结果指出冰雹和霰通常会并存在同一区域；图(i)图(j)为雪晶的仿真结果，显示指出雪晶多存在于风暴外围区域，形状的不规则性导致差分反射率因子数值变化范围较大；图(k)和图(l)为冰晶的仿真结果，此次气象场景仿真，冰晶存在方位与冰雹和霰的区域较为一致，是冰雹形成的粒子基础。由于数值模式仅是通过求解动力方程，不能完全与降水场景情况一模一样，当前仿真结果只可以在宏观上符合实际降水粒子的分布，在局部仿真结果与实测数据仍有细微偏差存在。
[0092]
对于本领域的技术人员来说，可以根据以上的技术方案和构思，给出各种相应的改变和变形，而所有的这些改变和变形，都应该包括在本发明权利要求的保护范围之内。