多普勒天气雷达是目前探测中小尺度对流天气系统的重要观测工具。美国自20世纪80年代开始建立多普勒天气雷达观测网, 为监测对流天气发挥了重要作用, 大大提高了灾害天气预报预警的准确率。我国从20世纪90年代末开始建设新一代多普勒天气雷达观测网, 为提高我国气象业务水平发挥了重要作用。但是, 目前国内雷达都是单站观测, 产品也基于单站, 预警范围和临近天气预报外推时间都受到较大的限制。国内已开发了基于雷达二次产品的拼图软件^[1], 并投入业务应用, 但其拼图产品单一, 预报员只能被动接收, 无法进行二次开发, 无法得到天气系统的三维结构等重要信息。为了提高雷达对中尺度灾害性天气预报预警能力, 有效利用雷达观测网的优势, 有必要多部雷达联合观测, 并进行多部雷达观测结果的拼接。雷达三维组网数据不仅具有很高的分辨率, 而且覆盖范围广, 可以用于强天气分析算法、区域降水估计、中小尺度数值天气模式中雷达数据同化等^[2-5]。美国国家强风暴实验室开发了全国范围内高分辨率的雷达数据三维组网系统, 于2004年投入业务试运行^[6-7]。该系统在强天气数值模式数据同化和降水估计等方面均发挥了重要作用。美国的人工影响天气工程中也纳入了雷达三维组网结果, 为人工影响天气提供支持。经多年试验和完善, 该系统已于2006年正式交给相关业务部门应用。

刘黎平等^[8]对雷达资料质量控制进行了深入研究, 取得了很好的效果。肖艳娇等^[9]探讨了雷达三维数字组网的多种方法, 分析了其业务应用能力和可行性。基于这些科研成果, 中国气象科学研究院灾害天气国家重点实验室有关人员在国内率先开发了多普勒天气雷达三维数字组网软件系统, 对区域内多部天气雷达的反射率因子数据进行三维数字组网, 并以组网数据为基础, 进一步处理得到一些雷达二次产品。

多普勒天气雷达三维数字组网系统只负责数据处理和产品生成, 为方便应用, 同时还开发了雷达三维数据显示系统。为避免相互争抢系统资源, 根据微机性能, 可将两个软件安装在同一台微机上, 或者分别运行于不同的微机上, 通过网络实现组网结果数据的共享。组网处理系统的核心功能就是雷达数据组网, 输出结果数据, 不需要用户的频繁干预, 因此界面功能简单便捷。组网过程中数据处理量大, 实时性要求强, 因此, 数据处理任务采用多线程实现。操作界面在必要时才显示, 操作完成后可隐藏, 实现系统后台运行。

多普勒天气雷达网由多型号(SA,SB,CB,CD,SC,CC)、多波段(S,C)的多普勒天气雷达组成, 且数据格式不统一, 因此, 本系统提供设置雷达站的位置、高度和雷达型号等参数接口, 并在预处理过程中首先将不同型号雷达的体扫数据转换为统一的格式, 然后对统一格式的数据进行质量控制等处理。

系统兼顾科研和业务两种需求, 既能实时运行, 又能分析处理历史资料; 还可选择输出一些中间数据, 如质量控制后的体扫数据, 单站三维格点反射率因子数据; 输出的组网结果可以多种格式保存。

本软件系统采用C++和Fortran语言开发, 定位于15部左右雷达的区域组网, 在以Windows XP为操作系统的微机上运行。在雷达资料控制的基础上, 系统完成省或区域中心的多普勒天气雷达基数据组网, 实时提供经、纬度及海拔高度网格的回波强度资料, 并生成部分二次产品, 为中小尺度天气分析、临近天气预报、人工影响天气等应用提供基础。

系统包括控制层和后台数据处理层两大部分。其中, 控制层为软件界面, 是系统和用户实现交互的平台; 后台数据处理层为软件的核心部分, 由一系列的功能模块构成:雷达基数据监视和时间匹配模块, 基数据质量控制模块, 体扫数据从雷达极坐标向以经、纬度、海拔高度为单位的笛卡儿坐标转换模块, 组网模块以及多种二次产品生成模块。各种功能模块分别由不同的线程来完成, 线程间通过消息传递数据, 通过同步控制来协调各功能模块的执行。

