一种基于星载三频毫米波雷达的云雨微物理参数反演方法与流程

本发明属于遥感技术领域，特别是一种基于星载三频毫米波雷达的云雨微物理参数反演方法。

背景技术：

云雨微物理参数的准确测量是气象准确预报的必要保证，是分析各种天气现象发生的关键决策辅助，对于人工影响天气作业也起到至关重要的作用。对云雨的探测主要分为直接测量和遥感探测两种。其中，直接测量方式包括升放无线电探空仪和机载实验测量云雨滴谱等；遥感探测方式是使用各个波段的雷达、激光雷达和辐射计等遥感仪器。其中，投放探空仪和机载实验测量实施物力和人力成本较大，所得结果较为精确，但是时空分辨率较低。用雷达进行探测能够连续不断地获得时空分辨率较高的回波数据，并可以基于相关特定算法来反演云雨物理参数。目前，地基探测主要利用天气雷达，以s波段和c波段雷达为主进行对云雨物理参数的探测。

相比较厘米波而言，毫米波雷达波长更接近云等小粒子尺度，且其较窄的波束宽度也更利于探测云和弱降水的结构特征。目前，已经研发出w波段(94ghz)测云雷达和ka(35ghz)/ku(13.6ghz)波段测雨雷达，并已投入使用，如目前在轨的cloudsat星载云廓线雷达和gpm星载测雨雷达。无论是测云雷达还是测雨雷达，现有的反演算法的基础都只是针对单种类型云或雨粒子，用简易化的粒子谱函数来拟合分布，并不能精细化地同时反映云雨粒子的实际分布。由于星载w，ka和ku三频毫米波雷达同时包括了w波段和ka/ku波段，所以在同时获取云雨微物理参数方面具有较大的优势。目前尚无针对星载w，ka和ku三频毫米波雷达反演云雨微物理参数的方法，因此，亟需寻找一种星载w，ka和ku三频毫米波雷达反演云雨微物理参数的方法。

技术实现要素：

本发明的目的在于提供一种时空分辨率高、反演结果精细、能够同时反演云雨粒子的实际分布的基于星载w、ka和ku三频毫米波雷达的云雨微物理参数反演方法。

实现本发明目的的技术解决方案为：一种基于星载三频毫米波雷达的云雨微物理参数反演方法，包括以下步骤：

步骤1、基于地基毫米波雷达和机载毫米波雷达的云雨探测数据，统计分析云雨滴粒子时空分布特征，拟合典型天气状况下的云雨滴谱函数并进行归纳，典型天气状况包括对流云降水、层状云降水；用广义gamma分布函数表征云雨滴粒子分布，并进行参数化；

步骤2、对星载三频毫米波雷达数据和地基毫米波雷达数据进行时空匹配；

步骤3、建立贝叶斯估计反演模型，对广义gamma分布谱函数中的任意一个谱参数预设先验概率分布，将星载w、ka和ku三频毫米波雷达反射率因子输入物理模型，得到其它三个谱参数；

步骤4、将地基毫米波雷达的时空匹配数据作为调整因子，输入贝叶斯估计反演模型进行后验，然后进行迭代计算，以修正先验概率分布，使代价函数最小，最后输出最优反演结果；

步骤5、进行云相态识别，计算在不同相态中每个距离库的云雨微物理参数，并进行线性化处理，得到整个廓线云雨微物理参数的反演结果；

步骤6、分析比对典型天气过程中的机载毫米波雷达云雨探测数据，将步骤4中反演得到的云雨微物理参数与所述云雨探测数据进行比对，计算偏差以验证结果的有效性。

进一步地，步骤1中所述的用广义gamma分布函数表征云雨滴粒子分布，并进行参数化，公式如下：

n(r)dr＝ar^μexp(-br^ν)dr

式中，n(r)为数密度分布函数，r为云雨滴粒子半径；n为粒子总数浓度，因式子中含有γ函数，故称为广义γ分布函数；该函数a是截距参数、b是尺度参数、μ是形状参数、ν是功率参数，四个参数均为正数，调节四个参数的值，拟合各种特征谱型的云雨滴谱。

