一种液态云微物理参数的反演方法与流程

本发明涉及一种利用毫米波云雷达获取的反射率因子反演液态云微物理参数的方法，特别涉及一种结合最优估计理论和经验公式的液态云微物理参数的反演方法。

背景技术：

云的全球分布及其微物理结构影响全球气候与环境变化，且对地球-大气系统的辐射收支平衡和水汽循环有着重要调节作用，通过观测全球云层特征，得到从薄云到浓厚云的垂直剖面特征，反演得到测量云水含量、云粒子尺寸以及其他参量，对气候研究、精细化天气过程分析、天气预测、人工影响天气等方面具有重要意义。

随着测云技术的快速发展，国内外学者不断深入研究云微物理参数的反演方法，取得了一些应用成果。传统的云物理参数反演主要借助统计分析得到的经验关系算法来完成，即以大量的实验数据为基础，利用函数拟合的方式建立雷达反射率与云物理参数的关系。

例如Atlas、Sauvageot、Kropfli和Kelly、Fox和Illingworth等人结合35GHz雷达和飞机穿云实测谱参数数据，得到雷达反射率与粒子有效半径、液水含量之间的函数关系，为指数形式；

Neil等人利用8mm地基云雷达探究了层云的云水含量和有效粒子半径大小，验证雷达反射率因子、粒子有效半径、云内液态水含量三者之间关系，结果与Atlas和Sauvageot的结果吻合。

Austin和Stephens基于最优估计理论，开发了一种基于最优估计的液态云微物理参数的业务化反演算法，利用先验数据信息，在状态向量(待反演参量)和先验数据差分以及已知测量向量(雷达反射率因子)和前向模型值差分的加权之和最小时，得到液态云微物理参数。

现有技术中常用的经验反演的优点在于计算简单，需要大量的探测数据进行统计分析，且根据观测地点、时间和云类的区别关键反演系数有所不同，但算法适用性和扩展性较差，例如，当云内存在毛毛雨降水时，由于降水粒 子对云内液水含量的贡献很小，利用Atlas经验公式反演得到的结果偏大，需要统计出适用于毛毛雨和小雨的关系式。最优估计反演精度受到先验数据的影响，目前采用观测统计值为先验数据，存在一定的局限性和不足。

专利方面，安徽师范大学提出“一种利用A-Train系列卫星数据协同反演云相态和云参量的新方法”(公开号：CN102707336A)与中国气象局气象探测中心提出了“一种云雷达与卫星探测数据融合方法及系统”(公开号：CN105445816A)，系统采用神经网络训练，建立亮温与云顶高度、云底高度和反射率因子的关系；北京无线电测量研究室提出了“一种毫米波云雷达的数据融合方法及系统”(公开号：CN104345312A)，系统融合雷达多模式工作方式得到的探测数据，提高数据质量和探测效率。此外，在探测手段方面，南京信息工程大学申请了“一种云粒子探测方法及探测器”(公开号：CN105115862A)的专利，利用云粒子的散射信号探测得到云粒子相态和云粒子大小；中国空间技术研究院申请了“一种用于测云的地面太赫兹雷达系统”(公开号：CN104569980A)以及“一种基于太赫兹主动测云雷达的测云实验装置及方法”(公开号：CN104597454A)。在所检索到的相关专利中，国内外的相关研究主要集中在云探测设备或仪器，部门单位提出了多源数据融合处理方法，云反演方法的研究较少。

技术实现要素：

本发明的目的根据毫米波云雷达提供的雷达反射率因子，基于最优估计理论，采用经验公式计算值为先验值，假设粒子谱服从对数正态分布，建立雷达反射率因子与液态云物理参数的函数关系，在待反演参量和先验数据以及雷达反射率因子和函数计算值的差值加权取得最小值的情况下求得反演最优解。

为了实现以上目的，本发明是通过以下技术方案实现的：

一种液态云微物理参数的反演方法，反演算法在代价函数取得最小值的情况下求得最优解，所述代价函数为状态向量和先验数据差分以及已知测量向量和前向模型值差分的加权之和，

式中，x_a为根据雷达反射率因子计算得到的先验数据，S_a为先验数据的协方差矩阵，F(x)为雷达反射率因子的前向物理模式，S_y是雷达反射率因子 测量误差的协方差矩阵，y为测量向量即雷达反射率因子，x为未知的状态向量即待反演参量，所述待反演参量包括云粒子数密度、几何平均半径和分布宽度参数。

将先验数据x_a作为迭代初值，通过连续最小化代价函数D，最终求得待反演参量x的迭代解。

优选地，假设大气中的云粒子的尺度分布符合对数正态分布，依此建立雷达反射率因子与液态云物理参数即云粒子数密度，云粒子有效半径，云粒子分布宽度参数的函数关系；

