1 引言
随着模式网格分辨率的增加， 单个格点上降水时间尺度的准确预报变得越来越困难（Lorenz， 1969； Davis et al， 2006a， 2006b）。一些主观上认为在空间位置和发生时间很好的降水预报， 可能会因为缺少观测站点而导致采用传统点对点的检验方法判定为“空报”或者“漏报”的双罚现象（潘留杰等， 2013， 2014a， 2014b； 张宏芳等， 2014）。
面向对象检验MODE（Method for Object-Based Diagnostic Evaluation）是将降水场中超过给定阈值的降水连续区域定义为降水对象， 进而在观测和预报场中进行对象识别和匹配， 检验时通过比较位置、 强度和形状等对象属性， 计算匹配对象的相似度。以往的工作中， Davis et al（2006b）利用MODE方法评估了网格4 km的对流可分辨的WRF中尺度模式在美国中部的预报表现。Johnson et al（2011a， 2011b）利用MODE方法和聚类分析评估集合预报系统降水预报的相似性和差异性。Duda and Gallus （2013）利用MODE方法诊断3 km网格间距WRF模式的对流初生和中尺度对流系统MCS的模拟表现。薛春芳和潘留杰（2016）基于MODE方法研究了日本细网格模式对2013 -2014年中国降水的预报能力， 并与传统检验方法（肖柳斯等， 2019； 杨显玉等， 2020； 张丹琦等， 2019； 石彦军等， 2012； 王颖等， 2010； 张利红和何光碧， 2014）进行了对比。一些工作还基于概率预报系统， 利用MODE方法识别降水对象来对定量降水预报进行修正（Johnson and Wang， 2012）； 比较集合预报不同成员的降水预报技巧（Johnson et al， 2013）。然而， 在已有的MODE方法研究应用中， 仅仅考虑的是两维空间的降水对象， 事实上， 降水对象不仅有二维平面属性， 而且有时间演变特征， 包括时间维的三维空间降水对象， 一方面可以反映二维空间对象的各种属性， 另一方面在表现降水对象的生命周期、 开始和结束时间、 传播速度和演变方面更为直观和有效。
MTD（MODE time-domain）是最近美国DTC（http： //www.dtcenter.org/met/users/downloads/） 在MODE检验方法基础上发展起来的包含时间的三维空间模式降水检验诊断分析工具， 不仅可以检验降水在给定时刻的二维属性特征， 也可以在较长时段内追踪降水对象的时间演变， 目前在降水检验分析中鲜有应用和研究。本文基于ECMWF高分辨率模式网格降水预报资料， 采用MTD方法评估模式在三维降水空间的预报表现， 以期为用户提供新的模式检验分析视角和应用思路。
2 资料来源和方法介绍
2.1 资料来源
一方面考虑到模式可能的调整和升级， 另一方面关注的重点是模式在天气尺度上的预报表现， 因此不选择过长的模式预报资料作为分析对象。资料选择时段为2018年和2019年的6 -8月。ECMWF模式资料选择为每日20:00（北京时， 下同）起报的降水预报， 原始分辨率为0.125° 0.125°。观测资料采用中国气象科学数据信息网CMORPH卫星与全国3万个自动观测站的逐时降水量融合产品（http： //data.cma.cn/site/index.html）， 时段与模式预报一致， 空间分辨率0.05° 0.05°。该产品是采用“PDF+BMA+DS+OI”技术方案研制的， 研究表明融合降水能够较好地表现我国夏季降水的实际分布情况， 可以作为参照物来评估模式的预报性能（潘旸等， 2012， 2018； 沈艳等， 2012）。此外， 为了方便对比， 采用双线性插值法将模式降水预报降尺度变换到0.05° 0.05°。研究范围为秦岭及周边区域（31°N -40°N， 103°E -113°E）， 研究区域的地形及观测站点如图1所示。文中涉及的地图是基于国家测绘地理信息局标准地图服务网站下载的审图号为 GS（2016）1606的中国地图制作， 底图无修改。
2.2 方法介绍
在MTD降水对象识别中， 平滑及降水阈值设定和二维空间检验MODE方法一致（潘留杰等， 2016， 2017）， 同时加入时间维。在给定的降水阈值条件下， 每个时步按照一定的尺度进行空间卷积， 提取对象的空间边界， 其目的是使获得的降水对象所包含的区域更加连续， 滤掉小的或者弱降水对象［图2（a）］。对于三维空间来说， 假定用三维空间内的椭球体降水对象进行拟合［图2（b）］。对于椭球体的三维卷积对象为：
式中： Rs为空间卷积半径； Rt为时间卷积半径； x， y， t分别对应椭球体的 轴。
椭球体内的所有格点降水值通过式（2） 改写：
当降水格点值在椭球体内时值为1， 反之为0。Zx， y， t是时刻为t时， 格点x， y上的降水值。给定卷积半径后， 通过式（2） 得到三维的降水场， 然后用给定的降水阈值， 检查每一个格点逐个时刻的降水值， 如果该格点的降水值达到或超过阈值， 再检查前一个时刻该格点是否属于某个降水对象， 如果属于则该格点和前一时步标记为相同的对象索引。用同样的方法， 判别该格点周围的每个格点， 即可以判别降水三维降水对象的移动、 合并和分裂。时步为时间步长Δt， 本文将降水每3 h的演变定义为一个时步［图2（a）］。
识别对象后， 还采用模糊逻辑来综合各种属性的总体表现， 从而计算预报对象的总收益函数T（a）。式（3） 给出了收益函数T（a）的计算表达式：
式中： ai表示对象的第i个属性； wi表示对象第i个属性的权重系数； Ci（a）表示第i个属性的置信水平； Ii（ai）表示预报场对象的第i个属性的收益函数， 是ai的函数。与二维平面空间对象识别不同， 三维时域追踪还考虑了对象的合并、 分裂以及开始、 结束时间和对象的持续时间。
3 研究结果
3.1 典型降水个例
图3给出了一个基于ECMWF全球模式的三维降水识别追踪个例， 研究的范围为（31°N -40°N， 103°E -113°E）。预报和观测的对应时段为ECMWF模式2019年6月19日20:00起报， 时间间隔为3 h的未来72 h降水， 即24个间隔时次。为了便于表达较强降水的时间演变属性， 降水的强度阈值设置为3 h降水量≥5 mm， 卷积半径设为4个格点， 即10 km。空间三维降水对象的体积最小设置为20个格点。为了表述方便， 图2概念示意图中采用椭球体来拟合降水对象， 但具体计算时采用多边形方法更为精确， 为此， 后续分析中统一采用多边形提取降水对象。从图3可以看出， 模式提取的降水对象相对观测显著偏弱， 在初始时刻（23:00）， 观测场识别了4个降水对象， 分别用数字编号1、 2、 3、 4表示［图3（a）］， 而预报场仅识别了3个降水对象， 也分别表示为降水对象1、 2、 3［图3（b）］。按照匹配规则， 预报场中对象2、 3分别与观测场中对象1、 2匹配。在第二个时刻（02:00）， 观测场中的对象2发生了分裂， 预报场中的对象1减弱消失。到第三个时刻（05:00）， 预报和观测场中同时有新生对象出现。从整体对象的匹配情况来看， 预报场中的对象2相对于观测场中的匹配对象位置偏西， 预报场中的对象3相对于观测场中的匹配对象位置显著偏北， 同时面积显著偏小， 位置少动。
图3 2019年6月19日20:00起报的降水追踪识别的20日23:00（a~b）、 02:00（c~d）、 05:00（e~f）的观测降水（左）和与之对应的预报降水（右）多边形表示降水场中的识别对象, 数字编号标记不同的降水对象 Fig.3 Observed precipitation (left) and corresponding forecast precipitation (right) at 23:00 (a~b), 02:00 (c~d) and 05:00 (e~f) on 20 of the precipitation tracking identification case started at 20:00 on 19 June 2019.The polygons represent the objects in the precipitation field, and the numbers indicate different precipitation objects |
表1给出了预报对象3和观测对象2在整个生命期8个时步（24 h）内的对象属性， 对应的空间位置分布如图4所示， 可以看出预报对象［图4（a）］在所有时步内的位置变化较小， 强降水对象的位置相对稳定， 而观测对象［图4（b）］的位置变化较大， 在质心位置上存在明显的向北推， 然后南压的过程。同时， 观测对象的降水面积也显著大于预报对象。从模式预报相对于观测的系统性时间偏差来看， 前四个时次预报对象的面积相对观测明显偏小， 后四个时次偏大。观测对象在第8个时次已经消散， 预报场中仍然有清楚的表现。从质心位置上来看， 预报的强降水对象中心较观测的明显偏东、 偏南。经向距离偏差接近一个经距， 纬向最大接近2个纬距。从轴角差异来看， 观测和预报对象都有从西南—东北方向逐渐向东西方向转变的趋势， 但预报对象的变化显著偏慢、 幅度偏小。
表1 两维匹配对象的属性特征Table 1 Attributes of two-dimensional matching objects |
| 时步 | 预报对象3 | 观测对象2 | ||||
|---|---|---|---|---|---|---|
| 面积/km² | 质心位置 | 轴角/（°） | 面积/km² | 质心位置 | 轴角/（°） | |
| 时步1 | 301 | 31.64°N, 106.49°E | 48.48 | 881 | 32.71°N, 107.36°E | 61.39 |
| 时步2 | 559 | 31.97°N, 106.89°E | 42.96 | 1336 | 33.09°N, 107.66°E | 52.56 |
| 时步3 | 1362 | 32.04°N, 107.43°E | 39.60 | 1720 | 33.41°N, 108.14°E | 35.25 |
| 时步4 | 1485 | 32.11°N, 107.77°E | 37.89 | 1861 | 33.33°N, 108.50°E | 20.82 |
| 时步5 | 1555 | 32.05°N, 108.04°E | 28.21 | 1258 | 33.07°N, 108.72°E | -3.67 |
| 时步6 | 1041 | 31.91°N, 108.35°E | 10.13 | 455 | 32.70°N, 109.09°E | -10.89 |
| 时步7 | 631 | 31.86°N, 108.45°E | 7.73 | 80 | 32.65°N, 109.01°E | -21.92 |
| 时步8 | 124 | 31.66°N, 108.92°E | 4.76 | / | / | / |
/表示无数据 |
与二维匹配对象属性不同， 三维对象计算整个持续时间内的对象属性， 并计算预报相对于观测的偏差。具体的计算量包括： （1）三维质心距离， 表示不同时步下预报和观测匹配对象质心位置平均值的距离差值； （2）轴角偏差和速度偏差。对不同时步下预报对象采用最小二乘法拟合出对象移动的路径和方向， 使得对象包含的所有空间点到拟合线的距离最小。进而计算出在对象整个时间持续期， 预报相对于观测的速度和方向偏差； （3）体积比率。与传统体积有所不同， 时域追踪时体积的计算方式为所有时步下格点的总和， 所以没有传统体积的量纲单位。体积比率为预报对象格点总数相对于观测的比率； （4）开始结束时间偏差。计算预报对象相对于观测对象在不同时步上的偏差， 譬如预报在第零个时步上识别出对象， 观测在第二个时步上识别出对象， 且这两个对象是匹配的， 则开始时间偏差为-2， 观测滞后预报两个时步， 反之则相反。（5）持续时间偏差。预报对象持续时间相对观测的差值， 负值表示持续时间偏短， 反之偏长。
表2给出了2019年6月19日个例中预报对象（F1~F9）和观测对象（O1~O20）总收益函数大于0.6的对象属性表。从三维匹配对象属性来看， 质心距离偏差最小为12.93 km， 最大为185.51 km， 中位数47.85 km。轴角偏差中位数23.57°， 表明预报和观测对象整体方位一致。从对象移动速度和方向来看， 预报对象平均移动速度较观测偏慢， 平均移速偏差为-0.23。体积比率平均值为0.92， 显示预报相对于观测整体偏小。方向整体和观测一致， 10个对象中有9个预报与观测的方向角偏差小于45°， 仅有预报对象F6和观测对象O9方向偏差158.52°， 对象移动方向相反， 同时该匹配对象的总收益函数也相对较低， 仅为0.614。对象的开始与结束时间显示， 预报与观测开始时间一致的有4个， 偏早的有3个， 偏晚的有2个， 整体比较吻合。结束时间中位数为 -1， 整体表现为结束时间偏早。10个匹配对象中， 2个持续时间完全一致， 6个持续时间偏差为负值， 表明预报对象的生命期相对于观测整体偏短。
表2 三维匹配对象的属性特征Table 2 Attributes of 3D matching objects |
| 匹配 对象 | 质心距离 /km | 轴角偏差 /（°） | 速度偏差 /[km·(3h)-1] | 方向偏差 /（°） | 体积 比率 | 开始与结束时间 的时步偏差/个 | 持续时间的 时步偏差/个 | 总收益 函数 |
|---|---|---|---|---|---|---|---|---|
| F1-O6 | 33.92 | 23.74 | -0.14 | 0.0 | 0.92 | -2, -4 | -2 | 0.839 |
| F2-O1 | 12.93 | 15.67 | 5.69 | 36.34 | 1.62 | 0, -1 | -1 | 0.766 |
| F3-O2 | 24.51 | 9.24 | -0.77 | 4.58 | 0.93 | 0, 1 | 1 | 0.899 |
| F4-O6 | 44.93 | 23.574 | -0.14 | 0.0 | 0.65 | 0, -2 | -2 | 0.881 |
| F5-O20 | 50.70 | 27.19 | 0.09 | 41.04 | 1.14 | -12, -12 | 0 | 0.736 |
| F6-O9 | 80.30 | 1.70 | -3.20 | 158.52 | 0.25 | 0, -1 | -1 | 0.614 |
| F7-O8 | 185.51 | 27.10 | 0.60 | 11.06 | 0.93 | 2, 2 | 0 | 0.725 |
| F8-O12 | 64.07 | 29.37 | 0.46 | 28.98 | 0.71 | -2, -3 | -1 | 0.746 |
| F9-O10 | 21.89 | 23.28 | -4.69 | 0.27 | 0.78 | 1, -1 | -2 | 0.862 |
3.2 统计结果
单个个例分析显示， ECMWF模式预报相对于观测， 无论是二维对象或者是三维匹配对象的属性偏差， 尽管不存在清晰的、 一致的、 精确的统计指标， 但可以发现一些相对较为明显的定性统计特征。例如， 2019年6月19日起报的个例中， 二维匹配对象经纬度整体特征表明模式预报对象位置较观测显著偏南、 偏东； 在三维对象上体积显著偏小， 结束时间偏早， 持续时间偏短、 移速偏慢。但单个个例不具有普遍性， 为此， 利用2018年和2019年6 -8月降水预报资料， 分析模式预报的整体统计特征。
从不同阈值下识别的二维降水对象的空间分布特征（图5）可知， 整体来看， 不同阈值下降水对象数量的纬向型特征分布非常明显， 呈现出明显的南高北低的特征， 在所有阈值下， 33°N以南为一致的降水对象高频次发生中心， 然后向北逐渐递减， 但低降水阈值［图5（a）， （b）］与高阈值降水［图5（c）， （d）］对象的发生频次在34°N附近的秦岭山脉附近存在显著差异， 0.1 mm和0.5 mm的低降水阈值在秦岭山脉附近有较为显著的高频次降水对象触发中心， 但在高阈值， 特别是5 mm降水阈值［图5（d）］情况下， 秦岭山脉为一致的低频次降水对象触发中心。对比预报和观测对象频次的空间分布， 可以看到在低阈值降水情况下［图5（a）］， 两者发生频次的高低值中心比较一致， 陕西西南部、 山西北部的高值中心以及关中盆地、 宁夏西北部的低值中心预报和观测都非常吻合， 但随着降水阈值的增大， 预报和观测对象发生频次的差异显著增大。降水阈值0.5 mm时［图5（b）］， 观测在甘肃西南部呈现出降水发生频次低值中心， 但预报频次的低值中心表现并不明显。随着降水阈值的增大， 两者的空间分布形态差异更为显著， 位于陕西关中的降水对象低频次发生中心与地形吻合较好， 但从预报来看， 低值中心显著偏北、 范围偏大。此外， 观测场上5 mm降水阈值在宁夏东部、 陕西西部的高阈值中心在预报场上整体没有表现。总体来看， 模式预报的强降水对象的质心位置比观测偏大， 弱降水对象与观测基本一致。
图6给出了不同降水阈值和卷积半径下降水对象的生命周期。整体来看， 预报对象较好地表现了观测对象的生命周期变化特征。在不同卷积半径条件下， 无论是低降水阈值或高降水阈值， 降水对象的生命周期呈明显的偏态分布， 超过80%降水对象的生命周期为2~5个时步， 即持续时间为6~15 h， 预报和观测表现基本一致， 但在较短生命周期条件下， 预报对象的数量较观测明显偏少。对比同一卷积半径下， 不同降水阈值下对象的持续时间， 可以发现当卷积半径较小时， 观测降水的低降水阈值［图6（a）］对象数量低于较高阈值［图6（b）］， 降水对象数量偏少， 即中等降水阈值下识别的降水对象的数量最多， 然后当降水阈值增大时［图6（c）］对象数量显著下降， 模式预报识别降水对象整体较观测偏少， 在中等降水阈值识别的偏少最为明显。
以往通过传统检验方法 （张宏芳等， 2014）认为， 模式对小雨预报显著偏多， 但本文研究结果， 模式预报对象较观测明显偏少， 分析发现， 这两者并不矛盾， 主要原因在于在低阈值条件下， 模式由于小雨预报偏多， 降水对象往往表现出面积大的特征， 整体使得降水对象的数量下降。
三维对象的统计属性如图7所示， 从对象质心距离来看， 6月降水个例质心距离整体偏大， 最大超过180 km， 平均为55 km。进入7月以来， 质心距离显著减小， 最大为144 km， 平均减小为44 km， 表明从平均空间位置来看， 7月和8月预报的比6月更好。持续时间偏差［图7（b）］显示， 6月主要为负偏差， 预报对象持续时间比观测显著偏短， 平均偏短0.5个时步（1.5 h）。进入7月， 预报对象持续时间整体转好， 最大相差3个时步（9 h）， 平均偏短0.15个时步， 约偏短0.45 h。8月预报对象持续时间仍然较观测偏短1.3个时步（3.9 h）。因此， 整体看来， 预报对象的持续时间较观测偏短。
从速度偏差来看［图7（c）］， 6月到7月上旬， 预报对象的移动速度较观测整体偏慢， 从统计结果上来看， 最大偏慢量为7.17格点·（3h）-1， 在MTD中， 移动速度偏差是通过不同时步下预报和观测识别对象质心位置的线性拟合， 进而根据空间分辨率和时步而计算出来的。在本研究中， 网格降水的空间分辨率为5 km， 3 h为一个时步， 可以得到移动速度偏差最大偏慢11.95 km·h-1。从7月中下旬的移动速度来看， 预报和观测对象整体一致， 平均移动速度量值为0.0668， 预报较观测对象移动速度略偏快， 为0.11 km·h-1。8月预报对象的移动速度又转为偏慢。从体积比率［图7（d）］来看， 模式预报对象的整体比观测对象偏大， 平均比率为观测的1.2倍。事实上， 由于预报对象相对观测的持续时间整体偏短［图7（b）］， 而MTD时域追踪时体积的计算方式为所有时步下格点的总和， 所以单个时步下， 预报对象相对观测对象大于5 mm的降水面积可能会更大一些。
开始结束时间偏差显示， 两者整体具有一致性， 从6月至7月上旬， 预报对象与观测对象的时间偏差整体为正值， 即预报对象相对观测的开始和结束时间整体偏晚， 到7月中下旬转为偏早， 至8月又转为偏晚。平均来看， 预报对象开始时间较观测偏晚0.23个时步（0.68 h）， 结束时间较观测偏早0.05个时步（0.15 h）。
4 结论和讨论
基于ECMWF模式2018年和2019年降水预报资料， 中国气象科学数据信息网CMORPH卫星与全国3万个自动观测站逐时降水量融合的实况分析产品， 利用三维时域追踪分析方法MTD， 研究ECMWF模式对秦岭及周边区域降水的预报表现， 主要结论如下：
（1） 逐3 h降水个例追踪显示， MTD能够很好地从降水场中识别降水对象， 并追踪不同时步下降水对象的属性特征， ECMWF模式预报较好地表现了观测降水的强度、 面积和时间偏差特征， 但模式对强降水预报位置整体偏南， 预报面积偏小， 强度偏低， 对象持续时间相对观测偏短。
（2） 降水对象的空间分布有明显的地形影响， 这种情况在高降水阈值条件下更为显著。模式预报很好地描述了降水对象的空间分布， 低阈值条件下预报和观测对象触发频次的高低值中心比较一致， 但随着降水阈值的增大， 预报和观测对象发生频次的差异显著增大。随着降水阈值的增大， 预报与观测降水对象的空间频次呈现出显著差异， 表明模式对强降水的位置预报仍然需要改进。
（3） 预报对象较好地表现了观测对象的生命周期变化特征， 在不同卷积半径条件下， 无论是低降水阈值或高降水阈值， 降水对象的生命周期呈明显的偏态分布， 超过80%的降水对象的生命周期为2~5个时步， 即持续时间为6~15 h， 预报和观测表现基本一致。但在较短生命周期条件下， 预报对象的数量较观测明显偏少， 分析发现， 主要原因在于在低阈值条件下， 模式出现过多微量降水， 降水对象往往表现出面积大的特征， 整体使得降水对象的数量下降。
（4） 三维对象追踪显示， 无论是质心距离、 体积比率、 持续时间， 还是开始和结束时间， 在较长时段内模式都没有一致的时间偏差预报表现。两年的平均统计结果显示， 预报对象的持续时间较观测偏短， 移动速度较观测整体偏慢， 预报对象开始时间较观测偏晚0.23个时步（0.68 h）， 结束时间较观测偏早0.05个时步（0.15 h）， 但在不同的时段呈现出显著差异， 6月降水预报质心距离整体大于7月和8月， 速度偏慢， 持续时间偏短， 开始时间偏早， 结束时间偏晚。
事实上， 在长时段内预报对象相对于观测没有一致的时间偏差， 这一方面是由于模式本身的升级和调整， 另一方面， 模式在不同月季尺度上的预报表现显著不同， 因此， 过长时段的统计结果对把握天气尺度的降水预报意义不大， 但在相对较短时间内， 预报对象的时间偏差具有明显的倾向性， 利用MTD检验方法获得定量化的时间偏差， 这对订正降水的时间预报偏差具有重要意义。本文的时步设置为3 h， 时步设置对检验结果可能会有一定影响， 后期将通过调整时步的设置， 开展更为细致的检验。