CN115657013B - 基于激光雷达和云雷达的冰云内冰晶粒子数浓度估算方法 - Google Patents

Abstract

本发明涉及基于激光雷达和云雷达的冰云内冰晶粒子数浓度估算方法，属于气象雷达探测技术领域，考虑了激光和毫米波的大气辐射传输物理过程和两种遥感设备高时空分辨率观测的特点，在给定冰晶粒子数浓度、有效粒子半径和多种典型形状的冰晶粒子权重函数等初值的基础上，利用云解析雷达仿真模型计算激光和毫米波雷达的模拟观测值，通过对比模拟器输出的模拟观测值和两部雷达实际的观测值的差异，利用最优估计方法对冰晶粒子数浓度、有效粒子半径的初始设置值进行调整，基于价值函数最小化原则来判断反演是否收敛，进而实现混合态冰晶粒子数浓度及不同形状冰晶粒子数浓度的估算，最终获得的冰晶粒子数浓度的估算精度较高。

Description

基于激光雷达和云雷达的冰云内冰晶粒子数浓度估算方法

技术领域

本发明涉及气象雷达探测技术领域，尤其涉及一种基于激光雷达和云雷达的冰云内冰晶粒子数浓度估算方法。

背景技术

由多种不同形状的冰晶粒子组成的冰云对地球-大气***的辐射收支平衡和水汽循环有着重要的调节作用；因此，在大气气候模式的改进与优化过程中，冰云中冰晶粒子的微物理参量，如冰晶粒子的形状、尺寸和数浓度等的准确测量或估算显得越来越重要，这些微物理参量会直接影响冰云散射辐射特性的准确认识及数学建模，进而会影响气候模式输出结果的准确性；在冰晶微物理参数中，冰晶粒子数浓度（Ni）的准确测量显得尤为重要。Ni将气溶胶的变化特征与促进云形成的垂直对流的动力效应联系起来，即Ni的变化直接影响冰云的电磁波辐射特性和降水、云的发展及生命周期。因此，Ni的准确测量或估算是十分必要且有意义的。

然而，目前对Ni的最直接探测手段为利用飞机搭载测量仪器对冰云内部进行直接测量，但这种方法主要存在以下问题：1、因飞机测量实验本身观测时间的限制，其数据采样的时空分辨率较低，数据量小；2、冰云内冰晶粒子会受到飞机飞行扰动的影响而产生破碎现象，这就破坏了冰云内部冰晶原有的状态，导致Ni测量出现较大的偏差。此外，在观测非降水云时，Ka波段的毫米雷达是一种常用的仪器，同时基于这种仪器也开发了较多的云微物理反演算法。然而，受限于毫米波云雷达的灵敏度，仅仅依靠它是无法实现对冰云内尺度很小冰晶粒子的数浓度进行准确的探测。同时，已有一些研究是利用卫星的被动遥感资料对Ni进行反演，这种方法存在的最大缺点是卫星的空间分辨率很低，远不能实现对冰云高空间分辨率精细化探测的目的。而且现有对冰云内冰晶粒子数浓度估算的研究均没有考虑给出不同形状冰晶粒子混合状态下每种形状冰晶粒子数浓度的反演结果，而最后的冰晶粒子数浓度估算结果只体现了不同形状冰晶粒子群的平均分布状态，这种估算策略无法实现对每种典型形状冰晶粒子数浓度的精细化观测。

需要说明的是，在上述背景技术部分公开的信息只用于加强对本公开的背景的理解，因此可以包括不构成对本领域普通技术人员已知的现有技术的信息。

发明内容

本发明的目的在于克服现有技术的缺点，提供了一种基于激光雷达和云雷达的冰云内冰晶粒子数浓度估算方法，解决了现有冰晶粒子数浓度估算方法存在的不足。

本发明的目的通过以下技术方案来实现：基于激光雷达和云雷达的冰云内冰晶粒子数浓度估算方法，所述估算方法包括：

S1、获取激光雷达和毫米波云雷达数据进行处理，形成激光雷达和云雷达观测的多参数矢量，并获取卫星的激光雷达和毫米波云雷达联合反演的DARDAR云参数和NCEP的大气温度和湿度的大气廓线场；

