海雾是发生在海上大气边界层内而导致近海面大气水平能见度低于1km的天气现象, 它严重影响海上作业及沿海地区的交通安全。黄海是我国近海海雾高发区之一(王彬华, 1983; Gao et al, 2007; 张苏平等, 2008), 且黄海海雾具有突发性特点。近年来, 随着计算能力的迅速提升, 数值模拟已成为研究黄海海雾的有力工具(傅刚等, 2016; 史得道等, 2016); 海雾数值预报尤其是短时临近预报, 需求迫切。

已有研究表明, 海雾数值预报对初始场很敏感(Nicholls, 1984; Findlater et al, 1989; Ballard et al, 1991; Koračin et al, 2001, 2005a, b; Lewis et al, 2003; Fu et al, 2006; Gao et al, 2007; 高山红等, 2010)。因此, 利用数据同化手段改进海雾数值预报中的初始场质量是一个关键问题。高山红等(2010)基于WRF(Weather Research and Forecasting)模式及其3DVAR(3-Dimen sional Variational)模块, 设计了一个循环3DVAR同化方案, 通过扩大同化窗而引入更多的常规观测, 显著改进了海雾模拟效果。在此基础上, 李冉等(2012)进行了卫星辐射数据的直接同化, 王静菊等(2017)同化AIRS(Atmospheric Infra-Red Sounder)反演温度廓线数据, 改善了海上大气边界层(Marine Atmospheric Boundary Layer, MABL)的温度层结。此外, 同化卫星反演海面风(刘宇迪等, 2011)与雷达径向风数据(王永明等, 2016)能够不同程度地改善控制海雾发展的MABL风场结构。

相比MABL的温度层结与风场结构, 初始场的湿度状况对黄海海雾短时临近预报更为重要。针对海上湿度观测严重不足的问题, Wang等(2014)提出了一种利用静止卫星MTSAT(Multifunctional Transport Satellite)数据反演湿度的同化方案, 将黄海海雾短时临近预报的ETS(Equitable Threat Score)评分提高了15%−20%。该方案的思路为:在同化窗(譬如预报前6h)内海雾已经发生的前提下, 利用静止卫星多通道数据反演海雾的三维信息(水平雾区与厚度), 假定雾体内部相对湿度(Relative Humidity, RH)为100%, 借助循环3DVAR手段同化RH; 在3DVAR更新过程中, RH被转变为水汽混合比进行同化。然而, 模式微物理方案本身的不完善, 会使得同化RH后的初始场, 未必能出现准确的海雾雾区。

当初始场存在海雾时, 雾顶的长波辐射将会导致雾顶冷却, 产生自上而下的湍流与加强雾顶卷夹作用(Oliver et al, 1978; Findlater et al, 1989; 杨悦等, 2016; Yang et al, 2019), 这会影响海雾的后续发展。因此, 直接同化海雾雾体内的云水含量信息, 将观测到的海雾在模拟初始场中尽可能地准确再现出来, 无疑将会改进海雾的短时临近预报。受技术手段限制, 目前尚无法利用静止卫星反演得到海雾内部的云水含量。然而, 星载微波辐射计对在海洋上空大气水含量的观测和反演工作已有近30年的历史(Mears et al, 2015), 整层大气云水含量(即云水路径, Cloud Water Path)的反演算法已经相当成熟(王颖等, 2016)。

因此, 本文尝试同化云水路径, 试图在初始场中再现已经发生的海雾, 从而改进海雾短时临近预报, 类似的工作还未曾见过。由于云水路径数据来源于极轨卫星, 其时空分布范围有限, 为了将此有限的观测信息合理扩散从而提高同化效果, 拟采用流依赖的同化手段——EnKF(Ensemble Kalman Filter; Evensen et al, 1996)。Gao等(2018)的研究已经表明, EnKF在海雾预报中的同化效果显著优于3DVAR。本文将基于GSI(Gridpoint Statistical Interpolation system)- EnKF同化系统(Shao et al, 2016)和WRF模式, 开展云水路径数据同化工作, 探究其对黄海海雾短时临近数值预报的改进效果。

