基于数据同化的降雨数值空间分布模拟研究

王 洋，刘 佳，于福亮，李传哲，田济扬，邱庆泰

（中国水利水电科学研究院 流域水循环模拟与调控国家重点实验室，北京 100038）

1 研究背景

水文预报是解决和预防暴雨、洪水等灾害的有效科学手段，通过准确、及时的预报，可以提高政府部门应对洪水灾害的主动性，防患于未然。水文预报首先需要解决的问题是延长预见期和提高预报的精度，就某个流域来说，特别是中小尺度流域，由于流域面积小，汇流时间短，如果以落地雨作为水文预报的降雨数据来源，则预见期无法达到水文预报的要求[1]。目前常用的解决方法是通过数值大气模式和流域水文模型进行耦合预报。数值大气模式可以提供具有一定预见期的降雨数据，将该数据输入流域水文模型可以延长预报时段[2]。Lu等[3]在淮河流域进行了陆气耦合实验研究，在保证模拟效果的前提下发现陆气耦合的方法可获得比传统洪水预报更长的预见期，并将该成果应用于流域洪水预警。于鑫等[4]将WRF（Weather Research and Forecasting Model）模式和HEC-HMS水文模型进行单向耦合，对研究区进行模拟研究，发现可通过滚动预报延长洪水预见期。

大气数值模式在进行降雨模拟时，除了可以提供未来的降雨数据，还可以通过多层嵌套的降尺度方法，提供从几百公里到一公里范围的高分辨率降雨信息，该降雨信息对降雨空间分布刻画的准确与否直接关系到水文预报的预报精度。特别是在采用分布式水文模型建立陆气耦合系统时，若大气数值模式输出降雨数据和实际的降雨空间分布不一致，则会给水文模型提供错误的降雨空间分布信息，从而影响水文模型模拟的精度。Hamill等[5]运用WRF模式对美国某地区的一次降雨过程进行了模拟，结果表明WRF模式的降雨输入结果可以较好地表征降雨的空间分布情况。Kryza等[6]在波兰对大气数值模式开展研究，结果表明当降雨强度较小时，WRF模式的模拟效果好于降雨强度大时的结果。Jasper等[7]运用WaSIM-ETH水文模型和大气模式进行耦合研究，结果显示提高数值大气模式的降雨模拟精度可以降低陆气耦合系统的模拟误差。

降雨数值模拟中驱动数据所提供的初始场和边界场条件和大气的实际状态并不是完全吻合，导致模拟结果存在误差。数据同化方法是提高数值大气模式输出降雨数据精度的有效方法，通过将数值模型拟合结果与观测数据相融合，不断更新系统状态与参数，最终实现改进降雨预报精度的效果。数据同化技术早在1950年代就被应用于数值天气预报，常用的数据同化方法主要有顺序同化和变分同化两大类，其中变分同化中的三维变分同化是使用较为广泛的同化方法。Routray等[8]是最早使用三维变分同化方法改善降雨预报数据的研究者。范水勇等[9]通过对北京地区数据同化前后降雨预报值的比较，验证了三维变分数据同化方法可以提高降雨的模拟精度。黄兵等[10]选用三维变分同化方法同化观测资料，数值模拟结果显示对局地暴雨强度和落区都有较好的反映。

综上所述，在洪水突发性强的中小尺度流域，数值大气模拟的降雨空间分布情况直接影响水文预报的结果。由于大气数值模式驱动数据的初始场和侧边界条件与大气的实际状态间存在差异，导致模拟的降雨数据空间分布和实际降雨分布之间存在误差，数据同化是解决该问题的有效手段。本文采用WRF模式开展基于数据同化的降雨数值空间分布模拟研究，探求数据同化对降雨模拟数据空间分布的改进效果，为进一步提高陆气耦合水文预报能力提供依据。

2 研究区域概况

本文选取位于大清河南支的阜平流域作为研究区。该流域位于东经113°45′至114°32′，北纬38°39′至39°8′之间，流域总面积为2219 km2，流域控制站为阜平水文站，流域位于大清河南支的支沙河上游。阜平流域全流域为山区，高程自西北往东南递减，高程变化大，高程最大值约为2290 m，最小值约为200 m，区域内坡度较陡。该流域属于暖温带半湿润地区，气候为北温带大陆性季风气候，雨热同期，且夏季多暴雨。

选取该流域作为研究区主要有以下几点原因：