式中：N_T是云粒子数密度；r是云粒子有效半径；r_g是几何平均半径；σ_log是分布宽度参数，为无量纲变量；σ_g为几何标准差；ln代表自然对数变换；横线表示求算术平均值。

在未降雨或有小雨，或有毛毛雨的情况下，认为云粒子尺度足够小，满足瑞利散射条件，根据云内液态水含量C_LW、云粒子有效半径r_e和雷达反射率因子Z的定义，推导得到下列公式。

式中：ρ_ω代表水密度。

优选地，所述先验数据x_a与先验数据的协方差矩阵S_a根据由历史统计数据得到粒子数密度N_T和分布宽度σ_log来计算得到，

式中，距离库z_i；z₁和z_n分别代表雷达廓线中云底和云顶的距离库，下标a表示状态向量x的先验值，r_ga、N_Ta和σ_loga是分别代表云粒子有效半径、粒子数浓度和分布宽度的先验值。

优选地，通过云粒子尺寸即云粒子有效半径r_e、粒子数浓度N_T和分布宽度σ_log，推算出所述雷达反射率因子的前向物理模式F(x)；

通过前向物理模式F(x)建立测量向量y与状态向量x两者之间的关系，前向物理模式可表示为

y＝F(x)+ε_y

式中ε_y表示测量误差，Z_FM(z_i)为各距离库z_i的雷达反射率因子；z₁和z_n分别代表雷达廓线中云底和云顶的距离库。

优选地，通过雷达反射率因子的定义公式推导计算各距离库z_i的雷达反射率因子Z_FM(z_i)；

K＝(m²-1)/(m²+2)

式中，m表示复折射指数。优选地，所述测量向量y与未知的状态向量x为

式中：Z_dB(z_i)、r_g(z_i)、N_T(z_i)和σ_log(z_i)分别代表雷达反射率因子、几何平均半径、云粒子数密度和分布宽度参数在距离库z_i处的值；n代表廓线内的云距离库数。

优选地，对于每条观测廓线，状态向量x为每个距离库内的待反演参量，已知测量向量y仅有对应距离库内的雷达反射率因子，需要借助先验数据x_a对反演进行约束，保证迭代收敛以及反演结果的可靠性。

优选地，所述将先验数据x_a作为迭代初值，通过连续最小化代价函数D，求得待反演参量x的迭代解，且在进行迭代计算时满足迭代的收敛条件；

式中：上标i和i+1表示迭代次数，L代表前向物理模式对状态向量x的灵敏度。

优选地，所述迭代的收敛条件为：

式中：S_x是迭代状态向量的误差协方差矩阵，表示三个待反演物理量的方差以及各参量之间的协方差，S_y是雷达反射率因子测量误差的协方差矩阵。

与现有技术相比，本发明有益效果是：

本发明可弥补传统的经验公式适用性差的缺点，而且除了常见的云粒子半径和云水含量，还可得到粒子数浓度和分布宽度，反演的到的结果更为全 面；由于采用了实时观测的雷达反射率因子来计算先验值，提高了反演结果精确性。

附图说明

图1是本发明一种液态云微物理参数的反演方法的流程图；

图2是本发明一种液态云微物理参数的反演方法的一个实施例中采用的云粒子半径先验数据示意图；

图3是本发明一种液态云微物理参数的反演方法的一个实施例中反演得到的液态云微物理参数示意图。

具体实施方式

以下结合附图，通过详细说明一个较佳的具体实施例，对本发明做进一步阐述。

毫米波云雷达利用云粒子对电磁波的散射特性，通过对云的雷达回波分析了解云的宏观和微观特性，即可探测直径从几微米到弱降水粒子，连续检测云的垂直剖面。回波强度反映了云中粒子的大小和浓度，回波强度在时间和空间上的变化反映了云微物理过程演变特征。对雷达所测的反射率因子进行预处理，即可得到有云区域内的连续、有效的廓线和剖面回波数据，及其误差信息。

如图1所示，一种液态云微物理参数的反演方法，

反演算法在代价函数取得最小值的情况下求得最优解，所述代价函数为状态向量和先验数据差分以及已知测量向量和前向模型值差分的加权之和，

式中，x_a为根据雷达反射率因子计算得到的先验数据，S_a为先验数据的协方差矩阵，F(x)为雷达反射率因子的前向物理模式，S_y是雷达测量误差的协方差矩阵，y为测量向量即雷达反射率因子，x为未知的状态向量即待反演参量，所述待反演参量包括云粒子数密度、几何平均半径和分布宽度参数。

将先验数据x_a作为迭代初值，通过连续最小化代价函数D，最终求得待反演参量x的迭代解。

假设大气中的云粒子的尺度分布符合对数正态分布，依此建立雷达反射率因子与液态云物理参数即云粒子数密度，云粒子有效半径，云粒子分布宽 度参数的函数关系；