WRF模式的背景场和时变侧边界条件来源于FNL(Final Analysis of NCEP, National Centers for Environmental Prediction)全球再分析数据, 水平分辨率1°×1°, 垂直共有26个气压层, 下载网址为https://rda.ucar.edu/datasets/ds083.2。底边界的海表面温度(Sea Surface Temperature, SST)来自NEAR- GOOS(North- East Asian Regional Global Ocean Observing System)发布的逐日融合产品, 分辨率0.25°×0.25°, 下载网址为http://ds.data.jma.go.jp/gmd/goos/data 。

用于同化的观测数据包括:常规地面和探空数据, 下载网址为http://222.195.136.24/gts_obs/gts_data; 卫星搭载的AMSU-A/B(Advanced Microwave Sound­ ing Unit A/B), HIRS-3/4(High Resolution Infrared Radiation Sounder 3/4)和MHS(Microwave Humidity Sounder)探测仪获取的辐射数据, 下载网址为https://rda.ucar.edu/datasets/ds735.0)。

云水路径数据(记为Cw)来自DMSP(Defense Meteorological Satellite Program)系列卫星搭载的SSMI(Special Sensor Microwave Imager)辐射计数据所反演的网格化数据, 分辨率为0.25°×0.25°, 下载网址为ftp://ftp.ssmi.com/ssmi。由于Cw来自极轨卫星, 黄海海域每天只有2次数据。图 1给出了2011年5月16日20:48 UTC黄海上空Cw的水平分布。

2011年5月16−17日, 黄海发生了一次大范围的海雾过程。图 2展示了16日21 UTC−17日09 UTC的海雾分布; 其中, 图 2a为利用Wang等(2014)方法得到夜间海雾反演结果, 图 2b−图 2d为MTSAT可见光图片。可见, 海雾最初形成于黄海中部, 受偏南风影响(见图 4中的海面风场)而向北发展; 17日上午, 海雾覆盖了长江入海口以北、成山头以南的大片海域, 且雾区的边缘已接触到山东半岛和江苏省沿岸地区; 17日下午, 受东移的陆地低压影响, 海上风速增大, 海雾自西向东逐渐消散。对比图 1和图 2可知, 在黄海中部的雾区范围内观测到了Cw, 体现了海雾的云水含量信息; 而黄海南部也存在大片的Cw分布, 与高云分布相对应(见图 2b、图 2c的下部)。

背景场、观测数据和同化分析场分别为xb, y和xa, 模式空间转换到观测空间由H实现, 同化更新过程可以表达为

其中, K为增益矩阵, 它根据观测和背景场误差的大小, 将观测增量y–Hxb转化为分析增量xa–xb。设观测误差为R, 背景误差为Pb, K由公式(2)决定:

在EnKF中, 背景场xb随机扰动后, 产生m个集合成员: 。WRF模式分别将每个成员积分至同化时刻, 统计集合成员预报结果的协方差作为背景误差。实际计算中, 为了减小计算代价, 集合同化不会直接计算庞大的Pb矩阵, 而是计算转换后的PbHT和HPbHT:

将每个集合成员分解为集合平均与扰动, 则用于同化更新的公式(1)被分解为(Whitaker et al, 2002):

受样本误差的影响, 随着预报的进行, 将会出现对背景误差的系统性低估, 使同化结果逼近背景场, 观测信息难以吸收, 导致滤波发散。样本误差来自有限集合成员统计特征的随机变化, 无法从根本上消除。为了削弱其影响, 使用Whitaker等(2012)提出的集合膨胀技术来增大分析场的离散度:

式中, 和分别为背景场和分析场在每个格点上的标准差。β参数控制膨胀程度:若, 则分析场不被膨胀; 若, 则分析场集合成员的离散度被膨胀至与背景场相同。

样本误差还会导致协方差估算中出现虚假的相关性。为了消除它, 将观测数据的同化限制在一定的空间范围内, 需要对协方差进行局地化处理(Houtekamer et al, 2001)。本文使用Gaspari等(1999)提出的五阶多项式函数来实现水平局地化, 将背景误差乘以一个系数γ :

代表对每个元素的舒尔积, γ是背景场格点与观测点的水平距离z与局地化尺度L的函数:

云水路径Cw不是模式预报变量, 在被同化之前必须建立它与背景场之间的相关关系。在WRF模式中, 云水混合比(记为qc)被用于表征云中液态水含量。当控制变量只包括温度T, 相对湿度RH, 表面气压Ps和风速U, V时, 背景场可以表示为

以公式(13)为基础, 增加qc的统计和更新, 背景场变为:

将Cw的观测值记为, 由背景场得到的Cw记为。某一水平点的由垂直积分得到:

式中, 为液态水密度, z1和zn分别为模式底层和顶层的高度, 为高度z上的空气密度; 根据理想气体状态方程, 可以由气压和温度计算:

式中, R为比气体常数。公式(15)与公式(16)将背景场变量转化为Cw, 起到了转换矩阵H的作用, 记为Hcw。

EnKF顺序同化每个观测数据, 将同化一个观测产生的分析场作为同化下一观测的背景场。于是, y只有一个自由度, 矩阵PbHT和HPbHT分别退化为矢量和标量。对于某一云水观测, 根据公式(15)与公式(16)计算每个集合成员在观测点的云水路径 (i=1, 2, …, m), 及其平均值, 观测增量即为, 而背景误差可以表示为:

值得注意的是, 公式(17)计算了与所有控制变量之间的协方差。根据的计算公式, 它无疑与qc关系最为密切, 同化反演云水路径将直接影响背景场的qc分布。另一方面, 还包含着水物质相变(例如水汽凝结, 云水蒸发等过程)与天气形势和温湿分布的相关性, 使对qc以外的其他变量产生影响。

图 3是某一时刻的EnKF同化示意图, 可具体分为3步:

(1) 初始背景场随机扰动, 生成m个集合成员, 并利用WRF模式将它们分别积分到同化时刻, 得到背景场集合体: (i=1, 2, …, m)。

(2) 同化时, 公式(15)与公式(16)计算 (i=1, 2, …, m), 公式(17)与公式(18)计算背景误差与, 进行局地化, 从而产生增益矩阵K, 更新集合成员, 进行协方差膨胀, 得到 (i=1, 2, …, m)。

(3) 同化非的其他观测时, 利用式(3)与公式(4)计算背景误差, 进行局地化, 计算增益矩阵; 利用公式(5)−公式(9)更新集合成员, 并进行协方差膨胀。

本文采用WRF模式(V3.8.1)进行海雾的数值预报, 模式区域为双重嵌套(图 4)。分辨率与物理选项等其它具体设置见表 1。陆雪等(2014)检验了不同边界层方案和云微物理方案模拟海雾的表现, 发现YSU(Yonsei University)方案与PLin(Purdue Lin)方案的组合最佳, 本文沿用这一组合。FNL数据与NEAR-GOOS SST为模式提供初边值条件。

借助GSI(V3.5)及EnKF(V1.1)模块进行数据同化。图 5给出了EnKF的循环同化过程, 同化窗为5月16日18 UTC—17日00 UTC, 间隔为3h, 循环同化3次。obs代表非的观测数据(常规地面和探空数据, AMSU-A/B、HIRS-3/4和MHS辐射数据); 由于数据时间为20:48 UTC, 它仅在第2次循环被同化。在同化起点(18 UTC), 利用随机扰动法(Barker et al, 2004; Wang et al, 2008a, b)生成初始背景场集合体 (i=1, 2, …, 40)。EnKF同化中, 膨胀参数β与局地化尺度L是2个重要参数, 它们会影响海雾的同化效果(Gao et al, 2018)。