多普勒天气雷达三维数字组网系统的主要工作是基于后台的数据处理, 要求界面简洁, 操作方便。因此, 软件采用了托盘图标结合窗口界面的方式, 在完成必要的设置和操作后, 可关闭界面窗口, 只在操作系统任务栏上显示一个托盘图标, 托盘图标是动态的, 不同图标表明系统处于不同状态。进行常用操作时不需要打开界面窗口, 只需通过托盘图标上的弹出菜单来选择。操作界面含有多属性页对话框, 其功能是用于设置软件系统必要的各类参数。

组网系统有两种数据处理模式:实时模式和事后分析模式。在实时模式下, 系统监视指定位置的实时观测数据, 检索可参与当前时次组网的数据; 在事后分析模式下, 系统一次性载入一定数量的各站体扫资料, 依据预先设置的时间间隔, 逐时次处理这些数据。体扫数据经过网络传输到指定中心, 网络传输速度可能造成数据延时到达区域中心, 为此, 实时处理时, 可设置时间延迟, 使系统处理当前时间之前体扫数据; 因不同区域雷达观测资料获取方式有差异, 系统或其附加工具提供多种载入雷达基数据的途径, 如:访问本地磁盘上的数据、通过ftp下载等, 还提供基数据实时解压功能。

系统监视到基数据后, 首先要对其进行时间匹配, 判断该数据是否参与本时次的组网。例如, 设置组网时次间隔为6 min, 匹配时间为组网时刻的前后2 min, 则各组网时次依次为00:06,00:12,00:18, ……。其中, 00:06这个时次可参与组网的观测数据的时间范围为00:04至00:08, 这个时间范围内的各站数据逐个处理后再组网, 形成00:06的组网结果, 其他时次依此类推。在本系统中, 组网时次间隔、匹配时间等都作为可调参数, 可以通过界面修改。

地物是影响雷达数据质量的重要因素, 雷达反射率因子数据在输入其他算法之前, 有必要对其进行质量控制。本系统在读入雷达基数据后, 采用基于模糊逻辑的分步式超折射地物回波识别方法进行地物识别和剔除, 同时对孤立回波等进行剔除, 保证数字组网的质量。图 1给出了2005年6月21日11:54(北京时, 下同), 天津SA雷达体扫数据进行地物回波和孤立回波识别处理前后, 0.5 °仰角的回波对比。结果表明:该方法对地物回波和降水回波区分效果较好, 大片地物被有效滤除, 甚至某些距离模糊区, 速度不可用的位置, 地物也被识别出来^[8]。

图 1

图 1. 0.5°仰角质量控制前(a)、后(b)雷达反射率因子和速度对比 Fig 1. Contrast of reflectivity and velocity without quality control (a) and with quality control (b) at elevation of 0.5°

雷达体扫数据经过质量控制后, 再从极坐标系转换到笛卡儿坐标系。格点插值有4种方法可选^[10-14], 本系统的单站雷达数据坐标转换采用NVI方法 (径向、方位上的最近邻居法结合垂直线性内插法)。根据肖艳娇等^[9]综合分析结果, 采用NVI方法, 处理得到的反射率因子格点数据不仅空间连续性好, 而且很好地保留了体扫资料中原有的结构特征。

笛卡儿坐标的高度层和经纬度格点分辨率可以通过操作界面进行设置, 垂直方向最多可设40个高度层。系统默认的水平分辨率为0.01°×0.01°, 在中纬度地区约1 km左右; 根据天气系统的高度分布特征, 垂直分辨率预设为随高度增加而降低, 默认的分辨率低层为500 m, 中层为1 km, 高层为2 km。

