基于变分分析的星载Ka/W双频云雨雷达粒子谱反演方法与流程

基于变分分析的星载ka/w双频云雨雷达粒子谱反演方法
技术领域
1.本发明涉及一种星载双频雷达粒子谱反演方法，尤其涉及一种基于变分分析的星载ka/w双频雷达的粒子谱反演方法。


背景技术：

2.星载ka/w双频云雨测量雷达可以测量ka和w两个频段的反射率因子，即 z
ka
和zw。由于ka和w两个不同频段的云降水粒子散射特性存在差异，相比与常规云雷达和降水雷达而言，双频段反射率因子联合观测可以获取更多的云雨粒子信息，为云、降水研究提供重要的数据基础。其中，自然界云降水中，双波长比dfr与粒子中值直径存在正相关。在实际应用中一般需要通过反演得到云雨粒子的谱分布，为天气研究提供云降水微物理特性的定量信息。我国已经开始初步研制星载ka和w双频云雨雷达，因此，反演算法的研制有利于未来的实际应用。
3.在ka与w两个频段上，雷达能量受云、降水的衰减作用明显，会给粒子谱的反演带来显著的不确定性。星载雷达具有独特的自上而下的观测模式，现有的降水雷达(如热带测雨卫星雷达)一般会装载ku波段雷达。ku波段的衰减相对于ka和w波段而言较小，一般能够观测到地表或海表面的反射率因子，在衰减订正时可以将地表/海表的反射率因子作为参考，其精度一般能满足研究需求。然而对于ka/w双频雷达而言，由于w波段的强衰减作用，一般在降水稍强时就很难观测到地表/海表的信号，给云雨粒子的反演带来难点。


技术实现要素：