设计并实施了2个数值试验: Exp-A仅同化常规观测和卫星辐射数据, Exp-B在Exp-A的基础上同化了。试验的详细设置见表 2。

雾区是海雾预报主要关注的要素。使用Wang等(2014)中所用的方法, 基于MTSAT红外和可见光数据反演出海雾观测雾区; 而模拟雾区依据模式预报的云水含量来诊断获取。采用高山红等(2010)提出的“鸟瞰”方式确定模拟雾区, 对于海上的某一水平模拟点, 在其垂直层中自上而下找到qc ≥0.016g/kg (此值对应约1km的大气水平能见度)所对应的高度, 作为雾顶; Zhou等(2010)指出, 海雾的雾顶高度极少超过400m。因此, 若上述雾顶≤400m, 则认为该水平点属于雾区。

图 6给出了几个时刻的观测和模拟雾区的对比。与观测(图 6a—6d)相比, Exp-A的雾区始终偏小(图 6e—6h): 17日01 UTC, 雾区大致呈三角形分布, 其西侧和北侧边缘距离陆地较远, 随后雾区增大, 但始终没有接触山东半岛; 从17日06 UTC开始, 雾区西侧逐渐消散, 与观测雾区的变化趋势较为一致, 但雾区消散较快。与Exp-A相比, Exp-B的雾区明显改善(图 6i—6l): 17日01 UTC, 雾区的西部和北部明显增大, 距离黄海西岸更近; 随后其分别西向与东北向扩展, 06 UTC时雾区已经到达山东半岛, 雾区整体的形状和大小与此时的观测十分吻合(比较图 6c、图 6g、图 6k; 比较图 6d、图 6h、图 6l)。

为了定量检验预报结果, 利用点对点的比较方法统计模拟雾区的正报率(Probability Of Detection, POD), 误报率(False Alarm Ratio, FAR), 雾区偏差Bias和公正预兆评分ETS (Doswell III et al, 1989; Zhou et al, 2010; Wang et al, 2014), 结果见表 3。Exp-A的POD与ETS分别为0.513与0.361, 而它们在Exp-B中分别提高到0.617和0.445;相比于Exp-A, Exp-B的综合评分ETS改进率超过20%。另外, Exp-B的Bias和FAR的改进率分别为31.5%和0.7%, 说明同化在增大正报雾区的同时, 并未导致虚假雾区的增长。

为了进一步评估海面气温与湿度预报结果, 利用船舶观测(位置如图 4, 分别记为S1, S2, …, S7)进行了检验。图 7给出了偏差Bias统计结果。总体上看, 2个试验所预报的海面气温偏低(图 7 a), 水汽偏少(图 7 b)。在船舶S2, S5和S6, 同化后的偏差几乎没有改变; 而在船舶S1, S3和S7, Exp-B较Exp-A的偏差明显减小:气温升高, 水汽增加。其中, 在位于观测雾区中的S1站点处(见图 4的站点位置与图 6的观测雾区), Exp-A和Exp-B的温度偏差分别为–1.41℃和–1.04℃, 水汽偏差分别为–0.83g/kg和–0.62g/kg; Exp-B有了很大改进, 特别是水汽的增加有利于雾区的增长。7个观测点的检验结果平均而言, Exp-A和Exp-B预报气温的均方根误差(RMSE)分别为1.38℃和1.25℃, 水汽混合比的RMSE分别为0.56g/kg和0.47g/kg, Exp-B的气温与水汽分别改进了7.4%和16.1%。这说明, 同化改善了海面气温与水汽的分布, 特别是后者。