利用2005年5月1日08:02的福州雷达体扫资料, 分别用SA雷达产品生成软件和本组网软件处理得到组合反射率因子及反射率因子垂直剖面, 图 2是两个软件处理结果的对比。SA雷达产品生成软件处理结果是直接从原体扫数据计算得到的, 可以认为是原始数据信息; 组网软件系统处理结果, 是根据坐标转换后的三维格点数据处理得到的。比较两个软件生成的组合反射率因子图 2a和图 2b, 可以看出, 回波强度的分布基本一致, 但由于插值处理对反射率因子有一定的平滑作用, 图 2b相对于图 2a, 有些位置的回波强度稍有减弱, 但差别不大, 在5 dBz以内, 不影响系统的定性分析; 比较图 2c、图 2d两个剖面图, 可以看出, 回波的基本结构和分布是基本一致的, 且三维格点数据因为经过了插值处理, 图 2d较图 2c平滑。对比结果表明, 经过坐标转换后的三维格点数据很好地保留了原始数据的信息, 能够体现系统的三维结构。

图 2

图 2. 由SA雷达软件RPG和三维组网软件得到的组合反射率因子及垂直剖面对比(R=240 km) (a)SA雷达软件RPG根据体扫资料生成的组合反射率因子(黑线标记处为剖面位置), (b)三维组网软件根据三维格点雷达资料生成的组合反射率因子, (c)SA雷达软件RPG根据体扫资料生成的反射率因子垂直剖面, (d)根据三维格点雷达资料生成的反射率因子垂直剖面 Fig 2. Contrast of compositer eflectivity and vertical cross section from volume scan and 3-D mosaic radar data(R=240 km) (a) composite reflectivity from RPG (black line is the position of cross section) , (b) composite reflectivity from 3-D digital mosaic system, (c) cross section from RPG, (d) cross section from 3-D mosaic data

当前时次可参与组网的全部单站数据都完成了坐标转换后, 启动组网任务, 将各单站的格点数据拼到一个大范围的笛卡儿坐标网格中。这个过程采用一种指数权重的算法^[9]。

肖艳娇等^[9]分析指出, 在多部雷达组网处理的多种算法中, 指数权重法针对雷达资料具有明显优势, 指数权重函数的形状易于调整以便达到随距离快速下降的效果, 同时又保证了权重值是一个正值。指数权重随距离快速下降, 可以保留原始雷达资料中近距离处的高分辨特征。对于远距离来说, 正的指数权重又可以确保雷达能影响到它所覆盖的整个区域。

2007年6月9日广东省有一次大范围降水过程, 图 3是该过程中08:06的组网反射率因子和单站格点反射率因子。参与组网的有广州、汕头、梅州、阳江、韶关、深圳6部雷达。图 3a是5km高度上雷达反射率因子的组网回波图, 从图中可见一片大范围的降水回波; 图 3b, c,d,e分别是梅州、深圳、韶关、汕头单站5 km高度回波图。首先, 任何单站都只能看到有限范围的雷达回波, 组网后, 除看到当地雷达回波, 还可以看到其他雷达回波, 这可为中小尺度天气分析提供便利。此外, 图 3b, d,e中, 雷达站附近几公里范围位于体扫的盲区, 这些位置都存在回波缺失的情况, 但在组网处理后, 由于邻近雷达可以观测到这些盲区, 使得这些区域回波得到填补; 图 3b中, 梅州站东南方向, 距离雷达站大约190 km以外的回波受地形阻挡影响严重, 看不到回波, 而从图 3a中黄色圆圈标记的梅州站范围看, 被阻挡处的回波得到填补, 这是因为邻近的汕头雷达站可以观测到该处, 通过组网处理, 汕头雷达的观测资料可以填补该处空白; 图 3c中, 深圳站东北方向, 距离雷达站大约200 km的回波也受地形阻挡, 从图 3b, e可看到, 梅州、汕头两部雷达都可以探测到该阻挡处的回波, 组网处理时, 这两个雷达站的观测资料对该处有补偿作用, 在图 3a中橙色圆圈标记的深圳站范围内, 可以看到, 该站受阻挡处被其他站的回波填补。

图 3

图 3. 5 km高度雷达反射率因子数据组网前后效果对照 (a)组网, (b)梅州雷达, (c)深圳雷达, (d)韶关雷达, (e) 汕头雷达 Fig 3. Contrast of reflectivity at the height of 5 km between single radar and mosaic (a) mosaic reflectivity, (b) reflectivity of Meizhou radar, (c) reflectivity of Shenzhen radar, (d) reflectivity of Shaoguan radar, (e)reflectivity of Shantou radar