（1）研究区内的降雨有突发性强、历时短、雨量大的特征，特别是在山区发生该类型降雨时，极易引起洪涝灾害。由于流域面积小，预见期短，预报难度较大，通过在该流域开展降雨数值空间分布模拟研究，可寻找提高流域水文预报能力的新思路。

图1 阜平流域地理位置示意

（2）研究区在大清河水系的南支上，对于分析整个大清河水系的降雨产流机理具有重要意义。雄安新区位于大清河流域的腹地，通过该区域开展降雨数值空间分布模拟研究，对调整完善流域防洪预案、保障雄安新区及下游防洪安全、配合雄安新区规划建设具有重要的现实与长远意义。

3 研究数据

3.1 全球模式产品WRF模式运行需要初始边界场作为驱动，本文采用NCEP的全球预报数据GFS作为WRF模式运行的初始场。GFS数据分辨率为1°×1°，其预报数据可预报未来8 d共192 h的天气，数据每隔6 h更新一次。目前全球模式产品在东亚地区的模拟结果相较于欧洲和美国还存在一定的差异，但模式中温度、风速、压力等要素的描述在东亚地区还是可以刻画出一定的趋势和规律性，故实际上GFS数据基本可以满足模拟我国中小尺度区域气象要素变化的要求[11]，GFS数据来源于（http：//nomads.ncep.noaa.gov/）。

3.2 数据同化资料本文选择多普勒天气雷达的反射率资料和NCAR提供的GTS全球观测资料作为数据同化资料。GTS数据是气压、温度、湿度、风向与风速等传统的地面与高空气象数据的集合，包括飞机报、船舶报、卫星资料和地面观测站数据等[12]。研究所采用的多普勒雷达是位于石家庄新乐市境内的SA波段雷达，雷达数据由国家气象局提供。该多普勒气象雷达扫描区域半径为250 km，可完全覆盖研究区。

3.3 实测降雨资料实测降雨资料主要来自于阜平流域的8个雨量站，分别为：冉庄、庄旺、不老台、下关、砂窝、龙泉关、桥南沟和阜平。通过对历史资料的分析，选取3场降雨过程作为研究对象，各场降雨雨量大小和历时见表1。

表1 降雨场次信息

4 研究方法

4.1 WRF模式WRF模式全称为Weather Research Forecast Model，译为天气预报模式，被誉为新一代的中尺度天气预报模式。该模式是一种完全可压缩及非静力模式，模式系统可以进行数值天气的预报，还可进行天气模拟和数据同化。WRF模式有两个版本，一个是由MM5模式发展而来，采用Arakawa C网格，另一个是由ETA模式发展而来，采用Arakawa E网格。相比于其他大气数值模式，WRF模式支持高精度的网格分辨率，最高可生成1km×1km的网格数据，采用了比ETA更好的地形模拟数据，保证了在长步长条件下运算的稳定性。WRF模式分为ARW和NMM两种，即研究用和业务用两种形式，本文中使用的是WRF ARW。WRF的运行需要进行物理参数化方案的设置，针对各场降雨过程的具体特点，研究选取的影响降雨过程的主要物理参数化方案如表2所示。

表2 各降雨场次选取的主要物理参数化方案情况

4.2 三维变分数据同化WRF模式中的WRFDA系统是WRF模式自身的数据同化系统，研究选取其中的3DVAR即三维变分同化方法。三维变分同化的基本理论是将同化问题归结为一个表征分析场与观测场和分析场与背景场偏差的二次最小泛函极小值问题，该问题可以通过以下公式表述：

式中：X为分析变量；Xb为背景场；Y0为观测向量；B为背景场误差协方差；H为观测算子，将模式变量由模式空间投影到观测空间；R为观测误差协方差。

则三维变分同化问题最终归结为求解上式给出的目标函数最优解问题。由于在求解最优值的过程中控制变量是一组非线性方程，其求解形式多样，无法直接求出控制变量的闭合解，在WAF 3DVAR中采用外循环和内循环共同迭代的方法得到最优解对应的控制变量。

WRF模式运行需进行嵌套网格的设置，网格的设置一方面需要考虑研究区的位置，另一方面需要匹配同化数据的精度。同化数据的精度和嵌套网格的精度越接近，同化效果往往越好。由于GTS数据的网格精度较粗，适合在外层网格进行同化，雷达数据的空间分辨率较高，适合在内层网格进行同化。参考文献[13]在大清河流域的研究成果，同时考虑计算效率问题，研究中采用双层嵌套网格设置，在内层网格同化雷达反射率数据，外层网格同化GTS数据，将最终输出的内层网格数据（Domain 2）作为模式的模拟结果。WRF模式运行的基本设置见表3。

