引言

云中液态水含量是极为重要的大气物理学参数(雷恒池等，2003)。云在地球能量平衡中的重要作用多来自其中的液态水含量(Turner et al，2007)，其反馈作用对全球气候变化产生十分重要的影响(Somerville and Remer, 1985；Stephens and Greenwald, 1991)。

微波辐射计是测量云中液态水的一种有效手段且应用广泛(崔雅琴等，2019；许皓琳等，2020；李金辉等；2022； 辛悦等，2023)。与地基微波辐射计相比，机载微波辐射计不受空间限制，可通过改变飞机位置探测云中不同区域的液态水，有利于探测人工影响降水作业潜力区——过冷水区。利用机载微波辐射计探测云中液态水开始于20世纪80年代。Warner and Drake(1988)研制了9.5 mm波长双向天线机载微波辐射计，在反演算法上采用层析技术原理，将目标云系离散化，根据辐射传输原理求得云中液态水含量的二维分布(Drake and Warner, 1998)。2000年初，雷恒池等(2003)在我国首次研制成功单波段单天线对空机载微波辐射计，其工作频率为31.65 GHz(波长为9.5 mm)。江芳等(2004)针对该辐射计获得的亮温数据采用统计方法进行反演，获取垂直累积液态水。金德镇等(2004)、梁谷等(2007)采用该反演算法对层状云进行了液态水观测。周珺等(2008)在上述统计方法基础上，对云液态水廓线产生方法进行改进，采用一维层状云模式产生云样本进行垂直累积液态水反演计算，探测误差估计表明其反演精度在所有高度上都优于原方法，之后又采用层析法建立了机载微波辐射计探测云中液态水二维分布反演算法(周珺等，2010)。

2004年底，美国Prosensing公司研制了中心工作频率在183.31 GHz的双边四通道对空微波辐射计GVR(G-band water vapor radiometer)，四个通道分别为183.31(±1，±3，±7，±14) GHz(Pazmany，2006；2007)。183.31 GHz位于水汽强吸收带，具有更高的灵敏度，适于探测较低含量水汽，同时，该吸收带对液态水测量也极为灵敏(周秀骥，1982；Pazmany，2007)，其测量结果受云中固态水影响很小，但对液态水非常敏感(Wang et al，1998)，因此183 GHz机载微波辐射计在测量低含量水汽和液态水特别是过冷水含量上具有优势。Pazmany(2007)针对GVR采用神经网络算法对积分水汽含量(PWV)和液态水路径(LWP)进行反演，输入变量为四个通道亮温及大气温度，训练数据选取美国中纬度冬季和北极地区的探空资料。2016年，天津市人工影响办公室引进该仪器并采用其自带算法(Pazmany算法)进行水汽和云中液态水观测，作者在数据质量控制基础上对该算法在本地的适用性进行研究(王婉等，2018a；2018b；2021)，结果表明与其他类型机载微波辐射计反演精度及GVR在北极地区地基观测反演精度相比，GVR自带反演算法在天津地区的反演精度相对偏差。分析认为造成这种误差的原因主要是未考虑飞机高度变化，云模型构建过于简单等。为更好使用该仪器，本文将建立新的积分水汽含量和云中液态水反演算法，通过数值模拟、飞机外场观测试验及仪器观测误差传递方法验证新算法及GVR自带算法反演精度并进行比较。

1 仪器及数据资料

GVR安装于飞机机翼下方，通过天线向上垂直观测，接收飞机上方云水信号。GVR天线是一口面直径10 cm、90°的抛物面金属镜，置于辐射计前端，伸出机翼外侧，天线正上方有一天线罩，来自辐射计上方的辐射能量透过天线罩被天线接收并反射进入辐射计接收机，从而测得四个通道辐射能量。由四个通道测量的亮温可以反演得到从飞机所在高度直到大气上界这一柱体中的积分水汽含量和积分液态云水含量(液态水路径)。