除组网功能外, 本软件系统还包含了一些二次产品处理模块, 这些模块以组网基数据为输入, 经过相应的算法处理, 输出雷达气象产品, 包括回波顶算法及产品生成模块, 垂直累积液态水含量算法及产品生成模块, 组合反射率因子算法及产品生成模块, 回波外推算法^[15-17]及产品生成模块等。

系统可输出雷达二次产品以及各级数据。主要产品为三维组网反射率因子数据, 以及在其基础上生成的组合反射率因子、回波顶高度、垂直累积液态水含量、回波移向移速产品和回波位置外推产品; 中间数据有质量控制后的基数据、单站三维反射率因子数据。

本文中组网图像及回波垂直结构由多普勒天气雷达三维显示系统绘制输出。它是组网系统的配套系统, 主要功能包括:实时监视组网系统新生成的数据和产品、显示各层组网回波、显示组合反射率因子、回波顶等各种二次产品、显示组网反射率因子的三维结构、显示任意位置的垂直剖面、背景地图叠加显示、获取指定回波的位置信息、多时次产品及回波动画、图像自动实时更新、图像放大、缩小、图像存储等。张志强等^[18]对显示系统的结构和功能、系统开发等情况有更详细的描述。

目前, 本系统已经在北京市气象局2008年奥运气象保障中进行了近半年的连续试运行, 运行稳定, 其输出的雷达三维组网数据供北京市气象局业务应用, 组网结果可靠。同时, 作为“我国南方致洪暴雨监测与预测的理论和方法研究”的重要成果推广, 该系统也被用于广东省气象局业务运行, 其三维组网结果及二次产品输出到统一的数据存储服务器, 作为中尺度灾害天气分析与预报平台产品数据源的一部分。此外, 本系统及其配套的雷达三维显示系统还在河北省气象台进行业务试用, 提供的大范围雷达回波以及二次产品, 避免了以前预报人员需要同时打开多个雷达显示软件查看周边天气信息的不足, 为天气分析提供了便利条件。

测试时采用的计算机为双核处理器、主频2.13 GHz、内存2G。三维组网系统在该机器上模拟业务运行, 区域内有12部雷达参与组网, 平均5 min完成一次组网处理及二次产品生成, CPU占有率平均不超过50%, 最大内存占用量102 MB。在北京市气象局服务器上业务试运行, 同步观测的4部雷达参与组网 (3部SA雷达, 1部CB雷达)。根据体扫周期, 组网频率为每6 min一次, 连续运行1个月以上没有出现系统退出现象。模拟测试和业务试运行结果表明:对于包含15部左右雷达的区域组网, 只需1台性能适合的微机, 可以实现周期为5~6 min (新一代多普勒天气雷达的体扫周期为5~6 min) 的三维数字组网, 且可以长期不间断地稳定运行, 满足业务需求。

2007年8月19日圣帕台风在福建登陆, 10:40台风已经减弱为热带风暴, 将该时间的雷达组网结果和全国拼图进行比较, 参与组网的为福州、厦门、建阳、龙岩、梅州、汕头和赣州7部雷达。

目前的雷达业务应用中, 针对热带风暴及更大尺度的天气系统, 除全国业务拼图产品外, 还没有更精细的产品能体现其结构。全国雷达拼图系统可以提供的雷达回波强度的产品是0. 5°仰角基本反射率因子拼图和组合反射率因子拼图, 0.5°仰角基本反射率因子拼图的回波不在同一高度上, 可比性差, 因此, 这里选取组合反射率因子进行比较 (图 4)。图 4a是全国雷达拼图业务系统对单站组合反射率因子产品的拼图结果, 图 4b是本系统在雷达三维组网格点数据基础上生成的组合反射率因子, 排除色标差异, 两幅图中的回波分布一致。组网系统提供的主要结果是高分辨率的反射率因子三维数据, 在这个数据的基础上, 不仅可以得到组合反射率因子、多种高度层上的基本反射率因子, 还可以显示回波的垂直分布、三维结构。图 4b中左上角的小窗口中是沿图中红色直线的反射率因子剖面, 显示出该处有发展很强的降水回波; 左下角小窗口中是红色矩形区的三维回波切面, 显示沿该矩形两个边界的回波垂直剖面和不同高度层的回波平面分布。因此, 雷达三维组网数据不仅覆盖范围广, 能反映大中尺度系统信息, 而且空间分辨率高, 并包含垂直回波分布信息, 有利于进行中小尺度天气系统分析。此外由于全国业务拼图结果是基于图像的拼接, 数据只有16个级别, 而三维组网数据具有与基数据同样的精度, 如果用来进行二次开发, 前者精度太差, 只能勉强使用, 甚至不可用。