S2、提取DARDAR和NCEP数据的经度、纬度和观测时间信息，结合激光雷达和云雷达的经度、纬度和观测时间信息，利用线性插值方法对DARDAR和NCEP数据进行时间和空间插值，获取激光雷达和云雷达观测位置的DARDAR和NCEP插值数据，进而生成与雷达观测时空匹配的大气温度、大气湿度、有效粒子半径和数浓度的数据，并将这些数据作为反演过程的初值；

S3、对组成冰晶粒子的各冰晶形状设置一权重初值，并设置冰晶粒子的有效半径、数浓度、形状和雷达频率，利用物理光学近似法计算激光雷达波段冰晶粒子散射特性，利用T矩阵模型获得毫米波段的冰晶粒子后向散射特性；

S4、将获取的与雷达时空匹配的NCEP大气温度和湿度廓线，以及计算的冰晶散射特性、激光雷达和毫米波云雷达的频率、功率、脉冲宽度、信噪比参数写入模型控制文件，将控制文件输入到云解析雷达仿真模型中，模拟计算得到激光雷达和毫米波云雷达的观测值；

S5、对比模拟和真实激光雷达与毫米波云雷达观测的差异，根据设定有效粒子半径和数浓度初值，以及与云解析雷达仿真模拟建立冰晶状态与雷达观测值的关系，并通过贝叶斯理论和Levenberg-Marquardt迭代法计算得到粒子数浓度值。

所述步骤S1具体包括以下内容：

S11、收集获取激光雷达和毫米波云雷达的数据，利用阈值方法对激光雷达和毫米波云雷达的有雨数据进行剔除，之后将激光雷达和毫米波云雷达进行时间和空间匹配，形成激光雷达和毫米波云雷达观测的多参数矢量；

S12、根据步骤S1获取的激光雷达和毫米波云雷达观测匹配数据集的时间和位置信息，获取卫星的激光雷达和毫米波云雷达联合反演的DARDAR云参数和NCEP的大气温度和湿度的大气廓线场。

所述物理光学近似法包括：根据粒子不同形状和激光束的轨迹，通过对麦克斯韦方程组的求解计算出粒子对光的后向散射特征，其计算公式为：，式中S₁是目标入射照明区，和s分别是入射方向的单位矢量，r表示距离，k表示波数，为坐标原点处的入射电场，代表点处的外向单位法向矢量，e为自动对数，j表示虚部，且j² =﹣1。

所述利用T矩阵模型获得毫米波段的冰晶粒子后向散射特性包括：将电磁波的入射场和散射场描述为矢量球形波函数，利用T矩阵关联函数展开项将冰晶粒子有效粒子半径、数浓度、形状、粒子谱分布和电磁波频率输入到T矩阵模型中计算出特定频率的冰晶粒子的后向散射特征参数。

所述冰晶状态与雷达观测值的关系为：y=F(x,b)+ε，其中x是待反演参数的状态向量，b是雷达仿真模型所需的其他参数，ε是雷达的观测误差，y是雷达观测值向量，F是雷达观测值向量y和雷达仿真模型所需的其他参数b用于模拟雷达观测的前向仿真模型。

所述步骤S5中通过贝叶斯理论和Levenberg-Marquardt迭代法计算得到粒子数浓度值具体包括：通过贝叶斯理论求解状态向量x转化为让价值函数最小化问题：，其中，是状态向量初值，是状态向量先验协方差矩阵，S_ε是观测误差协方差矩阵，y是雷达观测值向量，是价值函数；

通过Levenberg-Marquardt迭代法以逐步逼近冰晶粒子数浓度真值，并通过公式计算得到粒子数浓度值，其中，是第次迭代的状态向量结果，是第i次迭代的状态向量结果，是状态向量迭代第i次雷达仿真模型模拟的雷达观测向量，是Levenberg-Marquardt参数，是第i次迭代时计算的雅可比矩阵，经过多次迭代，实现价值函数最小化，最终实现冰晶粒子数浓度的估算。