本文利用历史探空数据，采用BP神经网络算法建立数学模型来反演飞机上方垂直气柱内的积分水汽含量和液态水路径。文章用到两类数据，历史探空数据和小时雨量数据。历史探空数据分为两部分：一部分用作建立数学模型，选取北京探空站2009—2018年08时和20时(北京时，下同)两个时次的探空资料；另一部分用作反演精度数值模拟，取北京站2019年08时和20时两个时次的探空资料。由于反演垂直累积液态水需要在探空资料中筛选有云无降水样本，文章选用北京探空站2009—2019年每日08时和20时两个时次的小时雨量资料，剔除雨天探空数据。

2 研究方案

GVR自带算法采用的是Pazmany(2007)建立的反演算法，新算法将在自带算法基础上进行改进并将两种算法的反演精度进行比较。表 1为新算法与GVR自带算法反演方案对比。

2.1 建立云天样本集和构建云模型

在历史探空数据中筛选有云无降水样本，建立云天样本集。首先在所有探空资料中剔除数据缺失资料，以相对湿度85%作为有云判据，只要存在某高度层相对湿度大于85%则视为有云探空资料，然后在有云探空中剔除雨天探空数据，由于08时和20时两个时次的探空施放时间为07：15和19：15，因此选择08时和20时的小时雨量资料，只要该时次存在降水则剔除相应时次探空数据。本文是针对机载微波辐射计来反演垂直累积液态水，因此高度是不断变化的，为获取更多样本，将每个探空0~9 km随机取8个高度层作为样本初始高度(飞机观测高度)，这样每个探空资料可获取8个探空样本，从而建立云天样本集。

针对每个云天样本构建云模型。云的模型与实际云的差别会给反演结果带来较大误差，因此准确构建云模型是影响反演精度的关键因素。国外常用的两种云模型为Decker模型(Decker et al，1978)和Salonen模型(Salonen and Uppala, 1991)，国内许多云液态水廓线构建方法类似于Decker模型(黄润恒和邹寿祥，1987；刘亚亚等，2010；王旭等，2016)，本文参照该模型，同时考虑温度影响，参照Mattioli et al(2009)进行云模型构建。具体做法为：针对每个探空廓线样本，当相对湿度＜85%时，液态水密度取0 g·m^-3；当相对湿度 > 95%、温度 > 0℃时，液态水密度取0.5 g·m^-3；当相对湿度 > 95%、温度≤-30℃时，液态水密度取0 g·m^-3；当相对湿度 > 95%、温度处于-30~0℃时，液态水密度计算公式如下：

$ \rho_{\mathrm{L}}=0.5\left[1-(T / 30)^4\right] $

(1)

当相对湿度处于85%~95%、温度 > 0℃时，液态水密度计算公式如下：

$ \rho_{\mathrm{L}}=-4.25+0.05 H_{\mathrm{u}} $

(2)

当相对湿度处于85%~95%、温度 < -30℃时，液态水密度取值为0 g·m^-3；

当相对湿度处于85%~95%、温度在-30~0℃时，液态水密度计算公式如下：

$ \rho_{\mathrm{L}}=\left(-4.25+0.05 H_{\mathrm{u}}\right) \cdot\left[1-(T / 30)^4\right] $

(3)

液态水路径(LWP)计算公式如下：

$ \mathrm{LWP}=\int_{z_0}^{\infty} \rho_L \mathrm{~d} z $

(4)

式中：ρ_L (单位：g·m^-3)表示液态水密度，T(单位：K)表示温度，H_u(单位：%)表示相对湿度，z₀(单位：m)表示飞机观测高度。

2.2 积分水汽含量计算

利用探空资料计算积分水汽含量(PWV)的方法有多种，本文将采用如下公式(盛裴轩等，2003)：

$ \mathrm{PWV}=\int_{z_0}^{\infty} \rho_{\mathrm{H}_2 \mathrm{O}} \mathrm{d} z $

(5)