进一步地，步骤3所述的建立贝叶斯估计反演模型，对广义gamma分布谱函数中的任意一个谱参数预设先验概率分布，将星载w、ka和ku三频毫米波雷达反射率因子输入物理模型，得到其它三个谱参数，具体如下：

步骤3.1、建立贝叶斯估计反演模型为：

式中，x为a、b、μ和ν四个参数中的一个参数，设定为已知量，服从预设的先验概率分布；y为测量向量，是在x服从先验概率分布条件下，输入星载w、ka和ku三频毫米波雷达的三个频段雷达反射率因子；ppr表示预设的先验概率密度函数，pf表示条件概率密度函数；

步骤3.2、联立方程组求得其余三个谱参数：

式中ze1、ze2和ze3分别为星载w、ka和ku三频毫米波雷达的三个频段经过衰减订正的雷达反射率因子值，λ1、λ2和λ3分别为星载w、ka和ku三频毫米波雷达的三个频段工作波长，σb1、σb2和σb3分别为星载w、ka和ku三频毫米波雷达的三个频段测得的后向散射截面；k为与复折射指数m有关的常数，d为粒子直径。

进一步地，步骤4中所述的进行迭代计算，以修正先验概率分布，使代价函数最小，具体如下：

迭代公式为：

式中，n为设定的廓线内的云距离库个数，x1和xn分别为在云低和云高距离库内的先验状态参数，a1和an、b1和bn、ν1和νn分别为在先验条件下云低和云高距离库内求得的其余三个谱参数；

反演算法在代价函数φ取得最小值的条件下求得最优解，φ为与x差分以及与y差分的权重之和，其中，为迭代计算后的状态参数，为迭代计算后的物理模型输出的谱参数，sx为先验状态参数向量的协方差矩阵，为物理模型输出的谱参数误差协方差矩阵。

进一步地，步骤5中所述的进行线性化处理，公式如下：

式中，同名像点α为线性化因子，t为温度，m为常数，根据云层相态信息进行取值，n1,n2...nn分别为各个有云距离库内经过反演模型得到的云雨微物理参数，包括粒子数浓度、水含量、水路径、有效粒子半径、消光系数和雨强。

本发明与现有技术相比，其显著优点在于：(1)使用雷达不间断获得回波数据，时空分辨率高；(2)经过迭代计算不断修正先验概率分布，使代价函数最小，反演结果精细，能够同时反演云雨粒子的实际分布；(3)利用获得的w，ka和ku三频的雷达回波数据，基于相关算法来反演云雨微物理参数，成本低。

附图说明

图1是本发明基于星载三频毫米波雷达的云雨微物理参数反演方法的流程示意图。

图2是本发明实施例中云雨微物理参数的反演结果曲线图。

具体实施方式

本发明基于星载三频毫米波雷达的云雨微物理参数反演方法，包括以下步骤：

步骤1、基于地基毫米波雷达和机载毫米波雷达的云雨探测数据，统计分析云雨滴粒子时空分布特征，拟合典型天气状况下的云雨滴谱函数并进行归纳，典型天气状况包括对流云降水、层状云降水；用广义gamma分布函数表征云雨滴粒子分布，并进行参数化；

步骤2、对星载三频毫米波雷达数据和地基毫米波雷达数据进行时空匹配；

步骤3、建立贝叶斯估计反演模型，对广义gamma分布谱函数中的任意一个谱参数预设先验概率分布，将星载w、ka和ku三频毫米波雷达反射率因子输入物理模型，得到其它三个谱参数；

步骤4、将地基毫米波雷达的时空匹配数据作为调整因子，输入贝叶斯估计反演模型进行后验，然后进行迭代计算，以修正先验概率分布，使代价函数最小，最后输出最优反演结果；

步骤5、进行云相态识别，计算在不同相态中每个距离库的云雨微物理参数，并进行线性化处理，得到整个廓线云雨微物理参数的反演结果；