预报效果检验表明, Exp-B的海雾雾区明显优于Exp-A, 这应归功于其初始场的改进。图 8给出了在WRF模式底层, 它们初始场之间差异(Exp-B–Exp-A)的水平分布。图 8中粗黑线与绿线分别为Exp-A与Exp-B初始场中雾区轮廓; 可见在预报起点, Exp-B雾区就已经大于Exp-A, 其云水含量(图 8a)也清楚显示如此。图 8b和8c中温度与水汽混合比的差异分布可以给出答案:与Exp-A相比, Exp-B在雾区周围温度显著降低(最大降幅接近1K), 水汽混合比总体偏大(雾区北侧增大了0.2—0.3g/kg); 温度降低促使水汽凝结, 使Exp-B雾区在初始时刻就生成范围较大的雾区(对比图 8a—8c中黑、绿粗线所示的范围)。图 8d显示, 相比于Exp-A的海面风场, Exp-B存在一个雾区附近呈辐散特征的风场增量(最大风速增量近1m/s), 这个增量是造成Exp-B中海雾逐渐向西和东北方向扩展的主要动力因素(对比图 6中Exp-A与Exp-B的预报雾区)。

Exp-A与Exp-B的初始场差异是由同化过程所导致的, 而同化结果具体体现在分析增量上。下面通过分析增量来解释Exp-A与Exp-B的初始场差异。

Exp-B仅比Exp-A多同化了 (表 2), 而只在EnKF循环同化的第2步被同化(图 5)。图 9为第3步同化之前, Exp-B − Exp-A的模式底层结果, 差异主要集中在海上的雾区所在海域。图 10a、图 10b为同化第2步(5月16日21 UTC)的在模式底层的分析增量。分别对比图 9为第3步同化之前, Exp-B−Exp-A的模式底层结果, 差异主要集中在海上的雾区所在海域。图 10a和图 10b为同化第2步(5月16日21 UTC)的在模式底层的分析增量。分别对比图 9a和图 9b与图 10a和图 10b, 发现变量(风、温度、qc与水汽混合比)的空间分布十分相似。在初始场时刻(即EnKF循环同化的第3步)仅同化了obs(探空、地面观测和卫星辐射数据); 第3步同化之前, Exp-B与Exp-A的差异是第2步同化造成的。所以, 第2步中 (图 10c)的同化是造成Exp-B初始场明显优于Exp-A的关键原因。另外, 图 8和图 9的变量分布形势大致相同, 说明在第3步同化过程中, obs未对2个试验的差异产生影响, Exp-B相对于Exp-A的变化主要来源于。

如图 10c, 黄海中部的超过了0.06mm, 由于没有中高云覆盖(根据MTSAT红外数据判断, 图略), 所表征的是海雾中的云水含量的垂直累积。虽然集合成员预报的也在黄海中部显现(图 10d), 但它却呈东西条状分布, 其值为0.04mm左右。与空间分布的差异, 造成了较大的观测增量, 使EnKF同化后产生了较显著的分析增量(图 10a、图 10b)。导致在模式底层, qc明显增加, 高达0.32g/kg(图 10a); 雾区北部出现1—2m/s的西北风增量, 南部出现1m/s的西南风增量。而且, 还导致了温度的降低和水汽的增大(图 10b):温度降低了1—1.4K, 水汽增大了1—1.5g/kg。

图 10的分析增量展示:虽然同化的是, 但模式中的风、温度与湿度皆有反应, 而且它们增量的空间分布十分相近。这应归功于EnKF同化过程中的流依赖背景误差。

EnKF背景误差的“流依赖”主要体现在2个方面:分析增量的空间分布, 控制变量之间的相关性。为了清晰地看到这种“流依赖”特征, 设计了一个单点同化试验:以Exp-B为基础, 在黄海中部某一点(图 11中的黑点P), 将其所同化的替换成一个较大0.01mm的观测值。