式中：$\rho_{\mathrm{H}_2 \mathrm{O}}=\frac{\varepsilon e_{\mathrm{s}} \mathrm{D}}{R_{\mathrm{d}} T}, \rho_{\mathrm{H}_2 \mathrm{O}}$ (单位：kg·m^-3)表示水汽密度，z₀表示飞机飞行高度；ε=0.622，e_s表示水面饱和水汽压(单位：Pa)，D表示相对湿度，R_d表示干空气比气体常数，取值287.05 J·kg^-1·K^-1，T(单位：K)表示某高度对应的温度。

2.3 大气下行辐射传输方程

GVR在飞机飞行高度进行观测，接收来自其上方的微波辐射，不考虑降雨时，其微波辐射传输方程可采用如下形式(江芳等，2004)：

$ \begin{aligned} T_{\mathrm{b}}\left(z_0, f\right)= & T_{\mathrm{b}}(\infty) \mathrm{e}^{-\int_{z_0}^{\infty} a(z, f) \mathrm{~d} z}+\int_{z_0}^{\infty} a(z, f) T(z) \cdot \\ & \mathrm{e}^{-\int_{z_0}^z a\left(z^{\prime}, f\right) \mathrm{d} z^{\prime}} \mathrm{d} z \end{aligned} $

(6)

式中：T_b(单位：K)表示辐射计接收到的大气向下的辐射亮温，T_b(∞)(单位：K)表示宇宙背景亮温，一般取常数2.7 K，z₀(单位：m)表示飞机飞行高度，z(单位：m)表示飞机至宇宙背景的某一高度，f(单位：GHz)表示辐射计工作频率，T(z)(单位：K)表示高度z处的大气温度，z′(单位：m)表示飞机至高度z之间的某一高度，a(单位：K·m^-1)表示大气吸收系数，一般为氧气(a_O2)、水汽(a_H2O)和云中液态水吸收系数(a_cloud)之和：

$ a=a_{\mathrm{O}_2}+a_{\mathrm{H}_2 \mathrm{O}}+a_{\text {cloud }} $

(7)

氧气、水汽和液态水吸收系数参照Ulaby et al(1981)公式。

氧气的微波吸收谱线是由许多分布在50~70 GHz频率范围上的吸收谱线和在118.75 GHz的附加谱线所组成，氧气的吸收系数为频率、大气压力和大气温度的函数：

$ a_{\mathrm{O}_2}=1.61 \times 10^{-2} f^2\left(\frac{p}{1013}\right)\left(\frac{300}{T}\right) F^{\prime} $

(8)

式中：F′是谱线线型函数，f(单位: GHz)为频率，p(单位: hPa)为大气压力，T(单位: K)为大气温度。

水汽吸收系数为频率、大气压力和大气温度的函数，表达式为：

$ a_{\mathrm{H}_2 \mathrm{O}}=a_{\mathrm{H}_2 \mathrm{O}}^{\prime}+\Delta a $

(9)

$ \begin{aligned} a_{\mathrm{H}_2 \mathrm{O}}^{\prime}= & 2 f^2 \rho_\nu\left(\frac{300}{T}\right)^{5 / 2} \sum\limits_{i=1}^{10} A_i \mathrm{e}^{-\varepsilon_i / T} \cdot \\ & {\left[\frac{\gamma_i}{\left(f_i^2-f^2\right)^2+4 f^2 \gamma_i^2}\right] } \end{aligned} $

(10)

γ_i为线宽函数，计算公式为

$ \gamma_i=\gamma_{i 0}\left(\frac{p}{1013}\right)\left(\frac{300}{T}\right)^x\left[1+10^{-2} a_i \frac{\rho_v T}{p}\right] $

(11)

Δa为100~300 GHz频率范围内使用的订正项，计算公式为：

$ \Delta a=4.69 \times 10^{-6} \rho_v\left(\frac{300}{T}\right)^{2.1}\left(\frac{p}{1000}\right) f^2 $

(12)

式中: ρ_v(单位：g·m^-3)为水汽密度，γ_i0，a_i，x, A_i，ε_i和f_i均为给出参数值。