本发明具有以下优点：一种基于激光雷达和云雷达的冰云内冰晶粒子数浓度估算方法，考虑了激光和毫米波的大气辐射传输物理过程和两种遥感设备高时空分辨率观测的特点，在给定冰晶粒子数浓度、有效粒子半径和多种典型形状的冰晶粒子权重函数等初值的基础上，利用云解析雷达仿真模型计算激光和毫米波云雷达的模拟观测值，通过对比模拟器输出的模拟观测值和两部雷达实际的观测值的差异，利用最优估计方法对冰晶粒子数浓度、有效粒子半径的初始设置值进行调整，基于价值函数最小化原则来判断反演是否收敛，进而实现混合态冰晶粒子数浓度及不同形状冰晶粒子数浓度的估算。估算过程中使用了激光雷达和双偏振毫米波云雷达的联合观测，因此获得的冰晶粒子数浓度的估算精度较高。此外，算法同时给出多种常见形状的冰晶粒子权重分布，更加接近实际冰云中冰晶粒子是由多种形状冰晶混合的情况。

附图说明

图1 为本发明的流程示意图；

图2 为本发明通过物理光学近似法后向散射计算方法的流程示意图；

图3 为本发明通过T矩阵方法后向散射计算方法的流程示意图。

实施方式

为使本申请实施例的目的、技术方案和优点更加清楚，下面将结合本申请实施例中附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本申请实施例的组件可以以各种不同的配置来布置和设计。因此，以下结合附图中提供的本申请的实施例的详细描述并非旨在限制要求保护的本申请的保护范围，而是仅仅表示本申请的选定实施例。基于本申请的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本申请保护的范围。下面结合附图对本发明做进一步的描述。

本发明具体涉及一种基于激光雷达和双偏振毫米波云雷达的冰云内冰晶粒子数浓度估算方法，该方法包括：收集激光和毫米波双偏振云雷达的数据，对激光和双偏振毫米波云雷达的有雨数据进行剔除，将激光和毫米波云雷达进行时间和空间匹配，形成激光和毫米波云雷达观测的多参数矢量。根据激光和毫米波云雷达观测匹配数据集的时间和位置信息，获取卫星的激光和毫米波云雷达联合反演的DARDAR(liDAR/raDAR)云参数和美国国家环境预报中心（NCEP）的大气温度、湿度大气廓线场；对DARDAR和NCEP数据进行时间和空间的二维线性插进行时空匹配，生成与雷达观测时空匹配的大气温度、大气湿度、有效粒子半径和数浓度的数据，将这些数据作为反演过程的初值；假定冰晶粒子由多种冰晶形状组成，设置各冰晶粒子形状的权重初值；利用物理光学近似计算激光雷达波段冰晶粒子散射特性，利用T矩阵方法计算毫米波波段的冰晶粒子散射特性；利用云解析雷达模拟模型（CR-SIM）计算激光雷达和双偏振毫米波云雷达的观测；对比模拟和真实激光雷达和毫米波双偏振云雷达观测的差异，利用最优估计模型调整有效粒子半径和数浓度的初值，同时计算价值函数，通过牛顿迭代方法对有效粒子半径和数浓度进行多次迭代，实现价值函数最小化，最终实现冰晶粒子数浓度的估算。

如图1-图3所示，具体包括以下内容：

步骤1：收集激光和双偏振毫米波云雷达的数据，利用阈值方法对激光和双偏振毫米波云雷达的有雨数据进行剔除，之后将激光和毫米波云雷达进行时间和空间匹配，形成激光和毫米波云雷达观测的多参数矢量。

步骤2：根据步骤1获取的激光和毫米波云雷达观测匹配数据集的时间和位置信息，获取卫星的激光和毫米波云雷达联合反演的DARDAR(liDAR/raDAR)云参数和美国国家环境预报中心（NCEP）的大气温度、湿度大气廓线场。

步骤3：对DARDAR和NCEP数据进行时间和空间的二维线性插进行时空匹配，生成与步骤1获取的雷达观测时空匹配的大气温度、大气湿度、有效粒子半径和数浓度的数据，将这些数据作为反演过程的初值。