图 4

图 4. 雷达三维数字组网结果和全国业务拼图结果对比 (a) 三维组网系统生成的组合反射率因子, (b) 来自全国雷达拼图业务系统的组合反射率因子 Fig 4. Results of 3-D digital mosaic system and national radar mosaic system (a) composite reflectivity generated by the 3-D digital mosaic system, (b) composite reflectivity from national radar mosaic system

雷达定位的准确性不仅是关系到组网结果正确与否的重要因素, 也直接影响单站雷达产品的业务应用。组网软件将雷达体扫数据通过插值转换到以经度、纬度、海拔高度为单位的笛卡儿坐标系, 对于不同雷达探测的同一块回波, 经过插值后, 其经、纬度位置应该是一致的。通过分析多个雷达观测到的同一块回波在组网后的形态和单站回波之间的差异, 或者假设某一雷达在不同位置观测该处回波, 分析其在不同位置时, 与其他邻近雷达对同一位置回波的观测值之间的相关性, 就能判断雷达定位的准确性。组网系统输出的单站三维格点反射率因子数据和三维组网反射率因子数据为这种分析提供了便利。

2007年7月11日00:06, 北京、河北、天津部分地区有强降水过程, 在用SA雷达PUP软件对比分析多部雷达回波的过程中, 发现同一块强回波 (图 5a中圆圈标记位置), 北京站显示的位置与天津、秦皇岛站显示的位置不一致, 距离偏差至少15 km。为分析这个问题, 对这3部雷达的观测资料进行组网, 同时选择输出单站格点资料。图 5b显示的是组网后5 km高度的回波。首先, 从组网回波图可直观的看出, 组网后回波形状和每个单站显示的回波形状 (图 5c, d,e) 是一致的, 因此基本可以排除北京雷达回波的定位错误。其次, 取图 5b, c,d,e中矩形区域内的回波, 分析北京雷达与其他两部雷达的回波相关性。

图 5

图 5. 北京 、天津 、秦皇岛雷达回波相关性分析 (a)北京 、天津 、秦皇岛SA雷达软件显示的 1.5°仰角反射率因子, (b) 3 km 高度组网反射率因子, (c)北京 3 km 高度反射率因子, (d)天津 3 km 高度反射率因子, (e)秦皇岛 3 km 高度反射率因子 Fig 5. Relativity between the echo of Beijing , Tianjin and Qinhuangdao radar (a) reflectivity at 1.5 °elevation produced by SA radar software of Beijing, Tianjin and Qinhuangdao, (b) mosaic reflectivity at 3 km height, (c) reflectivity of Beijing at 3 km height, (d) reflectivity of Tianjin at 3 km height, (e) reflectivity of Qinhuangdao at 3 km height

表 1, 2分别是北京雷达与天津、秦皇岛雷达之间的相关分析结果。很明显, 在北京雷达的回波不进行任何移位情况下, 与其他两部雷达的相关性最大 (均大于0.8), 且远大于移位后的相关性。分析结果表明:北京雷达的回波定位是正确的, 应是北京站显示的地理信息有误, 导致图上回波相对于地理信息产生位置偏差。

表 1

表 1 北京与天津雷达回波相关分析 Table 1 Correlation of radar echo between Beijing and Tianjin

表 1 北京与天津雷达回波相关分析 Table 1 Correlation of radar echo between Beijing and Tianjin

表 2

表 2 北京与秦皇岛雷达回波相关分析 Table 2 Correlation of radar echo between Beijing and Qinhuangdao

表 2 北京与秦皇岛雷达回波相关分析 Table 2 Correlation of radar echo between Beijing and Qinhuangdao