点P处的约为0.031mm, 因此被设为0.041mm。在5月16日21 UTC, 同化分析增量如图 11所示。围绕点P, 温度与湿度增量呈椭圆形分布, 前者最大降低约0.3K, 后者增加达到了0.2g/kg; 点P以北区域存在0.1—0.2m/s的西北风增量, 以南区域呈现较弱的偏西风增量。这种椭圆形分布而不是圆形分布, 是流依赖的一个主要特点; 控制变量之间的相关性, 使得一点的的同化就激发了关联变量(温度、湿度与风)的改变。

根据EnKF的增益矩阵(公式(2))和更新算法(公式(5)), 增量的分布形势取决于协方差矩阵PbHT的结构。基于集合成员的背景场, 统计了点P处的与其它变量的协方差分布(图 12)。

在图 12中, 与温度呈负相关, 与湿度呈正相关, 且点P以西和以北区域的协方差较大, 以南的却较小; 与纬向风呈正相关, 与点P附近和以北的经向风呈负相关, 而与以南的经向风呈正相关。若不进行局地化处理, 单点试验中增量的分布将与图 12几乎保持一致(图略)。

值得注意的是, 由于不同集合成员的差异较大, 当统计某一点上的协方差时, 往往有部分成员的为0, 使得实际统计中起作用的样本数减少, 非0的成员影响力增大, 样本误差增大。云水对温湿变量的影响由较强的物理关系来约束, 不易受样本误差影响, 如气温越低, 空气越容易达到饱和。但是, 云水与风的关系较为复杂, 样本误差可能导致云水观测产生不合理的风速增量。譬如, 对本次海雾个例, 偏南风应更有利于雾区向北发展; 但在点P以北, 却与经向风呈负相关(图 12d), 同化虽然增大了雾区但却导致了不利于雾区发展的北风增量(图 11c)。

为了分析风场变化对雾区预报的影响, 进行一个敏感性同化试验:与Exp-B一样, 但控制变量中去除了风速, 即公式(14)变为

该试验记为Exp-S。

与Exp-B相比, Exp-S预报的海雾雾区向北扩展的更符合观测(图略)。对Exp-S预报雾区评分的统计结果表明, Exp-S的POD、FAR、Bias与ETS评分分别为0.650、0.124、0.743与0.497。与表 3给出的Exp-B评分相比, Exp-S的POD、FAR、Bias与ETS评分分别提高了5.3%、24.8%、5.4%与11.7%。这表明, 如果控制变量不包含风, 即不考虑云水与风的关系, 那么同化效果会更好。这同时也暗示, 如果EnKF中无法统计出云水含量与风之间的准确关系, 则要谨慎利用云水含量来改变风场。

在黄海海雾短时临近预报中, 利用数据同化手段改进数值模拟的初始场质量至关重要。Wang等(2014)提出了一种海雾信息的间接同化方法:基于MTSAT静止卫星多通道数据反演出海雾雾体, 认为其内部相对湿度为100%并加以同化。与此方法不同, 本文借助EnKF方法, 提出直接同化海雾雾体云水信息的思路。由于目前尚无法通过静止卫星数据获取海雾雾体内的云水信息, 本文基于极轨卫星反演的云水路径数据进行了EnKF同化试验。针对2011年5月一次大范围的黄海海雾展开同化研究, 结论如下:

(1) EnKF中将云水混合比增加到背景场和分析场的控制变量中, 构建云水观测数据与背景场之间的关系, 实现云水路径数据的直接同化, 是一条可行的途径。

(2) EnKF同化云水路径数据, 可使预报雾区的POD评分提高20.3%, 公正预兆评分ETS提高23.3%。船舶观测的检验结果表明, 海面气温与湿度状态均得到显著改善。

(3) 同化云水路径数据使海雾预报效果大幅改进的原因是, EnKF能动态提供基于集合体所统计的流依赖背景误差协方差。但受集合样本误差的影响, 云水含量与风之间的相关关系的准确性值得特别关注。