液态水的吸收系数与滴谱分布无关，是频率、云滴温度和云中含水量的函数：

$ a_{\text {cloud }}=a_1 \cdot m $

(13)

式中：m(单位：g·m^-3)为云水含量，a_l表示液态水的质量吸收系数，与波长和云滴温度有关，其表达式为：

$ a_1=0.434 \frac{6 \pi}{\lambda_0} \operatorname{Im}\{-K\} $

(14)

$ K=\frac{m^2-1}{m^2+2} $

(15)

式中：λ₀(单位：cm)为波长，m为水的复折射指数，可表示为：

$ m=n-i k $

(16)

式中：n，k与波长和云滴温度有关，计算采用Liebe et al(1991)建立的水复折射指数经验公式。

GVR是双边四通道微波辐射计，四个通道的带宽分别为0.5，1.0，1.4和2.0 GHz，在计算通道亮温时需要根据式(6)和频率响应函数计算加权后的亮温，由于缺少频率响应函数本文参照Zuidema et al(2012)做法，采用各通道中心频率来模拟通道亮温，同时假设每个通道的两个边带对测量的贡献是一致的，通道亮温取两者平均。

3 反演误差分析 3.1 积分水汽含量反演误差分析

本文分别采用数值模拟检验、外场观测试验和观测误差传递方法对积分水汽含量进行误差分析。

数值模拟检验指选择一定探空资料作为验证样本，以探空资料计算的物理量作为真实值，利用探空资料模拟仿真通道亮温作为反演方程的输入参数计算相应物理量，以此作为反演值，通过物理量的反演值和真实值对比获取反演误差。本文以相关系数(R)、统计均方根绝对误差(E)、统计均方根相对误差(S)三个统计变量来体现反演误差，其计算方法为：利用探空资料建立云天样本集，采用式(5)计算每个样本的积分水汽含量作为真实值，记作V_Zi，i=1，2，…，n，将所有样本真实值的平均记为V_Z，采用式(6)计算每个样本辐射亮温并读取对应高度及温度，通过反演算法计算积分水汽含量作为反演值，记作V_Ri，i=1，2，…，n，所有样本的积分水汽含量反演值的平均记为V_R，则相关系数、统计均方根绝对误差、统计均方根相对误差计算公式如下：

$ R=\frac{\sum\limits_{i=1}^n\left(V_{\mathrm{R} i}-\overline{V_{\mathrm{R}}}\right)\left(V_{\mathrm{Z} i}-\overline{V_{\mathrm{Z}}}\right)}{\sqrt{\sum\limits_{i=1}^n\left(V_{\mathrm{R} i}-\overline{V_{\mathrm{R}}}\right)^2 \cdot \sum\limits_{i=1}^n\left(V_{\mathrm{Zi}}-\overline{V_{\mathrm{Z}}}\right)^2}} $

(17)

$ E=\sqrt{\sum\limits_{i=1}^n\left(V_{\mathrm{Z} i}-V_{\mathrm{R} i}\right)^2 / n} $

(18)

$ S=\frac{E}{\left(\sum\limits_{i=1}^n V_{\mathrm{Zi}_i} / n\right)} $

(19)

取北京2019年探空资料作为验证样本，计算积分水汽含量真实值，并分别采用新反演算法(GVR-NEW)和GVR自带算法(GVR-OWN)计算每个积分水汽含量真实值对应的反演值。利用式(17)计算积分水汽含量真实值和反演值的相关系数并将对应值标于散点图中，结果见图 1。GVR自带算法计算的反演值与真实值相关系数为0.9929，新算法为0.9988。显然新算法计算的积分水汽含量反演值与真实值相关性高于GVR自带算法。图 1a为GVR自带算法计算的反演值与真实值散点图，图 1b为新算法计算的反演值与真实值散点图，图 1中斜线表示反演期望值，斜线及两侧的散点表示反演值偏离期望值(真实值)的程度，从图中可以看出，与新算法相比，GVR自带算法计算的积分水汽含量反演值偏离真实值较大且普遍低于真实值，两种算法计算的反演值偏离期望值的程度都随着水汽含量增大而增大，说明积分水汽含量越大反演误差越大。