步骤6、分析比对典型天气过程中的机载毫米波雷达云雨探测数据，将步骤4中反演得到的云雨微物理参数与所述云雨探测数据进行比对，计算偏差以验证结果的有效性。

进一步地，步骤1中所述的用广义gamma分布函数表征云雨滴粒子分布，并进行参数化，公式如下：

n(r)dr＝ar^μexp(-br^ν)dr

式中，n(r)为数密度分布函数，r为云雨滴粒子半径；n为粒子总数浓度，因式子中含有γ函数，故称为广义γ分布函数；该函数a是截距参数、b是尺度参数、μ是形状参数、ν是功率参数，四个参数均为正数，调节四个参数的值，拟合各种特征谱型的云雨滴谱。

进一步地，步骤3所述的建立贝叶斯估计反演模型，对广义gamma分布谱函数中的任意一个谱参数预设先验概率分布，将星载w、ka和ku三频毫米波雷达反射率因子输入物理模型，得到其它三个谱参数，具体如下：

步骤3.1、建立贝叶斯估计反演模型为：

式中，x为a、b、μ和ν四个参数中的一个参数，设定为已知量，服从预设的先验概率分布；y为测量向量，是在x服从先验概率分布条件下，输入星载w、ka和ku三频毫米波雷达的三个频段雷达反射率因子；ppr表示预设的先验概率密度函数，pf表示条件概率密度函数；

步骤3.2、联立方程组求得其余三个谱参数：

式中ze1、ze2和ze3分别为星载w、ka和ku三频毫米波雷达的三个频段经过衰减订正的雷达反射率因子值，λ1、λ2和λ3分别为星载w、ka和ku三频毫米波雷达的三个频段工作波长，σb1、σb2和σb3分别为星载w、ka和ku三频毫米波雷达的三个频段测得的后向散射截面；k为与复折射指数m有关的常数，d为粒子直径。

进一步地，步骤4中所述的进行迭代计算，以修正先验概率分布，使代价函数最小，具体如下：

迭代公式为：

式中，n为设定的廓线内的云距离库个数，x1和xn分别为在云低和云高距离库内的先验状态参数，a1和an、b1和bn、ν1和νn分别为在先验条件下云低和云高距离库内求得的其余三个谱参数；

反演算法在代价函数φ取得最小值的条件下求得最优解，φ为与x差分以及与y差分的权重之和，其中，为迭代计算后的状态参数，为迭代计算后的物理模型输出的谱参数，sx为先验状态参数向量的协方差矩阵，为物理模型输出的谱参数误差协方差矩阵。

进一步地，步骤5中所述的进行线性化处理，公式如下：

式中，同名像点α为线性化因子，t为温度，m为常数，根据云层相态信息进行取值，n1,n2...nn分别为各个有云距离库内经过反演模型得到的云雨微物理参数，包括粒子数浓度、水含量、水路径、有效粒子半径、消光系数和雨强。

下面结合附图及具体实施例对本发明作进一步详细说明。

实施例

结合图1，本发明基于星载三频毫米波雷达的云雨微物理参数反演方法，包括以下步骤：：

步骤1、基于大量地基和机载毫米波雷达云雨探测数据，对云雨滴粒子的时空分布特征进行统计分析，拟合各种典型天气状况下的云雨滴谱函数并进行归纳；

步骤2、对星地探测数据进行时空匹配，基于阈值法进行云检测；用广义gamma分布函数表征云雨滴粒子分布，并进行参数化，具体如下：

步骤2.1、对星载w、ka和ku三频毫米波雷达数据和地基毫米波雷达探测数据进行时间和空间上的匹配，要求每条廓线的时间差值不能超过5分钟，空间距离差不能超出2公里；利用星载微波辐射计获取得到的亮温数据，选择适当通道以设定阈值进行云检测；

步骤2.2、用广义gamma分布函数表征云雨滴粒子分布，并进行参数化，公式如下：

n(r)dr＝ar^μexp(-br^ν)dr(1)