表3 WRF模式运行参数设置

由于WRF模式在计算过程中误差存在累积效应，随着计算时长的增加，模式的误差会逐渐增大。为避免此类问题，数据同化计算中参考常用的循环同化方法，将前一个时段WRF模式运行得到的背景场和侧边界条件作为下一个时段的初始条件进行同化，从而达到减小模式计算误差的目的。假设某一场降雨开始于2012年7月21日03∶00，结束于2012年7月22日03∶00，降雨历时为24 h。数据同化的起始时间为2012年7月21日00∶00，同化时间间隔为6h，共运行6次。则该场次降雨数据同化后的输出数据来源如图2所示。

图2中，run1代表WRF模式第1次模拟过程，run2-run6代表WRF模式第1次至第5次数据同化过程，run1并没有提供同化后的输出数据，该过程主要作用是对模型进行预热，预热后模型在同化过程中的表现会更加稳定。

4.3 降雨空间分布评价指标采用临界成功率指标CSI和均方根误差指标RMSE对同化前后的降雨数据进行评价，均方根误差指标RMSE计算公式为：

图2 循环同化方案示意

式中：Pj为某一个网格点上的降雨量模拟值；Qj为某一个网格点上整个观测时段内累积雨量的实测值；M为网格个数。

临界成功率指标CSI表示正确模拟的降雨次数占所有可能发生降雨情况的比例，该指标最优值为1，CSI值越接近于1，表明模拟效果越好。该指标计算公式为：

根据对历史降雨资料的分析，当某一网格上的降雨量大于0.01 mm时，说明该网格在统计时段内发生了降雨。当某一个网格上的模拟值与实测值都显示了降雨的发生时，NA值加1；当某一网格上模拟值显示了降雨的发生，而实际没有发生降雨时，NB值加1；反之当网格上模拟值显示没有发生降雨，而实际发生了降雨，则NC值加1。

5 结果分析

5.1 降雨过程代表性分析从阜平流域历史降雨资料中选取3场降雨过程作为研究对象，各场次降雨过程柱状图见图3。

图3 典型降雨过程柱状图

图4 雨量站实测累计降雨空间分布

将雨量站的实测累计降雨数据通过泰森多边形法展布到流域面上，结果如图4所示。

从图4可以看出，3个降雨过程在空间分布的均匀程度并不相同。为进一步量化降雨的空间分布情况，选取变差系数为指标进行描述。空间分布的变差系数反映的是某一雨量站的降雨量距平数与各雨量站平均降水量之间的比值，变差系数Cv按下式计算：

式中：xi为某一雨量站的累积降雨量；为各雨量站的累积降雨量的平均值；n为网格个数；在计算降雨时间分布时，xi为第i个小时的面雨量，为i个小时面雨量的平均值，n为降雨小时数。

各降雨场次变差系数Cv计算结果如表3所示：

由表3计算可见，降雨场次2的Cv值最小，说明某一雨量站的降雨量的距平数和该雨量站平均降雨量之间的差值小，反映出该场降雨空间分布均匀，而降雨场次3的Cv值最大，表明该场降雨过程的降雨量空间分布差距较大，属于空间分布不均匀降雨，降雨场次2的Cv值介于场次1和场次3之间。综合降雨过程的空间Cv值和雨量分布图，可以看出3场降雨的空间分布均有一定的差异性，可以用来代表不同的降雨空间分布类型，代表性较好。

5.2 降雨空间分布的评价指标分析计算的各场降雨临界成功率指标CSI和均方根误差指标RMSE见表4。

表3 3场降雨的Cv值

表4 同化前后评价指标对比

CSI表示正确模拟的降雨频次占所有可能发生降雨情况的比例，RMSE则是可以更加直观的反映空间尺度下模拟误差的大小。从表4中可以看出，对于各场降雨，同化后的CSI值和RMSE值均优于同化前。降雨场次1同化前后CSI值提高了0.1353，降雨场次2同化前后CSI值提高了0.0659，降雨场次3同化前后CSI值提高了0.1031，CSI值的提高说明数据同化方法有助于模式捕捉不同降雨场次的降雨落区；同化前后3场降雨的RMSE指标分别下降了0.1644、0.1120和0.2407，说明通过数据同化的方法，对空间上各点时段累计降雨量模拟的正确率均有提高。综合来看，同化后的数据可以更加准确地刻画降雨的空间分布规律，模拟降雨空间精度的提高有助于在陆气耦合中为陆面水文模型提供更加准确的降雨信息，从而得到准确的预报结果。