利用北京2019年探空资料分别计算GVR自带算法和新算法反演的积分水汽含量统计均方根绝对误差和相对误差，其计算结果随高度分布见图 2。图 2中蓝色曲线表示GVR自带算法的反演误差，黑色曲线表示新算法反演误差，由图可见，新算法和GVR自带算法反演的积分水汽含量统计均方根绝对误差分别为0.05~1.30 mm、0.2~3.0 mm，相对误差分别为1%~10%、4%~65%。Cadeddu et al(2007)利用GVR在北极地区进行地基观测获得的积分水汽含量反演结果相对误差低于5%，新算法在6 km以下计算的积分水汽含量相对误差最小值为1%，整层高度上约75%相对误差低于5%, 而自带算法仅有不到1%的相对误差低于5%，因此与GVR自带算法相比，新算法反演的积分水汽含量误差更小，反演效果更好。

上述分析表明新算法反演效果优于GVR自带算法，从表 1可看到新算法在自带算法基础上做了三处改进。为验证每一处改进对改善反演效果的贡献，本文在自带算法基础上每次做一处改进，可对应得到三个反演方案，这三个方案分别命名为M1(改变云模型)、M2(增加高度数据作为输入参数)和M3(补齐探空资料到30 km高度)。取北京2019年探空资料作为反演结果验证数据，为便于比较，直接计算所有验证样本的统计均方根绝对误差，计算结果显示GVR自带算法、方案M1、M2和M3反演的PWV统计均方根绝对误差分别为1.29、0.69、0.91、1.12 mm，方案M1、M2和M3对反演效果改进的贡献分别为52%、33%和15%，可见云模型对反演效果的影响最大。

利用数值模拟检验方法验证积分水汽含量反演精度时，由于输入亮温取自探空资料仿真值而非仪器实际探测结果，因此其计算误差仅能反映反演算法本身造成的误差，其中包括拟合反演方程时用到的历史资料背景场差异等因素。除反演误差外，仪器误差也会影响最后的水汽观测结果。本文结合2016年11月20日飞行探测个例来讨论仪器误差对水汽观测结果造成的影响。GVR定标误差低于1 K(Pazmany，2007)，本文以1 K作为仪器误差，分别对新算法和自带算法进行误差传递计算(经反演公式传递后造成的误差)，结果在图 3中分别以黑色和蓝色虚线表示。图 3中黑色和蓝色实线分别表示采用新算法和自带算法计算的本次积分水汽含量反演值，由图可见，与新算法相比，自带算法反演值偏低，这与图 1中自带算法模拟计算的反演值偏低是一致的。从图中虚线可以看到，两种算法经反演公式传递后造成的误差值都随着积分水汽含量值增大而增大，在3000 m飞行高度，自带算法与本文算法的传递误差分别约为0.06 mm和0.05 mm, 随水汽探测值增大，自带算法传递误差增大幅度明显高于新算法。

3.2 液态水路径反演误差分析

由于液态水无法直接由探空计算得到，其值与构建的云模型相关，特别是GVR自带算法对于每个探空样本都添加了数条液态水廓线，因此很难利用数值模拟检验法验证其反演精度。利用外场观测中不同仪器的探测结果进行对比分析是验证反演精度的最直观方式，但由于缺少其他有效的液态水路径探测手段，无法获取其真值，为外场观测验证带来一定困难。当飞机GVR上方无云存在时，可确定其上方的液态水路径真值为0 mm，这为外场观测对比提供了可能。要确认飞机上方无云需要必要的探测设备，由于目前携带GVR的增雨飞机缺少这种设备，故采用安装在飞机上的视频拍摄装置寻找无云时段，同时通过DMT粒子测量系统确认飞机是否处于云外。通过多个观测个例分析，选择2021年6月30日个例。2021年6月30日10:18—15:53，天津增雨飞机空中国王B300进行了一次探测飞行，根据安装在飞机上的视频显示，本次飞行有一段时间处于云层上方且上方完全无云。

