2019, Vol. 77
中国气象学会主办。
文章信息
- 邵长亮, 闵锦忠. 2019.
- SHAO Changliang, MIN Jinzhong. 2019.
- 基于地面资料集合均方根滤波同化方案的京津冀暴雨模拟研究
- A numerical study of the rainstorm in Beijing-Tianjin-Hebei region based on assimilation of surface AWS data using the Ensemble Square Root Filter
- 气象学报, 77(2): 233-242.
- Acta Meteorologica Sinica, 77(2): 233-242.
- http://dx.doi.org/10.11676/qxxb2019.008
-
文章历史
- 2017-12-20 收稿
- 2018-07-06 改回
2. 中国气象局气象探测中心, 北京, 100081
2. CMA Meteorological Observation Centre, Beijing 100081, China
地面观测资料是实际观测站地形高度上的观测量，不是模式地形上的观测量，且一般模式地形与实际观测站地形存在一定的高度差异，传统的插值方法并不适用。因此，如何解决模式地形与观测站地形高度差异(以下简称高度差)始终是地面观测资料同化研究工作的关键问题。
Miller等(1992)采用最优插值法(OI法)，通过设计背景误差空间相关模拟函数，合理选择分析变量，进而减小高度差和复杂地形的影响。Ruggiero等(1996)考虑到模式地形与实际观测站地形高度有一定的差异，利用近地层相似理论同化地面观测资料时，分3种情况进行处理：(1)测站地形高度高于模式最低层高度，则将地面观测资料作为高空资料进入模式；(2)模式最低层高度高出测站地形的高度超过100 m，则该站点资料剔除不用；(3)模式最低层高度比测站地形高度高，且模式最低层高度高出测站地形的高度小于100 m时，该站点资料则利用背景场信息将观测资料反演到模式最低层。这种方案最大的优点就是不需要观测算子和相应的切线性模式和伴随模式，但是它也将许多观测资料剔除了，并且高于模式最低层的测站资料被当作探空资料处理会导致变量的不连续(即与真实上层大气相比，风场会变弱，温度会变高或变低)，进入模式会造成不协调；Urban(1996)认为当高度差异很大时，不能通过简单的内插外延模式值而获得对近地面资料的最佳估计。并设计了对温度和雪深的观测算子，其思想方法是参数化一个物理过程的概率分布，然后通过求数学期望得到相应的模式值。这实际上是在统计意义上考虑了高度差异。为了解决Ruggiero等(1996)方案中的不足，Guo等(2002)设计的新方案则不考虑实际观测站地形与模式地形高度的差异，而是假定所有测站的资料(除地面气压)都是位于模式面，然后利用相似理论建立10 m高风场(U10, V10)和2 m高温度(T2)、湿度(q2)的观测算子及相应的切线性模式和伴随模式，同时在进行极小化运算前将地面气压(psfc)折算到模式最低层。实现该方案采用的同化分析模式是MM5-3DVAR，地面观测资料得到了充分利用，并取得了一定的效果。这种方案比前一种方案更充分地利用观测资料，但却没有考虑模式地形与实际测站地形的高度差异，对于地形分布复杂，模式地形与观测站地形高度差异较大的区域，由于温度、气压、风场等观测资料是随着海拔高度和地形分布而变化的，差异较大就会造成模式中各物理量梯度不协调。Lazarus等(2002)在ADAS系统(the ARPS (Advanced Regional Prediction System) Data Analysis System)客观分析复杂地形下地面观测和探空观测资料时，也考虑了高度差异，并在进行分析时所取的客观分析权重系数中不仅包含了水平距离权重系数，还增加了模式分析层高度与观测站地形高度差权重系数以及模式地形与实际观测地形高度差权重系数，有效地减轻了高度差异在资料客观分析过程可能带来的负面效果。Devenyi等(2003)和Benjamin等(2004)在NCEP的RUC(Rapid Update Cycle)系统中分别采用最优插值(OI)法和3维变分同化法，通过采用局地递减率将地面观测气压、温度和湿度由观测站地形订正到模式地形高度，解决高度差异问题。他们对气压和风的观测资料不考虑测站和模式高度差，仅对气压在模式最下面5层使用局地递减率订正。而温度和露点温度资料假定高度差大于70 hPa则被剔除，否则采用局地递减率由观测站地形订正到模式地形高度，然后利用近地层相似性理论将2 m高温度和湿度观测资料与10 m高风匹配成背景场。但引进地面观测资料同化后，改进效果不是很明显。Deng等(2005)采用Bratseth逐步订正法，先用Mother-Daugther方法对复杂地形下的地面温度资料进行客观分析，将山谷和山顶的观测资料分开处理，分别求得山谷和山顶的环流传播距离和分享因子，再将这两个参数用于各项异性的背景误差相关函数中，从而使山谷资料沿复杂地形下的山谷环流传播山谷地面的信息，减少对山脊上空的影响；山顶资料信息向邻近山顶传播。由于只同化了温度资料，产生了初始场不平衡，对风场的预报效果较差；Lee等(2005)设计了以莫宁-奥布霍夫相似理论为基础的近地面综合观测算子；Myrick等(2005)对复杂地形下的背景误差相关函数考虑格点和观测点地理位置作了局部调整。Fujita等(2006)在MM5模式中应用集合卡尔曼滤波方法，同化风场(U, V)、位温(θ)和露点(Td)，而不同化气压(p)，从而降低高度差异带来的影响。徐枝芳等(2007)认为，中国的地形比较复杂，模式地形与实际观测站地形在许多地区存在较大差异，地面观测资料同化方案设计中有必要考虑模式与实际观测站地形高度差异，而Guo等(2002)的地面观测资料同化方案没有考虑模式与实际观测站地形高度差异，因此不能将该方案直接应用到中国，于是其基于MM5_3DVAR系统对Guo等(2002)采用的方法又提出了改进：在地面观测误差中增加模式地形与观测站地形高度差异引起的地形代表性误差。这个同化方案有效合理地将地面资料同化到了数值模式中，改进了暴雨模拟结果。但当模式地形与观测站地形高度差异较大时，两种改进方案都存在一定敏感性。而且，增加地形代表性误差会导致地面观测误差不满足无偏假定，分析场也不能达到最优。Stensrud等(2009)在WRF-DART(Data Assimilation Research Testbed)系统中使用集合卡尔曼滤波方法直接同化10 m高风场(U10, V10)和2 m高位温(θ2)、露点(Td2)，但是未考虑高度差异。Pu等(2013)基于WRF模式分别应用集合卡尔曼滤波方法与3维变分同化方法在复杂地形下直接同化地面观测，结果表明在复杂地形下，集合卡尔曼滤波通过背景误差协方差的流依赖特性可以产生更加接近实际的分析增量。邵长亮等(2015)在WRF模式中应用集合均方根滤波直接同化10 m高风场(U10, V10)、2 m高位温(θ2)、2 m高露点(Td2)和地表气压(psfc)，能够改善18 h累计降水预报，但除地表气压外，其他要素的同化没有考虑高度差异。Chen等(2016)在同化雷达观测资料的基础上同化地面观测资料，进一步改善了分析和预报结果，其中地面观测资料可以改善低层的分析和预报，使得地面与低层的飑线特征更加接近实况观测。
文中根据邵长亮等(2015)研究工作，考虑模式地形与观测站地形高度差，设计对应的同化方案，并进行个例试验。减少模式地形与观测站地形高度差对地面资料同化的影响。
2 改进的自动气象站观测资料同化方案 2.1 地形代表性误差的引入文中考虑模式地形高度与实际地形高度的差异，在位温和露点观测误差中引入地形代表性误差，并由温度地形代表性误差求出位温地形代表性误差和露点地形代表性误差，进而研究其对地面自动气象站观测资料同化效果的影响。
温度和位温、露点的观测误差分别为Te=Ti+Tr+Tt和θe=θi+θr+θt、Tde=Tdi+Tdr+Tdt, 其中, Te和θe、Tde分别为温度和位温、露点的总观测误差, Ti和θi、Tdi为仪器误差, Tr和θr、Tdr为代表性误差, Tt和θt、Tdt为地形代表性误差，以上变量单位均为K。试验采用的地面观测误差在所有观测站点都不相关，但大小一致且均为正值。位温地形代表性误差(θt)和露点地形代表性误差(Tdt)由以下公式计算得到
|
(1) |
|
(2) |
|
(3) |
|
(4) |
式中，T、p、f为观测值，其中T为2 m高温度，单位为K；p为地表气压，单位为hPa，f为2 m高相对湿度；R=8.31 J/(mol·K)为理想气体常数，cp=29.1 J/(kg·K)为理想气体比定压热容，Rv=461.6 J/(kg·K)为水汽比气体常数，Lv=2501 J/kg为水的汽化潜热，
|
将式(7)由泰勒级数展开并取一阶近似得
|
(9) |
|
(10) |
由以上公式得
|
(11) |
地形代表性误差主要是在观测误差中加入合适的位温和露点地形代表性误差。温度地形代表性误差(Tt)由|Δh|·|γ|求得(Δh为模式与观测站地形高度差，单位为m，γ为温度垂直递减率，单位为K/m)，即地形高度差的绝对值、温度递减率的绝对值与温度地形代表性误差成正比。从物理意义上讲，地形高度差的绝对值、温度递减率的绝对值越大，则引入到观测误差中的值越大，从而观测误差值越大，该观测值的可信度越低，反之亦然。位温地形代表性误差(θt)和露点地形代表性误差(Tdt)是用温度地形代表性误差(Tt)代入式(1)得到，所以温度垂直递减率(γ)是引入地形代表性误差的重要参数。
2.3 温度垂直递减率(γ)的设计实际上，模式与观测站地形高度之间没有真实大气，因而不适合使用模式层计算γ。这里使用邻近探空站观测数据计算温度垂直递减率，即某一站点的γ由邻近探空站的温度垂直递减率决定。
3 试验设计 3.1 暴雨个例试验的个例为2016年7月19日00时至21日00时(世界时，下同)出现在京津冀地区的一次暴雨过程(图 3)。
|
| 图 3 2016年7月19日00时—20日00时(a)与20日00时—21日00时(b)24 h累计降水 (a1、b1.自动气象站观测；a2、b2.试验1，a3、b3.试验2，a4、b4.试验3；单位：mm) Fig. 3 24 h precipitation (mm) during 00:00 UTC 19 Jul-00:00 UTC 20 Jul 2016 (a) and during 00:00 UTC 20 Jul-00:00 UTC 21 Jul 2016 (b) (a1, b1: AWS; a2, b2: Expt 1; a3, b3: Expt 2; a4, b4: Expt 3) |
模式参数设置：采用WRF模式单重区域，水平格距10 km，试验区中心点为(37°N，117.5°E)，格点数取150×150，垂直方向共33层，模式层顶为100 hPa。物理方案采用Lin微物理方案和Kain-Fritsch积云对流参数化方案。对应的物理过程选用MYJ边界层方案，近地面层采用莫宁-奥布霍夫方案，长波辐射方案为RRTM(Rapid Radiative Transfer Model)，短波辐射方案采用Dudhia，陆面过程采用Noah方案，模式积分时间为2016年7月18日12时至7月21日00时，积分步长30 s。使用集合均方根滤波(EnSRF)同化方法。集合数40，初始扰动为全场为0的随机扰动，三维风速U、V、W的扰动标准差分别为2、2、0.5 m/s，T的扰动标准差为2 K，水汽混合比Qv的扰动标准差为0.002 kg/kg(此处U、V、W、T、Qv为背景场要素值)。使用常数膨胀算法进行协方差膨胀(Tong, et al, 2005)，膨胀系数为1.25；局地化使用五阶距离相关函数方案(Gaspari, et al, 1999)使用Schur算子。水平和垂直局地化距离分别为45和10 km；所有试验的模式起始时间为18日12时，16时加入扰动，17时开始同化第1次自动气象站观测资料，共同化8个时次至19日00时，然后开始预报至21日00时结束。
基本试验资料为2016年7月18日12时—7月21日00时，逐6 h的NCEP 1°×1° FNL背景场资料和18日17时—19日00时逐时的京津冀地面自动气象站观测资料，包括175个中国国家级自动气象站和2559个区域自动气象站，其中，中国国家级自动气象站无缺测，区域自动气象站包含很多缺测数据。经变量变换(邵长亮等，2015)后，每个时次可用要素个数分别为经向风746、纬向风746，位温429，露点443，地表气压429。18日12时和19日00时两次京津冀探空观测资料，共4个探空站数据。图 1为2016年7月18日17时175个中国国家级自动气象站观测和4个探空站分布及中国国家级自动气象站实际地形高度与模式地形高度差，首先将模式地形高度插值到中国国家级自动气象站位置得到模式地形高度值，再用中国国家级自动气象站高度值减去得到的模式地形高度。该图反映了模式地形高度与中国国家级自动气象站地形高度的差异。
|
| 图 1 2016年7月18日17时自动气象站和探空站观测分布和自动气象站实际地形高度与模式地形高度差(色阶，单位：m) (绿色十字代表自动气象站，蓝色圆圈代表探空站) Fig. 1 Distribution of AWS observations at 17:00 UTC 18 Jul 2016 and elevation differences (shded, unit: m) between surface AWS and numerical model topography, AWS location (green cross) and sounding location (blue circle) |
试验方案见表 1，试验1为对照试验，不同化地面资料；试验2同化地面资料但γ=0, 即不引入地形代表性误差；试验3引入地形代表性误差。各观测站点的γ由邻近探空站的温度垂直递减率决定，由于同化的自动气象站观测时次为18日17时—19日00时，因而取18日12时和19日00时两个时次的平均代表该时段内温度垂直递减率，从而得到53798、54401、54511和54539四个站点的温度垂直递减率分别为0.466、0.662、0.571和0.455 K/(100 m)，进而分析各试验结果。
| 试验 | 自动气象站观测资料 | 温度地形代表性误差(单位：K) | EnSRF同化处理 |
| 1 | 否 | — | 否 |
| 2 | 是 (U10, V10, θ2, Td2, psfc) |
0 | 是 |
| 3 | 是 (U10, V10, θ2, Td2, psfc) |
|Δh|·|γ| | 是 |
对应京津冀地区，将各时次背景场要素值插值到观测站位置，得到各试验RMSE_U10、RMSE_V10、RMSE_θ2、RMSE_Td2和RMSE_psfc随时间(时次)的变化(图 2)。试验2和试验3各要素均方根误差随时间(时次)均呈下降趋势，且在同化结束的时次均明显低于试验1；试验3与试验2相比，风场均方根误差在各时次均较小，RMSE_θ2在前5个时次较高，后3个时次较低，RMSE_Td2在第6和第8个时次较低，其余时次则较高。同时，从同化第1个时次的各试验的均方根误差可以发现，试验2和试验3在同化资料前，集合预报阶段结束时RMSE_U10、RMSE_θ2、RMSE_Td2和RMSE_psfc均高于对照试验，RMSE_V10低于对照试验。结果表明，在同化资料前，各要素的均方根误差除V风场外，均高于对照试验；在同化阶段，同化自动气象站观测资料能有效降低各要素的均方根误差；在观测误差中引入地形代表性误差与引入前相比，风场均方根误差得到整体改进；位温和露点的均方根误差在前期表现并不稳定，后期有所改进；地表气压均方根误差差别不大。
|
| 图 2 各变量均方根误差随时间(时次)变化 (a. RMSE_U10，b. RMSE_V10，c. RMSE_θ2，d. RMSE_Td2，e. RMSE_psfc) Fig. 2 RMSEs of various variables (a. RMSE_U10, b. RMSE_V10, c. RMSE_θ2, d. RMSE_Td2, e. RMSE_psfc) |
实况资料为自动气象站累计降水观测，降水分为两个阶段，第1个阶段为2016年7月19日00时—20日00时，降水中心在河北西南部地区；第2个阶段为2016年7月20日00时—21日00时，降水中心有两个，分别在北京南部和河北东部地区，图 3为实况降水和模拟累计降水。各试验模拟出的雨带和走势基本一致，均模拟出了降水中心；试验2和试验3与试验1相比，在第1个阶段降水量在京津冀大部分地区有整体增强，与实况更加接近；在第2个阶段北京地区降水中心范围扩大，消除了河北西部地区的虚假降水中心，并且河北东部地区降水增强，更加接近实况。对应京津冀地区，将降水模拟结果插值到观测站位置进行降水TS(Threat Score)评分(图 4)，可见第1阶段试验3累计降水评分均高于试验2和试验1；第2阶段24 h累计降水除小雨和暴雨量级外，试验3评分均高于试验2和试验1。
|
| 图 4 各试验降水TS评分 (a. 2016年7月19日00时—20日00时，b. 20日00时—21日00时) Fig. 4 TS scores of precipitation (a. 00:00 UTC 19 Jul-00:00 UTC 20 Jul 2016, b. 00:00 UTC 20 Jul-00:00 UTC 21 Jul 2016) |
将试验3与试验2的两个阶段累计降水相减，可以得到引入地形代表性误差后与引入前累计降水分布差异(图 5)，可见随着两个阶段降水中心位置由京津冀西南地区向东北地区移动，差值较大区域也随之移动，大部分集中在降水中心附近；同时，正、负差值基本上成对分布且量级相近，说明引入地形代表性误差对累计降水的改变主要表现为对降水位置的调整。
|
| 图 5 试验3与试验2的24 h累计降水模拟结果差值 (a. 2016年7月19日00时—20日00时，b. 20日00时—21日00时) Fig. 5 Differences in simulated 24 h accumulated rainfall between Expt 3 and Expt 2 (Expt 3-Expt 2) (a. 00:00 UTC 19 Jul-00:00 UTC 20 Jul 2016, b. 00:00 UTC 20 Jul-00:00 UTC 21 Jul 2016) |
预报阶段，对风场(U10, V10)、温度(T2)、相对湿度(f)和地表气压(psfc)进行评估，由图 6可见，试验2、试验3与试验1相比，RMSE_U10在前28个时次较小，之后较大；RMSE_V10在前18个时次较小，之后较大；RMSE_T2整体上较小；RMSE_f在前13个时次较小，之后基本一致；RMSE_psfc在前32 h基本一致，之后较大。试验2与试验3各要素的均方根误差较为一致。总体上，在预报阶段，同化自动气象站资料与对照试验相比，前13个时次各要素预报均有所改进，更加接近实际观测，而之后各要素预报并不稳定；引入地形代表性误差与不引入地形代表性误差相比，对要素预报并无明显改进。
|
| 图 6 各试验2016年7月19日00时—21日00时各要素均方根误差 (a. RMSE_U10，b. RMSE_V10，c. RMSE_T2，d. RMSE_f，e. RMSE_psfc) Fig. 6 RMSE of various fields simulated during 0000 UTC 19 Jul-0000 UTC 21 Jul 2016 (a. RMSE_U10, b. RMSE_V10, c. RMSE_T2, d. RMSE_f, e. RMSE_psfc) |
通过在位温和露点观测误差中分别引入位温和露点地形代表性误差，来解决模式地形与观测站地形存在的高度差对地面资料同化效果的影响，并对一次暴雨个例进行了数值试验。研究结果表明：同化与不同化地面观测资料相比，同化阶段的均方根误差、预报阶段的降水TS评分和前13个时次各要素预报均有整体改进；在观测误差中引入地形代表性误差与引入前相比，风场均方根误差得到整体改进；位温和露点的均方根误差在前期表现并不稳定，后期有所改进；在观测误差中引入地形代表性误差能够影响降水模拟结果，预报阶段前24 h累计降水与后24 h累计降水TS评分在整体上均有所提高；引入地形代表性误差对要素预报并无明显改进。综上所述，新方案能够减少模式地形与观测站地形高度差异对地面资料同化的影响。
文中对(U10, V10)进行了直接同化，但并未考虑地形差异，研究过程中曾使用风向、风速代替(U10, V10)进行直接同化(Pu, et al, 2013)，但该方案在文中个例中试验效果并不好，因此，文中未使用。
| 邵长亮, 闵锦忠. 2015. 集合均方根滤波同化地面自动站资料的技术研究. 大气科学, 39(1): 1–11. Shao C L, Min J Z. 2015. A study of the assimilation of surface automatic weather station data using the ensemble square root filter. Chinese J Atmos Sci, 39(1): 1–11. (in Chinese) |
| 徐枝芳, 龚建东, 王建捷, 等. 2007. 复杂地形下地面观测资料同化Ⅱ:模式地形与观测站地形高度差异代表性误差. 大气科学, 31(3): 449–458. Xu Z F, Gong J D, Wang J J, et al. 2007. A study of assimilation of surface observational data in complex terrain part Ⅱ:Representative error of the elevation difference between model surface and observation site. Chinese J Atmos Sci, 31(3): 449–458. (in Chinese) |
| Benjamin S G, Dévényi D, Weygandt S S, et al. 2004. An hourly assimilation-forecast cycle:The RUC. Mon Wea Rev, 132(2): 495–518. DOI:10.1175/1520-0493(2004)132<0495:AHACTR>2.0.CO;2 |
| Chen X C, Zhao K, Sun J Z, et al. 2016. Assimilating surface observations in a four-dimensional variational Doppler radar data assimilation system to improve the analysis and forecast of a squall line case. Adv Atmos Sci, 33(10): 1106–1119. DOI:10.1007/s00376-016-5290-0 |
| Deng X X, Stull R. 2005. A mesoscale analysis method for surface potential temperature in mountainous and coastal terrain. Mon Wea Rev, 133(2): 389–408. DOI:10.1175/MWR-2859.1 |
| Devenyi D, Benjamin S G. 2003. A variational assimilation technique in a hybrid isentropic-sigma coordinate. Meteor Atmos Phys, 82(1-4): 245–257. DOI:10.1007/s00703-001-0590-y |
| Fujita T, Stensrud D J, Dowell D C. 2006. Surface data assimilation using an ensemble Kalman filter approach with initial condition and model physics uncertainties. Mon Wea Rev, 135(5): 1846–1868. |
| Gaspari G, Cohn S E. 1999. Construction of correlation functions in two and three dimensions. Quart J Roy Meteoro Soc, 125(554): 723–757. DOI:10.1002/(ISSN)1477-870X |
| Guo Y R, Shin D H, Lee J H, et al. 2002. Application of the MM5 3DVAR system for a heavy rain case over the Korean peninsula//Papers Presented at the Twelfth PSU/NCAR Mesoscale Model Users' Workshop NCAR. Boulder |
| Lazarus S M, Ciliberti C M, Horel J D, et al. 2002. Near-real-time applications of a mesoscale analysis system to complex Terrain. Wea Forecasting, 17(5): 971–1000. DOI:10.1175/1520-0434(2002)017<0971:NRTAOA>2.0.CO;2 |
| Lee S J, Parrish D F, Wu W S. 2005. Near-surface data assimilation in the NCEP Gridpoint Statistical-Interpolation System: Use of land temperature data and a comprehensive forward model. Office Note, 446 |
| Miller P A, Benjamin S G. 1992. A system for the hourly assimilation of surface observations in mountainous and flat terrain. Mon Wea Rev, 120(10): 2342–2359. DOI:10.1175/1520-0493(1992)120<2342:ASFTHA>2.0.CO;2 |
| Myrick D T, Horel J D, Lazarus S M. 2005. Local adjustment of the background error correlation for surface analyses over complex terrain. Wea Forecasting, 20(2): 149–160. DOI:10.1175/WAF847.1 |
| Pu Z X, Zhang H L, Anderson J. 2013. Ensemble Kalman filter assimilation of near-surface observations over complex terrain:Comparison with 3DVAR for short-range forecasts. Tellus A:Dyn Meteor Oceanogr, 65(1): 19620. DOI:10.3402/tellusa.v65i0.19620 |
| Ruggiero F H, Sashegyi K D, Madala R V, et al. 1996. The use of surface observations in four-dimensional data assimilation using a mesoscale model. Mon Wea Rev, 124(5): 1018–1033. DOI:10.1175/1520-0493(1996)124<1018:TUOSOI>2.0.CO;2 |
| Stensrud D J, Yussouf N, Dowell D C, et al. 2009. Assimilating surface data into a mesoscale model ensemble:Cold pool analyses from spring 2007. Atmos Res, 93(1-3): 207–220. DOI:10.1016/j.atmosres.2008.10.009 |
| Tong M J, Xue M. 2005. Ensemble Kalman filter assimilation of Doppler radar data with a compressible nonhydrostatic model:OSS experiments. Mon Wea Rev, 133(7): 1789–1807. DOI:10.1175/MWR2898.1 |
| Urban B. 1996. Coherent observation operators for surface data assimilation with application to snow depth. J Appl Meteor, 35(2): 258–270. DOI:10.1175/1520-0450(1996)035<0258:COOFSD>2.0.CO;2 |