从图2中可以看出，降雨场次3属于短时强降雨过程，根据田济扬等[14-15]的研究表明，WRF模式对该类型降雨的模拟效果较差。从表4中也可以看出，虽然同化后CSI值和RMSE值均有优化，但总体程度差于降雨场次1和降雨场次2的模拟结果，在实际水文预报过程中需提供对该类型降雨的关注。

5.3 降雨模拟结果的空间分布通过网格中心点的坐标读取WRF模式同化前后内层嵌套网格的模拟降雨数据，并展布在3 km×3 km的网格图上，网格中心点位置和流域之间的相对位置如图4所示。

图5 网格中心点分布

图6 降雨场次1同化前后降雨空间分布

图5—图7分别给出了3个降雨过程在数据同化前和同化后的降雨空间分布情况，图中柱状图为雨量站实测雨量值。在降雨场次1中，根据实测资料显示，降雨主要集中在中部和南部部分区域，流域东、西、北部边界处的累计降雨量小于流域中心的累计降雨量，流域出口处的累计降雨量较小。未同化前的累计雨量分布表现为中部高，东西低，在流域的出口处，同化前模拟结果高估了累计降雨量，东部区域同化前的模拟结果则低估了累计降雨量。通过数据同化方法，模拟的降雨结果捕捉到东部地区累计降雨量较大的区域，对流域出口处累计降雨量偏小区域的模拟效果也更好，总体上更接近实际的降雨空间分布。

图7 降雨场次2同化前后降雨空间分布

图8 降雨场次3同化前后降雨空间分布

降雨场次2降雨主要集中在南部地区，但南北总体的累计降雨量差距不大，在流域面上的分布较为平均，降雨量分布整体呈现从东南往西北减少的趋势。同化前模拟降雨结果的空间分布也较平均，但同化后的降雨分布图相比于同化前更好地反映了该降雨分布从东南往西北减少的趋势，特别是对北部降雨较少地区的模拟更加准确。

降雨场次3同化前后CSI值和RMSE值虽然偏低，但数据同化方法仍然对降雨空间模拟结果具有改善作用。该场实际降雨空间分布呈现南北多，中部低的态势，同化前模拟结果中南部地区的降雨量相对较少，和实际分布有一定的差异，北部累计降雨量较大的区域和实际相比相对靠南，没有模拟出东部累计降雨量较少的区域。同化后的结果不仅更加准确的模拟出北部降雨量较大地区的位置，同时也刻画出了流域南部、特别是流域出口处累计降雨量大的特点。

通过三场降雨过程同化前后的模拟结果可以看出，通过数据同化的方法，改善了WRF模式模拟过程中初始场和边界条件的误差，修正了降雨模拟的空间分布结果，使得模拟结果更接近实际的降雨分布情况。

6 结论

由于雷达反射率可以反映出云雨造成的雷达回波功率的大小，直观的刻画大气中云雨的分布情况，有助于模拟大气中的降雨。通过对研究区3场降雨过程同化前后模拟效果的对比分析可以发现：

（1）同化后的模拟降雨数据在空间分布上更符合实测的降雨分布情况，相比于同化前，同化后的模拟降雨数据可以更好地捕捉空间尺度累计降雨量的分布情况，对同化前未模拟到的降雨集中区域有较好的改进作用。此外，CSI指标和RMSE均优于同化前的模拟结果，说明通过数据同化的方法可以提高数值降雨的模拟精度，更好的捕捉不同降雨场次的降雨落区，为大气模式和水文模型的耦合研究提供更加准确的降雨空间分布信息，提高陆气耦合模拟的精度。

（2）降雨场次1和降雨场次2的空间变差系数Cv明显小于降雨场次3，说明降雨场次1和降雨场次2在空间分布上属于相对均匀的降雨类型，而降雨场次3则属于空间分布相对不均匀的降雨类型。从模拟结果中可以看出，数据同化前后降雨场次1和降雨场次2的评价指标远好于降雨场次3的评价指标，从降雨过程可以看出降雨场次3属于短时强降雨，WRF模式和数据同化方法对该降雨类型的模拟效果较差，在今后的研究中可以通过利用多源雷达数据，如增加雷达径向风的同化来改善大气模式对风场的模拟，提高模式对不同类型降雨空间分布的模拟精度。