步骤4：假定冰晶粒子由实心柱状、子弹花瓣状、聚集物状、空心柱状和平板状形状等5种形状的冰晶组成，设置各形状冰晶粒子的权重初值，各类型权重初值均设置为0.2；将步骤3的冰晶粒子有效粒子半径、数浓度、形状和雷达频率、观测方式作为输入数据，利用物理光学近似计算激光雷达532nm波段冰晶粒子散射特性，利用T矩阵方法计算Ka波段毫米波云雷达的冰晶粒子后向散射特性。

所述的物理光学近似法是根据粒子不同形状和激光束的轨迹，通过麦克斯韦方程组的求解计算出粒子对光的后向散射特征，具体的计算公式如下：，式中S₁是目标入射照明区，和s分别是入射方向的单位矢量，r表示距离，k表示波数，为坐标原点处的入射电场，代表点处的外向单位法向矢量，e为自动对数，j表示虚部，且j² =﹣1。

T矩阵方法包括：将电磁波的入射场和散射场描述为矢量球形波函数，利用T矩阵关联函数展开项将冰晶粒子有效粒子半径、数浓度、形状、粒子谱分布和电磁波频率输入到T矩阵模型中计算出特定频率的冰晶粒子的后向散射特征参数。

其中，冰晶粒子后向散射特性包括：激光532nm和毫米波Ka波段的冰晶粒子散射相函数、散射截面、消光系数、吸收系数、散射系数和消光效率。

步骤5：将步骤2获取NCEP大气温度、湿度廓线、步骤4计算的的冰晶散射特性、激光和毫米波云雷达的频率、功率、脉冲宽度、信噪比参数输入到云解析雷达仿真模型(CR-SIM)，模拟计算得到激光和双偏振毫米波云雷达的观测值。其中激光雷达波长设置为532nm，毫米波云雷达设置为Ka波段，观测方式为垂直指向。

步骤6：对比模拟和真实532nm激光雷达和Ka波段双偏振毫米波云雷达观测的差异，假设给定有效粒子半径和数浓度的初值，根据云解析雷达仿真模型可以建立冰晶状态与雷达观测值的关系，如式y=F(x,b)+ε，其中x是待反演参数的状态向量，b是雷达仿真模型所需的其他参数，ε是雷达的观测误差，y是雷达观测值向量，F是雷达观测值向量y和雷达仿真模型所需的其他参数b用于模拟雷达观测的前向仿真模型。基于贝叶斯理论，求解状态向量x转化为让价值函数最小化的问题：，式中，是状态向量初值，是状态向量先验协方差矩阵，S_ε是观测误差协方差矩阵，y是雷达观测值向量，是价值函数。

进一步使用Levenberg-Marquardt迭代法以逐步逼近冰晶粒子数浓度真值。粒子数浓度值由以下公式迭代计算得到：，式中，是第次迭代的状态向量结果，是第i次迭代的状态向量结果，是状态向量迭代第i次雷达仿真模型模拟的雷达观测向量，是Levenberg-Marquardt参数，是第i次迭代时计算的雅可比矩阵，经过多次迭代，实现价值函数最小化，最终实现冰晶粒子数浓度的估算。

以上所述仅是本发明的优选实施方式，应当理解本发明并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本文所述构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离本发明的精神和范围，则都应在本发明所附权利要求的保护范围内。

Claims (6)

1.基于激光雷达和云雷达的冰云内冰晶粒子数浓度估算方法，其特征在于：所述估算方法包括：

S1、获取激光雷达和毫米波云雷达数据进行处理，形成激光雷达和毫米波云雷达观测的多参数矢量，并获取卫星的激光雷达和毫米波云雷达联合反演的DARDAR云参数和NCEP的大气温度和湿度的大气廓线场；

S2、提取DARDAR和NCEP数据的经度、纬度和观测时间信息，结合激光雷达和云雷达的经度、纬度和观测时间信息，利用线性插值方法对DARDAR和NCEP数据进行时间和空间插值，获取激光雷达和云雷达观测位置的DARDAR和NCEP插值数据，进而生成与雷达观测时空匹配的大气温度、大气湿度、有效粒子半径和数浓度的数据，并将这些数据作为反演过程的初值；

