1 引言
冰云作为气候系统的重要组成部分， 覆盖了地球上空20%~30%的范围， 其分布和演变对天气和气候变化研究具有重要影响， 冰云主要由各种球形、 非球形的冰晶粒子组成（Lynch， 2002； Heymsfield et al， 2017）， 其云高、 云厚等宏观参量， 以及云中冰晶粒子大小、 形状、 数浓度等微观特征的精度直接影响其散射辐射特性描述的准确性（Stephens， 2005； Mace et al， 2006； Loeb et al， 2018）， 对气候模式、 遥感探测、 大气环境等也有着不可忽视的影响（Sassen et al， 2001； 张琪等， 2018）， 其特性观测及反演研究有着非常重要的意义。
毫米波云雷达（35 GHz和95 GHz）是目前冰云的主要观测手段（Stephens et al， 2002； 仲凌志等， 2009； 高仲辉等， 2014； 郑晨雨等， 2020）， 为了能够更好的认识冰云特征， 许多学者开展了冰云物理参数的垂直分布及时空演变统计分析（李帅等， 2019； 陈纹锋等， 2019； 林彤等， 2018）， 探究高精度的云物理参数反演方法（Liu et al， 2000； Protat et al， 2007； 姚志刚等， 2018）。与毫米波相比， 太赫兹波的波长更接近冰云粒子尺寸， 在理论上具有更高的分辨力和灵敏度（Mead et al， 1990， 1989）， 特别是220 GHz大气窗口的电磁波， 适用于星载平台冰晶粒子以及含水量相对少高纬度地区或高空区域气象遥感探测， 能够实现冰云三维结构的高精度测量（Battaglia et al， 2014）。太赫兹云雷达的研发给冰云微物理参数反演方法的研发工作提出了新的研究课题， 以求实现太赫兹频段的冰云特征参量反演， 满足云物理研究、 气象保障的需求。因此， 针对太赫兹云雷达的技术发展需求， 本文开展了220 GHz冰云微物理参数反演方法的模拟研究。
受观测数据所限， 目前云反演研究集中在毫米波频段， 少有文献探讨太赫兹频段的云反演方法。基于最优估计理论的反演算法如今被广泛应用于毫米波雷达的冰云微物理参数反演研究， Benedetti et al（2003）提出将最优估计理论用于星载云雷达CloudSat/CPR的云反演算法， Austin et al（2009）在此基础上开发了一种应用于CloudSat/CPR的冰云微物理参数的业务化反演算法， 在该算法中， 假设粒子服从正态分布， 根据瑞利散射及其修正来建立前向模型。Delanoë et al（2008）、 Deng et al（2010）也基于该理论， 利用CloudSat和CALIPSO的探测数据联合反演得到冰云有效粒子半径和冰水含量。
本文利用最优估计方法， 进行太赫兹频段的球形及非球形冰晶粒子反演算法研究， 并验证其可靠性。在前向模型的建立过程中， 并没有像Austin et al（2009）的算法提出基于瑞利散射及修正来计算前向模型值， 而是采用Deng et al（2010）提出的算法， 利用220 GHz不同形状冰晶粒子的散射特性数据， 建立卷云冰晶粒子的散射特性库（Look Up Table， LUT）。在忽略冰晶粒子对太赫兹信号的衰减和多次散射等条件下， 基于最优估计理论， 反演得到云中冰水含量、 有效粒子半径， 并在Deng et al（2010）提出的算法基础上， 补充了粒子数密度的反演方法， 模拟验证算法的可靠性， 研究先验数据的偏离对于反演误差和计算效率的影响。
2 冰云微物理参数反演方法
2.1 太赫兹频段（220 GHz）冰晶粒子散射特性研究
云中冰晶粒子主要呈现盘状、 柱状、 针状、 树枝状、 星状以及子弹状， 且随着云高、 气压、 温度等因素的变化而变化， 粒子尺寸从几微米到上千微米不等。根据Baum et al（2005a）关于中纬卷云中冰晶粒子形状分布的统计结果可知， 当粒子最大尺度小于60 μm时冰晶粒子基本全为球形， 当粒子最大尺度为60~1000 μm时， 冰云由15%的子弹花、 50%的六棱柱和35%的六角板冰晶粒子组成， 本文选取球形（Sphere）、 六棱柱（Column）、 六角板（Plate）和子弹花（Rosette）形状的冰晶粒子作为研究对象。
采用离散偶极子近似法（Theoiserete-Dipole Approximation， DDA）计算球形、 六棱柱和六角板粒子在220 GHz的散射参数， 其中板状粒子轴比关系为l/d=0.2， 粒子最大尺度范围81~3246 μm， 柱状粒子的轴比关系为l/d=2， 最大尺度范围为83~3304 μm， 该方法已成功应用于各频段的云粒子散射特征研究（Draine et al， 1994； Hong， 2007； 吴举秀等， 2014）。子弹花形状的冰晶粒子散射参数来源于Liu（2008）的冰云粒子散射特性数据库， 该数据库同样由DDA计算得到， 粒子最大尺度范围50~10000 μm。忽略冰晶粒子对太赫兹信号的衰减和多次散射， 依次计算得到太赫兹频段下冰云雷达反射率因子、 冰水含量和有效粒子半径。
云物理中把云粒子浓度（单位体积中云滴的个数）随尺度的变化叫做云粒子谱， 在进行太赫兹频段冰云粒子特征反演研究时， 假设云粒子谱满足Gamma分布， 则有
式中： N T为单位体积内粒子数； D为有效粒子直径； D g为中值长度； α为伽马分布参数。雷达反射率因子Z e、 冰水含量（Ice Water Content， IWC）和有效粒子半径r e定义为
式中： ρi代表冰粒子密度； σ bk为冰云粒子的后向散射截面； K ice由粒子的复折射指数m计算得到。水和冰粒子的复折射指数Ray（1972）和Matzler（2006）总结的公式计算得到。
根据220 GHz频段四种形状粒子的后向散射特性模拟计算得到的r e、 N T与Z e、 IWC的关系， 建立云粒子散射特征数据表LUT， 用于冰云微物理参数r e、 IWC和N T的反演计算。
2.2 前向模型
对于太赫兹云雷达， 其实际观测值可计算如下
式中： σ ra为雷达消光系数， 但本文在计算过程中没有考虑云粒子分布随高度的变化， 利用单层数据进行反演算法验证， 因此该项可忽略。即不考虑雷达等效反射率因子随高度的衰减， 所得值为单个高度点上的观测值。该公式即为前向模型 F（ x）。
在反演过程中， 雷达反射率因子Z e为LUT中的模拟计算值， 对应的冰水含量IWC、 有效粒子半径r e和粒子数浓度N T将作为真值进行算法验证。
2.3 最优估计理论
测量向量 y为雷达反射率因子Z e， 未知的状态向量 x为待反演参数， 包含有效粒子半径r e和冰水含量IWC。通过前向物理模式 F（ x）建立两者之间的关系。前向物理模式可表示如下。
式中： ； ε y表示测量误差； z代表雷达廓线中单层的距离库。
反演算法在代价函数取得最小值的情况下求得最优解， 代价函数公式如下， 代表状态向量和先验数据差分以及已知测量向量和前向模型值差分的加权之和。
式中： x a为先验数据； S a为对应的协方差矩阵； S y为测量误差的协方差矩阵。
待反演参量x以先验数据x a（r ga IWC a）为初始值， 通过Gauss-Newton迭代， 并设置收敛条件， 求得距离库z上的待反演变量。
式中： 上标i和i+1表示迭代次数； Jacobian矩阵 L代表前向物理模式对状态向量的灵敏度， 由式（9） 计算得到矩阵 L， 不考虑衰减效应， 式中Z e对r e和IWC偏导可通过LUT查询得到。
2.4 粒子数密度 N T反演方法
上小节， 在假设N T=1的前提下进行有效粒子半径r e和冰水含量IWC反演验证， 本小节将进一步完成N T的反演。
由公式（9） 可知， 在迭代运算的过程中， 实现矩阵 L中Ze对r e的偏导计算时， 将Ze对IWC进行了归一化， 即 ， 由公式（2） 和（3）可知， N T与Z e、 IWC成正比， 即矩阵 L中Ze对r e的偏导项消除了N T的影响， 说明该算法并不局限于N T=1的情况， 可适用于N T为任意值的反演计算。在基于最优估计理论反演得到r ere和IWC re之后， 根据LUT得到r ere对应的IWC/N T值， IWC re已知， 可计算得到N Tre， 完成冰云微物理参数r ere、 IWC re和N Tre的反演计算。
2.5 先验数据和误差设置
对于每条云观测廓线， 测量变量为各距离库的雷达有效反射率因子Z e， 待反演变量为对应每个距离库内的冰水含量IWC和有效粒子半径r e， 反演变量多于测量变量， 需借助先验信息x a对反演过程进行约束， 保证迭代过程的收敛性。目前基于最优估计理论的反演算法所采用的先验数据多来源于云微物理参数测量数据集和经验关系， 云微物理参数测量数据集包含不同云类、 陆地云、 海洋云以及不同纬度的云微物理参数的统计分析结果。在反演过程中可通过调整x a和S a内各变量来实现对不同物理机制下各类云微物理参数的反演。
Austin et al（2001）在进行液态云粒子反演时， 采用Miles et al（2000）统计的云微物理参数数据集来设置先验数据x a（r ga、 N Ta、 σ loga）， 其中， r ga、 N Ta和σ loga分别为几何平均半径、 云粒子数密度和分布宽度的先验值， 在进行冰云粒子反演时， Austin et al（2009）利用σ log-T和r g-T的经验关系来得到r ga和σ loga， 随后利用Z-IWC的经验关系和前向模型设置N T先验值， 对应的协方差矩阵S a设为固定值。Deng et al（2010）则采用经验关系和CloudSat/CALIPSO的反演结果作为先验值x a（r ga， IWC a）， 先验数据协方差矩阵设置为2x a。
本文将LUT中雷达反射率因子Z e对应的计算得到的冰水含量IWC和有效粒子半径r e作为真值， 先验数据x a协方差矩阵设置为2x a， 为分析反演算法的可靠性， 将先验值不同程度的偏离真值， 分析先验值对反演精度和计算效率的影响。
3 算法性能分析
3.1 反演算法模拟验证
为模拟验证上述方法的可靠性， 对四种形状的冰晶粒子采用任意N T和r e值进行粒子散射特性运算， 得到对应的Z e和IWC的模拟值（表1）。以Z e作为测量向量进行最优估计反演运算， 在计算过程中， 先验数据x a（r ea， IWC a）设定为对应的有效粒子半径r e和冰水含量IWC的模拟值， 反演得到冰云微物理参数r ere、 IWC re和N Tre。在所选个例中， r e反演误差最大为2.51%， IWC反演误差最大值为3.26%， N T反演误差不超过12%， 证明了该方法可有效应用于冰云微物理参数r ere、 IWC re和N Tre的反演计算。
表1 220 GHz频段冰晶粒子随机模拟计算参数及其反演误差Table 1 Random simulated ice crystals properties at 220 GHz and their retrieval errors |
| 粒子形状 | 随机模拟计算结果 | 反演误差 | |||||||
|---|---|---|---|---|---|---|---|---|---|
| Z e/dBZ | r e/μm | IWC/(g·m-3) | N T/m-3 | r ere/% | IWC re/% | N Tre/% | |||
| 球形 | Case1 | -24.43 | 36 | 3.28×10-3 | 2000 | 0.24 | 1.26 | 2.20 | |
| Case2 | -32.19 | 20 | 3.13×10-3 | 20000 | 0.31 | 0.06 | 1.31 | ||
| 六棱柱 | Case1 | -1.69 | 40 | 1.25×10-2 | 5000 | 1.52 | 3.31 | 9.84 | |
| Case2 | 23.14 | 120 | 0.41 | 2000 | 0.03 | 1.68 | 1.56 | ||
| 六角板 | Case1 | -1.90 | 32 | 2.05×10-2 | 20000 | 1.47 | 3.26 | 9.57 | |
| Case2 | 27.64 | 160 | 1.28 | 2000 | 2.51 | 2.39 | 11.61 | ||
| 子弹花 | Case1 | 5.87 | 80 | 4.00×10-2 | 1000 | 1.44 | 1.16 | 7.23 | |
| Case2 | 2.06 | 52 | 3.58×10-2 | 5000 | 0.24 | 0.03 | 0.97 | ||
3.2 有效粒子半径 r e和冰水含量 IWC的反演误差分析
将LUT中的雷达反射率因子Ze作为测量向量， 对应的有效粒子半径r e和冰水含量IWC作为先验数据x a（r ea， IWC a）， 下标a表示先验值， 即先验数据等于真值， 假设N T为固定值， 按上述方法， 基于最优估计理论反演得到r ere和IWC re， 为区别与真值（r e和IWC）， 文中增加下标re表示反演值， 四种冰晶粒子的反演结果和真值如图3所示。反演结果与真值基本吻合， 说明在先验数据准确时， 该算法可准确反演球形和非球形冰晶粒子的微物理参数。
当先验数据x a（r ea， IWC a）中的r ea不同程度偏离真值时（IWC a等于真值）， r ere反演误差如图4（a）所示， 横坐标r ea/r e的比值从0.1增加到1， 表示先验数据逐步贴近真值。在r ea仅为真值r e的1/10时， 反演误差值最大， 其中球形粒子的反演误差的平均值为12.78%， 六棱柱的平均反演误差值为13.91%， 六角板的平均反演误差值为17.02%， 子弹花的平均反演误差值最小， 仅为1.39%。平均误差随r ea/r e的比值增加而减小， 当r ea等于真值时， 四种冰晶粒子的反演结果与真值基本相同， 误差趋于零。模拟研究证明， 准确的先验数据r ea可以提升r e的反演精度。
图4（b）展示了先验值r ea的准确性对计算效率的提升作用， 随着r ea逐渐接近真值， 四种形状的冰晶粒子反演计算的平均迭代次数基本呈现递减趋势， 说明准确的先验数据可以减少最优估计算法的迭代次数， 提高计算效率。
在反演计算过程中， 四种形状的冰晶粒子的平均迭代次数相近， 随着IWC a逐步贴近真值， 呈现递减的趋势， 当IWC a=0.001IWC时， 算法平均迭代次数约为7次， 当IWC a=IWC时， 迭代次数减少为1次。同样证明了准确的先验值对最优估计反演算法的准确性和计算效率的优化作用。
联合上述的反演结果分析可知， 在N T为定值的假设下， 当x a（r ea， IWC a）偏离真值时， 经过多次迭代运算后， 最优估计反演算法仍然可以得到贴近于真值的反演结果， 说明该算法不完全依赖先验值的可靠性， 这是由于观测值（此处为模拟值Ze）已经包含了待反演的r ere和IWC re信息。
3.3 粒子数密度 N T的反演误差分析
前一小节证明了先验数据x a（r ea， IWC a）设定为真值时， 利用最优估计算法反演r ere、 IWC re和N Tre的可靠性， 接下来分析当IWC a偏离真值时， 算法反演得到N Tre的误差变化情况。由公式（3） 和（4）可知， r e与N T的变化毫无关联， IWC与N T成正比， 因此， 假设r e为固定值（50 μm）， 分析先验值IWC a的变化对于N T反演精度的影响。
当先验数据等于真值时， 基于最优估计理论反演分别得到四种冰晶粒子的r ere、 IWC re和N Tre， r ere基本为定值， IWC re和N Tre反演结果如图6所示， 四种冰晶粒子反演结果均与真值基本吻合。
假设r ea等于真值， 分析IWC a不同程度偏离等于真值时， IWC re和N Tre的平均反演误差的变化， 以及反演计算过程中平均迭代次数的变化（图7）。随着IWC a逐步贴近真值， 四种冰晶粒子的平均反演误差和平均迭代次数均呈现递减的趋势。
由图7（a）可知， 子弹花粒子的N Tre反演误差最大， 六棱柱和六角板粒子的N Tre反演误差十分相近， 基本重合， 球形冰晶粒子的N Tre反演误差最小。当IWC a=0.001IWC时， 球形粒子的N Tre平均误差为1.94%， 六棱柱和六角板粒子的N Tre平均误差相似， 分别为2.64%和2.63%， 子弹花粒子的N Tre平均误差为2.71%， 与前两者相差不远， 而球形粒子的N Tre平均误差为1.94%， 当IWC a=IWC时， 球形粒子的N Tre平均误差减少到0.2%左右， 其他三种粒子的N Tre平均误差减少到1%左右。球形冰晶粒子的IWC re反演误差最大， 六棱柱和六角板粒子的反演误差最小， 且基本重合［图7（b）］， 当IWC a=0.001IWC时， 六棱柱和六角板粒子的IWC re误差分别为1.67%和1.66%， 子弹花粒子的IWC re误差为1.74%， 球形粒子的IWC re误差为1.91%。球形粒子的算法平均迭代次数最少， 计算效率最高， 其他三种粒子平均迭代次数相同， 当IWC a=IWC时， 平均迭代次数减少到1次。冰晶粒子的形状在一定程度上会影响算法反演误差和计算效率。
由此说明， 在r e为定值的假设下， 该算法仍旧能不完全依赖先验值， 在迭代初值严重偏离的情况下， 准确反演得到冰云微物理参数r ere、 IWC re和N Tre， 与真值贴近的先验值对反演算法的准确性和计算效率有着明显的提升作用。
4 结论与讨论
根据不同形状冰晶粒子在太赫兹频段（220 GHz）的散射特性， 假设冰晶粒子服从Gama分布， 建立云粒子散射特征数据表LUT， 基于最优估计理论， 利用雷达反射率因子Ze反演冰云粒子有效半径r e、 冰水含量IWC和粒子数密度N T， 模拟研究结果表明， 该算法适用于220 GHz云雷达探测数据的冰云微物理参数反演。
将LUT中的雷达反射率因子Ze作为测量向量， 对应的有效粒子半径r e和冰水含量IWC作为先验数据x a（r ea， IWC a）， 假设N T为定值， 基于最优估计理论反演得到r ere和IWC re。由于Ze包含了待反演的r ere和IWC re信息， 最优估计反演算法可不依赖先验数据， 反演得到贴近于真值的反演结果， 且随着先验数据x a逐步贴近真值， 四种形状的冰晶粒子的平均反演误差和迭代次数基本呈现递减的趋势。然后， 假设r e为固定值， 同样将LUT中的雷达反射率因子Ze作为测量向量， 基于最优估计理论反演得到r ere、 IWC re和N Tre， 且随着IWC a逐步贴近真值， 四种形状的冰晶粒子的平均反演误差和迭代次数逐渐减小， 直到IWC a取真值时， 反演结果与真值基本吻合。
上述研究证明， 基于最优估计理论的反演算法可有效应用于220 GHz球形及非球形冰晶粒子的微物理参数r e、 IWC和N T反演， 该算法不完全依赖先验值x a， 在先验数据严重偏离的情况下， 仍然可以准确反演得到卷云冰晶粒子的微物理参数， 而与真值贴近的x a对反演算法的准确性和计算效率有着显著的提升作用。
本文的首要研究目的是开发可应用于220 GHz冰云观测数据的球形及非球形冰晶粒子的微物理参数反演方法， 为未来地基或星载太赫兹云雷达参数业务化反演算法的研发提供可行的计算方法， 根据冰云实测数据的形状统计结果选择了四种常见的冰晶粒子形状开展模拟研究， 但不能确保与实际冰云没有区别， 因此需要进一步开展观测试验， 根据实测的云中冰晶粒子形状及尺度的真实分布情况， 在更为贴近真实冰云特征的情况下继续开展反演算法改进和验证工作。