式中n(r)为数密度分布函数，r为云雨滴粒子半径，n为粒子总数浓度，因式子中含有γ函数，故称为广义或修正γ分布函数；该函数a，b，μ和ν四个参数均为正数，且相互影响，调节四个参数的值，可以拟合各种特征谱型的云雨滴谱。

步骤3、建立贝叶斯估计反演模型，对其中一个谱参数预设先验概率分布，将星载w、ka和ku三频毫米波雷达反射率因子输入物理模型，得到其它三个谱参数，具体如下：

步骤3.1、针对星载三频毫米波雷达能够获得w、ka和ku三个频段雷达反射数据的特点，基于贝叶斯估计理论的思想，设定其中一个云雨滴谱参数服从预设的先验概率分布，建立贝叶斯估计反演模型，公式如下：

式中x为a、b、μ和ν四个参数中的一个参数，设定为已知量，服从预设的先验概率分布；y为测量向量，是在x服从先验概率分布条件下，输入星载w、ka和ku波段毫米波雷达的三个频段雷达反射率因子；

步骤3.2、联立方程组求得的其余三个谱参数，公式如下：

式中ze1、ze2和ze3分别为星载w、ka和ku三频毫米波雷达的三个频段经过衰减订正的雷达反射率因子值，λ1、λ2和λ3分别为星载w、ka和ku三频毫米波雷达的三个频段工作波长，σb1、σb2和σb3分别为星载w、ka和ku三频毫米波雷达的三个频段测得的后向散射截面，k为与复折射指数m有关的常数，d为粒子直径。

步骤4、将地基毫米波雷达探测云雨的时空匹配数据作为调整因子，输入贝叶斯估计反演模型进行后验，以修正先验概率分布，不断地调整反演模型；经过迭代计算之后，反演模型经过不断调整，使得代价函数最小，此时反演模型所得结果即为最优反演结果，公式如下：

式中，n为设定的廓线内的云距离库个数，x1和xn分别为在云低和云高距离库内的先验状态参数，a1和an、b1和bn、ν1和νn分别为在先验条件下云低和云高距离库内求得的其余被三个谱参数；

反演算法在代价函数φ取得最小值的条件下求得最优解，式中，φ为与x差分以及与y差分的权重之和，其中，为迭代计算后的状态参数，表示迭代计算后的物理模型输出的谱参数，sx为先验状态参数向量的协方差矩阵，为物理模型输出的谱参数误差协方差矩阵。

步骤5、根据在不同相态中每个云雨微物理参数的表达式，将最优反演结果谱函数代入贝叶斯估计反演模型，计算每个距离库的云雨微物理参数，以温度为判据进行云层的相态识别：当温度t<-20℃时，判定为冰云，反演结果全部采纳；当t>0℃时，判定为液态云，反演结果全部采纳；当-20℃<t<0℃时，判定为冰水混合相云，此时需要根据各个距离库内的相态变化引起的温度差异对每个库内的反演结果进行线性化处理，公式如下：

式中，同名像点α为线性化因子，t为温度，m为常数，根据云层相态信息进行取值，n1,n2...nn分别为各个有云距离库内经过反演模型得到的云雨微物理参数，包括粒子数浓度、水含量、水路径、有效粒子半径、消光系数和雨强。综合各相态的处理结果，得到整个廓线内云雨微物理参数的反演结果，如图2所示。

步骤6、分析比对典型天气过程中的机载毫米波雷达云雨探测数据，将步骤4中反演得到的云雨微物理参数与其进行比对，计算偏差，验证本研究方法的有效性，改进算法，具体如下：

由于机载实验更贴近深入云雨内部，所以获得的机载毫米波雷达回波数据更为准确，因此选择典型天气过程，将星地毫米波雷达与机载毫米波雷达探测云雨数据进行时空匹配，并以机载毫米波雷达数据作为真值，将基于贝叶斯估计理论建立的反演模型得到的云雨微物理参数与其进行比对析，计算偏差，验证本研究方法的有效性，改进算法。