组网区域内各雷达标定正确, 系统误差尽可能小, 才能保证雷达观测数据的可靠性。目前发现个别雷达观测误差较大, 观测值和组网区域内其他雷达相比有明显差异, 这些观测误差会影响组网结果的质量。为此, 利用组网系统的中间数据, 即单站三维格点反射率因子资料, 取相邻两个雷达之间距两个雷达相近的回波, 分析它们之间回波强度的差异。

这里取2007年8月19日13:00-19:00, 福建省福州、建阳、龙岩、厦门4部雷达的观测资料, 转换为三维格点资料。相邻雷达比较时, 距两部雷达的距离之差在5 km以内的回波, 近似认为位于两部雷达的中点位置, 取中点处5~10 km高度范围内的两部雷达的反射率因子值进行对比, 计算反射率因子差值的格点平均和时间平均。

(1)

式(1)中, D_ref是两个雷达反射率因子差的平均, r₁, r₂分别是两个雷达的反射率因子观测值, T是观测次数, N是格点数。

对两相邻雷达同一回波的观测值进行上述计算, 建阳站比厦门站小2.1 dBz, 比福州站小2.4 dBz, 比龙岩站小3.1 dBz, 其他站之间的观测差异不明显, 这说明建阳雷达相对其他3部雷达, 观测值偏低大约2~3 dBz。

鉴于上述情况, 这4部雷达进行组网, 在建阳雷达与其他雷达共同覆盖的区域, 组网后反射率因子值会偏低, 且建阳雷达单独覆盖的区域, 反射率因子值偏低程度最大。

多普勒天气雷达三维数字组网系统在国内首次开发并投入业务试运行。该系统应用了中国气象科学研究院灾害天气国家重点实验室近年来在天气雷达应用方面的最新科研成果。系统包括基数据时间匹配和载入、质量控制、单站体扫数据坐标转换、多个雷达站笛卡儿坐标数据组网、二次产品生成等一系列处理过程, 最终输出产品包括:质量控制后的基数据、单站三维格点反射率因子数据、三维组网反射率因子数据和基于组网反射率因子得到的二次产品。

1) 与现有的业务拼图相比, 多普勒天气雷达三维数字组网结果包含回波的三维信息, 数据精确度更高, 这既有利于天气分析, 还可为临近预报二次开发和应用、中小尺度数值天气模式中雷达数据同化等提供数据支持。

2) 模拟试运行和业务试运行结果表明:系统运行稳定, 数据的时间、空间分辨率高, 符合业务需求。

3) 实例分析和业务试用结果表明:组网结果具有单站雷达数据不可比拟的优势, 且结果可靠, 对中小尺度灾害天气系统业务分析提供一定的帮助。

4) 通过分析组网系统的输出结果可确定雷达定位的准确性。

5) 有些雷达观测值偏差明显, 可能影响组网效果。利用组网系统的数据结果, 通过雷达观测误差分析, 可确定不同雷达之间的观测差异。在系统的升级版本中有必要增加误差订正功能。

本组网系统目前处于业务试运行阶段, 结合业务试用情况和新的科研成果, 下一步准备从以下几个方面对软件进行完善。

1) 由于系统标定、强降水衰减程度以及不同波段雷达性能等差异, 导致不同雷达对同一块降水云的观测值可能存在差异。为此, 准备在软件系统中增加雷达观测值订正功能。

2) 三维组网雷达和自动雨量站降水估算方法已经过个例研究, 将在软件系统中增加降水产品, 并增加雨量计订正功能。

3) 有些雷达受阻挡严重, 经过组网, 邻近雷达在一定程度上可以弥补被阻挡区的反射率因子值, 但采用考虑权重插值方法, 相对于实际观测值, 这些弥补回波强度会不同程度地减弱, 进而导致根据反射率因子得到的降水估算值会产生较大偏差。针对这种情况, 将增加混合扫描处理方法, 弥补其不足。

4) 单站雷达数据坐标转换处理增加另外3种插值方法, 应用中可以根据实际情况选择不同算法。

5) 以高性能微机或工作站为平台, 进一步提高软件系统运行效率, 实现更大区域的雷达数据组网。