图 4绿色曲线为DMT粒子测量系统的云、气溶胶和降水粒子组合探头CIP测量的粒子数浓度随时间变化，图中可见，12:30—12:54粒子数浓度基本为0个·cm^-3，这和机上人员观测是一致的，飞机该时段处于云外飞行。选取机载微波辐射计GVR 12:30—13:10观测数据，按照GVR亮温数据中RFI识别和订正方法(王婉等，2018a; 2018b)进行质量控制，然后分别采用GVR自带算法和本研究算法反演计算液态水路径(LWP)，红色曲线表示自带算法反演结果，蓝色曲线表示新算法反演结果(图 4)。由图可见，两种反演结果在12:39—12:48时段内保持一致，基本趋于0 mm，这和机上观测及CIP探头观测结果都是一致的，可认为该时段内飞机在云外飞行且上方无云层，液态水路径值应为0 mm。图 5为该时段内的GVR自带算法和本研究算法反演结果，从图中可以看到，该时段内新算法获取的液态水路径反演值全部为0 mm，自带算法获取的反演值范围为0~0.006 mm，平均值约为0.001 mm，显然，新算法计算的液态水路径精度更高。12:48—12:55时段内，GVR自带算法反演出一定的垂直累积液态水，变化范围在0.001~0.020 mm，但新算法反演值和CIP粒子探头测量的粒子数浓度都趋于0个·cm^-3，此时段飞机飞行高度位于4000~4700 m，不排除飞机上方存在一定云层，但据机上视频显示，此时段飞机上方并无云层，因此该时段内新算法反演结果更为可信。对比本次飞行的液态水路径反演结果，新算法反演值普遍低于自带算法，当不存在液态水时，新算法反演精度高于自带算法，但当存在液态水时，由于暂时缺少对比观测，尚无法定量计算两种算法的反演误差。

4 结论

(1) 本文在以往研究基础上，建立新的积分水汽含量和云中液态水反演算法，与GVR自带算法相比，新算法增加了观测高度作为输入变量，考虑了温度对云水相态的影响建立新的云模型构建云液态水廓线，对历史探空资料进行补全，增加微波辐射传输方程积分高度至30 km。

(2) 利用数值模拟检验、外场观测试验和观测误差传递分析，对新算法和GVR自带算法反演的积分水汽含量反演误差进行验证。新算法和自带算法模拟计算的反演值与探空计算值相关系数分别为0.9988、0.9929，新算法计算的积分水汽含量反演值与真实值相关性高于GVR自带算法。两种算法计算的反演值与真实值散点图表明，与新算法相比，GVR自带算法计算的积分水汽含量反演值偏离真实值较大且普遍低于真实值，两种算法计算的反演值偏离期望值的程度都随着水汽含量增大而增大。以2016年11月20日飞行探测个例为例，讨论仪器误差对水汽观测结果造成的影响，结果表明两种算法经反演公式传递后造成的误差值都随着积分水汽含量值增大而增大，在3000 m飞行高度，自带算法与新算法的传递误差分别约为0.06 mm和0.05 mm，随水汽探测值增大，自带算法传递误差增大幅度明显高于新算法。为验证每一处改进对改善反演效果的贡献，本文在自带算法基础上每次做一处改进，结果表明云模型、增加观测高度作为输入变量和补齐探空资料到30 km高度对反演效果改进的贡献分别为52%、33%和15%。

(3) 选择2021年6月30日飞行个例，对两种算法反演的液态水路径反演误差进行验证，结果表明：在飞机上空无云时段，新算法和自带算法反演的液态水路径范围分别为0 mm、0~0.006 mm。误差分析表明，新算法可提升本地积分水汽含量和液态水路径反演精度。