4.本发明所要解决的技术问题是观测强降水时，ka和w频段雷达信号衰减较强导致观测不到地表信号，衰减订正困难，进而影响星载ka/w双频云雨雷达的粒子谱的反演。
5.为了解决上述技术问题，本发明的基于变分分析的星载ka/w双频云雨雷达粒子谱反演方法，包含如下步骤：
6.步骤1,接收星载ka/w双频雷达一个扫描径向上的ka反射率因子观测z
ka
和w波段反射率因子观测zw。
7.步骤2，根据星载ka/w双频雷达的扫描径向上各观测点的时间、经纬度和高度，利用数值模式数据进行三维线性插值,得到各观测点对应的温度t。
8.步骤3，根据温度t确定各观测点上的云降水粒子的相态、密度、形状、倾斜角方差、融化度和介电常数属性，并计算在ka和w波段云降水粒子的后向散射截面和消失截面。
9.步骤4，将云降水粒子谱分布假定为伽马分布，将云降水粒子谱分布模型参数中云降水粒子谱平均直径dm和浓度截距参数nw作为变分分析的状态向量，将z
ka
和zw作为观测量，构建变分分析的代价函数，其中观测算子中考虑云降水粒子谱和后向散射截面计算的反射率因子和单位距离衰减。
10.步骤5，利用lbfgsb算法最小化代价函数，并在其中引入状态变量为非负的条件，求得云降水粒子谱平均直径dm和浓度截距参数nw，以及云降水粒子谱。
11.上述技术方案中，当温度t≤0摄氏度时，云降水粒子假定为冰和空气的混合体，即冰晶和雪粒子，冰相粒子假定为扁球体，其形状由轴长比表示，轴长比为扁球体短轴比长轴；当温度t﹥0摄氏度时，云降水粒子假定为扁球体，其属性设定由高度h决定，高度单位为米。
12.上述技术方案中，当t≤-20摄氏度时，云降水粒子密度设定为0.4g/cm3，轴长比设定为0.5，云降水粒子倾斜角方差为20度；当温度t在-10﹤t≤0摄氏度时，云降水粒子密度设定为0.2g/cm3，轴长比设定为0.75，云降水粒子倾斜角方差为20度；当温度t在-20﹤t≤-10摄氏度之间时，云降水粒子密度设定为0.5 0.03
×
(t 20)g/cm3，轴长比设定为0.4 0.035
×
(t 20)，云降水粒子倾斜角方差为20度。
13.上述技术方案中，当温度t﹥0摄氏度时，将温度t＝0摄氏度的高度标记为 h0，单位为米，当h《h
0-1000时，云降水粒子为纯液态，密度为1g/cm3，其中轴长比ar假定与粒子的直径d的关系符合brandes等(2002)提出的关系 ar
brandes
(d)＝0.9951 0.02510d-0.03644d2 0.005030d
3-0.0002492d4，云降水粒子倾斜角为5度；而当h介于h0和h
0-1000时，粒子为半融化态，即冰、空气、液态水的混合体，当h0﹥h≥h
0-500时，云降水粒子的融化度为0.001
×
(h
0-h)，云降水粒子的密度为0.2 0.3/500
×
(h
0-h)g/cm3，直径为d的云降水粒子的轴长比设定为0.75 0.001
×
(h
0-h)
×
[ar
brandes
(d)-0.75]；当h
0-500﹥h≥h
0-1000 时，云降水粒子的融化度为0.001
×
(h
0-h)，云降水粒子的密度为 0.5 0.5/500
×
(h
0-h-500)g/cm3，直径为d的云降水粒子的轴长比设定为0.75 0.001
×
(h
0-h)
×
[ar
brandes
(d)-0.75]，云降水粒子倾斜角为20-0.015
ꢀ×
(h
0-h)。
[0014]
上述技术方案中，所述变分分析的代价函数为所述变分分析的代价函数为其中x是n个有效观测点上的状态变量dm和nw，x＝ [dm(1)，dm(2)，
…
，dm(n)，nw(1)，nw(2)，
…
，nw(n)]。和分别是有效观测点上的z
ka
和zw观测，和分别是z
ka
和zw对应的观测误差，和分别是从x计算z
ka
和zw的观测算子，f
dfr
是纯液态粒子时利用ka波段反射率因子z
ka
预测dfr的关系，of是n
×
n的对角阵，其中纯液态粒子观测点对应值设为 1，其他观测点上设为一个极大值，可设为10的8次方。
[0015]
上述技术方案中，所述数值模式数据选择中国气象局的全球模式数据或者最近时刻的欧洲中心第五代再分析数据。
[0016]
上述技术方案中，液态水的介电常数根据liebe等(1991)给出的方法计算，冰的介电常数根据hufford(1991)给出的方法计算。融化粒子的介电常数根据 maxwell-garnett(m-g)混合公式计算。粒子的后向散射截面和消失截面采用 t-矩阵算法计算。
[0017]
上述技术方案中，所述利用lbfgsb算法最小化代价函数时不采用线性化步骤，各迭代最小化过程中的梯度利用自动微分算法求得。
[0018]
本方案的基于变分分析的星载ka/w双频云雨雷达粒子谱反演方法解决了降水粒子对ka和w频段雷达信号较强所导致的观测不到地表信号时衰减订正困难的问题。
附图说明
[0019]
图1 2017年8月23日00:00华南hato台风个例wrf数值模式模拟对应的 ka(a)和w(b)波段星载雷达观测反射率因子观测，以及双波比dfr(c)。在仿真中采用wolfensberger和berne(2018)开发的仿真算子，雷达高度设定为 410km，经纬度分别为东经114.3度和北纬21.0度，雷达波束宽度0.14度，径向分辨率250m，交轨方向为方位角0度(南北)。降水对雷达电磁波的衰减已考虑，ka和w波段雷达的灵敏度假设为-12dbz。仿真中假定空气衰减已经经过提前订正。
[0020]
图2是图1对应的温度(a)和利用温度确定粒子密度的示意图(b)。
[0021]
图3是利用t矩阵算法计算的ka和w波段雪粒子的后向散射截面和粒子直径之间的关系，雪密度设置为0.2g/cm3，温度为0度。其中，介电常数是利用 maxwell-garnett混合方案计算。
[0022]
图4是与图1相对应的变分反演的结果，分别为z
ka
(a)和dfr(b)，以及反演得到的dm(c)。
[0023]
图5是基于变分分析的星载ka/w双频云雨雷达粒子谱反演方法。
具体实施方式
[0024]
下面结合附图对本发明的技术方案进行详细说明。通过仿真数据验证本方案的实现过程和性能。
[0025]
参见图5，具体步骤如下:
[0026]
步骤1，接收星载ka/w双频雷达的扫描数据，以如图1中利用数值模拟进行仿真的一次台风过程结果为例。其中仿真采用的方法见图注。
[0027]
步骤2，通过插值得到观测所对应的温度t，图2(a)为利用欧洲中心第五代再分析数据(era5)进行三维线性插值得到的各层的温度结果。
[0028]
步骤3，根据温度t确定各观测点上的云降水粒子的相态、密度、形状、倾斜角方差、融化度和介电常数属性，并计算在ka和w波段粒子的后向散射截面和消失截面。图2(b)给出了粒子的属性和温度与高度的关系，其中a-e分别为-20摄氏度、-10摄氏度、0摄氏度、0摄氏度以下500米、0摄氏度以下1000 米的节点，右侧分别为对应的密度、轴长比、粒子倾斜角方差；其余点分别是这些点关于高度的线性插值。图3给出了在节点c，利用t矩阵算法计算的后向散射截面与粒子直径的结果，黑色实线和虚线分别是ka和w波段的结果。
[0029]
步骤4，将粒子谱分布假定为伽马分布，构建变分分析的代价函数。该步骤中，dm和nw是伽马分布的两个参数，伽马分布的公式如下，其中μ设定为常数3。
[0030][0031][0032]
根据散射原理，z
e，ka
、z
e，w
、a
ka
和aw均写为关于各观测点的粒子谱n(d) [m-3
mm-1
]的函数，即
[0033][0034][0035][0036][0037]
其中λ是电磁波波长，d是降水物等效直径，k
w，ka
和k
w，w
是两个波段的介电因子(常数)，σ
b，ka
[mm2]和σ
b，w
[mm2]分别为ka和w波段的后向散射截面，σ
e，ka
[mm2] 和σ
e，w
[mm2]分别为ka和w波段的消失截面。n(d)是粒子谱，它表示单位取样体积范围内单位粒子直径档内的粒子数，是关于粒子的等效体积直径(d，单位为mm)的函数。由上述公式，我们可以得到观测前项算子
[0038][0039][0040]
同时，观测算子中引入了纯液态粒子时ka波段反射率因子z
ka
和dfr的关系作为约束，防止反演本身的数值不稳定，我们这里采用的关系为dfr(z
ka
)＝0.54z
ka
‑ꢀ
5.44。
[0041]
步骤5，利用lbfgs-b算法(limited-memorybroyden-fletcher-go1dfarb-shanno algorithm for bound optimization)计算步骤5中的代价函数最小时的x，同时在最小化时设置x中对应dm的元素大于0；在最小化时用到的代价函数的梯度我们利用autograd软件自动获取。
[0042]
步骤6，虽然图1中给出了整个剖面扫描图像，但反演是逐个雷达径向进行，直至扫描结束。图4给出了与图1相对应的变分反演的结果，图4(a)和图4 (b)为反演结果正演得到的z
ka
和dfr，图4(c)为反演得到的dm，结果相对比较合理。