S3、对组成冰晶粒子的各冰晶形状设置一权重初值，并设置冰晶粒子的有效半径、数浓度、形状和雷达频率，利用物理光学近似法计算激光雷达波段冰晶粒子散射特性，利用T矩阵模型获得毫米波段的冰晶粒子后向散射特性；

S4、将获取的与雷达时空匹配的NCEP大气温度和湿度廓线，以及计算的冰晶散射特性、激光雷达和毫米波云雷达的频率、功率、脉冲宽度、信噪比参数写入模型控制文件，将控制文件输入到云解析雷达仿真模型中，模拟计算得到激光雷达和毫米波云雷达的观测值；

S5、对比模拟和真实激光雷达与毫米波云雷达观测的差异，根据设定有效粒子半径和数浓度初值，以及与云解析雷达仿真模拟建立冰晶状态与雷达观测值的关系，并通过贝叶斯理论和Levenberg-Marquardt迭代法计算得到粒子数浓度值。

2.根据权利要求1所述的基于激光雷达和云雷达的冰云内冰晶粒子数浓度估算方法，其特征在于：所述S1的步骤具体包括以下内容：

S11、收集获取激光雷达和毫米波云雷达的数据，利用阈值方法对激光雷达和毫米波云雷达的有雨数据进行剔除，之后将激光雷达和毫米波云雷达进行时间和空间匹配，形成激光雷达和毫米波云雷达观测的多参数矢量；

S12、根据步骤S1获取的激光雷达和毫米波云雷达观测匹配数据集的时间和位置信息，获取卫星的激光雷达和毫米波云雷达联合反演的DARDAR云参数和NCEP的大气温度和湿度的大气廓线场。

3.根据权利要求1所述的基于激光雷达和云雷达的冰云内冰晶粒子数浓度估算方法，其特征在于：所述物理光学近似法包括：根据粒子不同形状和激光束的轨迹，通过对麦克斯韦方程组的求解计算出粒子对光的后向散射特征，其计算公式为：，式中S₁是目标入射照明区，和s分别是入射方向的单位矢量，r表示距离，k表示波数，为坐标原点处的入射电场，代表点处的外向单位法向矢量，e为自动对数，j表示虚部，且j² =﹣1。

4.根据权利要求1所述的基于激光雷达和云雷达的冰云内冰晶粒子数浓度估算方法，其特征在于：所述利用T矩阵模型获得毫米波段的冰晶粒子后向散射特性包括：将电磁波的入射场和散射场描述为矢量球形波函数，利用T矩阵关联函数展开项将冰晶粒子有效粒子半径、数浓度、形状、粒子谱分布和电磁波频率输入到T矩阵模型中计算出特定频率的冰晶粒子的后向散射特征参数。

5.根据权利要求1所述的基于激光雷达和云雷达的冰云内冰晶粒子数浓度估算方法，其特征在于：所述冰晶状态与雷达观测值的关系为：y=F(x,b)+ε，其中x是待反演参数的状态向量，b是雷达仿真模型所需的其他参数，ε是雷达的观测误差，y是雷达观测值向量，F是雷达观测值向量y和雷达仿真模型所需的其他参数b用于模拟雷达观测的前向仿真模型。

6.根据权利要求5所述的基于激光雷达和云雷达的冰云内冰晶粒子数浓度估算方法，其特征在于：所述S5的步骤中通过贝叶斯理论和Levenberg-Marquardt迭代法计算得到粒子数浓度值具体包括：

通过贝叶斯理论求解状态向量x转化为让价值函数最小化问题：，其中，是状态向量初值，是状态向量先验协方差矩阵，S_ε是观测误差协方差矩阵，y是雷达观测值向量，是价值函数；

通过Levenberg-Marquardt迭代法以逐步逼近冰晶粒子数浓度真值，并通过公式计算得到粒子数浓度值，其中，是第次迭代的状态向量结果，是第i次迭代的状态向量结果，是状态向量迭代第i次雷达仿真模型模拟的雷达观测向量，是Levenberg-Marquardt参数，是第i次迭代时计算的雅可比矩阵，经过多次迭代，实现价值函数最小化，最终实现冰晶粒子数浓度的估算。

Images