式中：N_T是云粒子数密度；r是云粒子有效半径；r_g是几何平均半径；σ_log是分布宽度参数，为无量纲变量；σ_g为几何标准差；ln代表自然对数变换；横线表示求算术平均值。

对于先验数据获取及计算，根据经验公式，由雷达有效反射率因子计算得到云粒子半径和液水含量。前人利用云雷达以及其它云粒子谱探测设备，测出云粒子的滴谱分布、粒子数浓度、粒子有效半径和云内液态水含量，进行大量数据的统计拟合得出云雷达反演粒子半径和液态水的经验关系，得出经典的常用经验公式，

Z＝138×10^-12D⁶ (4)

其中C_LW代表云内液态水含量，D为粒子直径，Z是毫米波云雷达测得的雷达反射率因子。

在未降雨或有小雨，或有毛毛雨的情况下，认为云粒子尺度足够小，满足瑞利散射条件，根据云内液态水含量C_LW、云粒子有效半径r_e和雷达反射率因子Z的定义，推导得到下列公式。

式中：ρ_ω代表水密度。

所述先验数据x_a与先验数据的协方差矩阵S_a根据由历史统计数据得到粒子数密度N_T和分布宽度σ_log来计算得粒子数浓度和分布宽度由历史统计数据得到，并依此计算先验数据的协方差矩阵。

式中，距离库z_i；z₁和z_n分别代表雷达廓线中云底和云顶的距离库，下标a表示先验值，r_ga代表云粒子有效半径的先验数据，同理，N_Ta和σ_lo_ga是分别代表粒子数浓度和分布宽度的先验值。

所述的雷达反射率因子的前向物理模式F(x)是通过已知的云粒子尺寸即云粒子有效半径r_e、粒子数浓度N_T和分布宽度σ_log，推算出雷达反射率因子的前向物理模式F(x)；即考虑毫米云雷达接收给定距离的云粒子散射能量时，所经路径上的双向衰减，得到经衰减订正后的雷达反射率因子Z_FM；通过雷达反射率因子Z的定义公式推导计算各距离库z_i的雷达反射率因子Z_FM(z_i)；

K＝(m²-1)/(m²+2) (12)

式中，z₁和z_n分别代表雷达廓线中云底和云顶的距离库，下标FM表示前向物理模式计算值，m表示复折射指数。

通过前向物理模式F(x)建立测量向量y与待反演参量x两者之间的关系，前向物理模式可表示为

y＝F(x)+ε_y (13)

式中ε_y表示测量误差。

所述测量向量y与未知的状态向量x为

式中：Z_dB(z_i)、r_g(z_i)、N_T(z_i)和σ_log(z_i)分别代表雷达反射率因子、几何平均半径、云粒子数密度和分布宽度参数在距离库z_i处的值；n代表廓线内的云距离库数。

对于每条观测廓线，状态向量x为每个距离库内的待反演参量，已知测量向量y仅有对应距离库内的雷达反射率因子，需要借助先验数据x_a对反演进行约束，保证迭代收敛以及反演结果的可靠性。

所述将先验数据x_a作为迭代初值，通过连续最小化代价函数D，求得待反演参量x的迭代解，且在进行迭代计算时满足迭代的收敛条件；

式中：上标i和i+1表示迭代次数，L代表前向物理模式对状态向量x的灵敏度。

所述迭代的收敛条件为：

式中：S_x是迭代状态向量的误差协方差矩阵，表示三个待反演物理量的方差以及各参量之间的协方差，S_y是云雷达测量误差的协方差矩阵。

反演误差主要来源于雷达反射率因子和先验数据获取算法，根据误差传递理论，计算出反演的待反演参数的标准差和误差，也称为不确定度。

根据所述液态云微物理参数的反演方法，采用CloudSat实测数据实现反演过程，并分析反演结果。所选反演个例为层积云，反演采用的云粒子半径先验数据如图2所示，与现有反演过程中采用的固定先验数据相比(即有效半径先验数据为固定值，且不随高度和经纬度变化)，由实测雷达反射率因子计算得到的有效粒子半径更能反映云物理参数的真实特性。

反演得到的液态云微物理参数结果如图3所示。云微物理参数包括有效粒子半径、液态水含量、粒子数浓度和分布宽度参数，有效粒子半径和云中液态水含量的最大值出现在云层中部雷达回波较强的区域，云粒子数浓度呈现随高度增加逐渐增大的趋势，分布宽度参数则随高度增加逐渐减小，各参数的空间分布和变化趋势与CloudSat发布结果基本一致，且符合国内外学者统计的晴天积云的云滴特征参数范围，反演结果可信。

尽管本发明的内容已经通过上述优选实施例作了详细介绍，但应当认识到上述的描述不应被认为是对本发明的限制。在本领域技术人员阅读了上述内容后，对于本发明的多种修改和替代都将是显而易见的。因此，本发明的保护范围应由所附的权利要求来限定。