引言

双偏振雷达是监测、预警突发性灾害天气最有效的手段之一，在监测雷暴、冰雹、台风、暴雪等强天气系统中发挥着重要作用^[1]。根本上讲，气象雷达探测到的数据是3D结构，但由于目前雷达产品呈现主要针对二维栅格图像的点、线、面矢量符号，缺乏立体空间的垂直连贯信息，无法展现垂直方向上气象要素的分布特征^[2-3]，故开展双偏振雷达的三维可视化研究具有重要意义。

国内外，针对雷达探测数据处理分析和立体可视化应用技术有很多研究。美国LAKSHMANAN et al.^[4]对多部天气雷达的三维可视化组网进行了初步研究与探索，RU^[5]对美国多普勒天气雷达进行了数据解析与回波立体绘制研究。国内孟昭林^[6]对天气雷达数据与地理信息立体结合进行了探索，范大伟等^[7]对机场多普勒天气雷达三维显示平台进行了研究开发。这些研究在立体回波透视程度上大都效果不佳，且仅对反射率产品进行了探索，未对双偏振产品展开深入研究。

本文在双偏振雷达所采集的反射率因子数据基础之上施以体绘制算法，首先对原始数据进行解析、质量控制与坐标变换等预处理，再通过立体建模形成三维数据场作为体绘制的插值初值场，最后通过插值法进行三维云体绘制，使之最终以三维透视可视化的产品形式展示，该立体回波模型可应用于基于地理信息系统(GIS)的3D雷达产品开发，也可为预报人员剖析云体内部特征提供重要依据。

1 天气雷达双偏振产品三维场建立

天气雷达双偏振产品三维场的建立主要包括双偏振数据地物杂波质量控制，坐标变换，三维网格建立与插值等数据场建立过程。双偏振雷达采集到的双偏振数据在距离雷达中心较近时会受到地物遮挡影响产生杂波^[8]，必须排除地物杂波的干扰，并且雷达回波数据是以极坐标的方式保存，而三维场是以直角坐标的形式保存^[9]，这就需要对双偏振原始数据进行坐标变换与三维插值处理。

1.1 地物杂波抑制

地形遮挡会严重影响回波质量，尤其近场回波干扰严重。在三维场绘制中，地物杂波干扰会影响近场数据的准确性并破坏整体数据纵向拉伸性^[10-11]。由于地物杂波的平均径向速度产品特征为一个接近零速度的区域中分布着零星非零速度值，利用地物杂波径向速度趋近于零的特性，抑制掉地物杂波的干扰^[12]。

1.2 坐标变换

双偏振产品原始数据为同一体扫内多个离散仰角的锥面扫描数据集，是以雷达馈源为中心的仰角、方位和距离为参数的极坐标数据^[13]。为了便于3D插值以及简化可视化的处理过程，需要将极坐标系下的双偏振原始数据转化为经度、纬度和高度的三维笛卡尔坐标系格点数据。在格点转换中，需要根据不同雷达型号以及不同数据类型的分辨率选择相匹配的格距^[14]。对于双偏振产品，其反射率因子分辨率为250 m，单一径向双偏振反射率因子为920个库。计算双偏振数据三维格点位置时，假设雷达天线位置为极坐标原点，则双偏振原始数据的格点坐标(X_o, Y_o, Z_o)与每一个格点的方位角(a)、仰角(e)和径向距离(l)的计算公式如下：

$ {X_o} = {\rm{cos}}\left( a \right)\cdot{\rm{sin}}\left( e \right)\cdot l $

(1)

$ {Y_o} = {\rm{cos}}\left( a \right)\cdot{\rm{cos}}\left( e \right)\cdot l $

(2)

$ {Z_o} = {\rm{cos}}\left( e \right)\cdot l $

(3)

1.3 网格插值

本文通过插值法将极坐标系中的雷达反射率因子数据内插到笛卡尔坐标系下的三维网格点^[15]。首先在笛卡尔坐标系下建立一个如图 1所示的三维网格模型，选取半径230 km数据，每个格点空间分辨率为0.5 km，故对称后确定该模型水平网格数量为920×920;为充分利用网格，减少无效填充，故只取地面以上10 km数据，垂直格点数为20，每个格点空间分辨率为0.5 km。然后根据笛卡尔坐标系中每一个格点的坐标反算出其在极坐标系中的仰角、方位和距离，方便插值中调用^[16-17]。为使雷达资料插值过程中尽可能地保留雷达资料的原始回波结构特征, 故使用径向和方位上的最邻近与垂直线性插值法。假设雷达馈源位置为(X_t, Y_t, Z_t)，则根据格点坐标(X_g, Y_g, Z_g)计算雷达方位角(a)、仰角(e)与距离(r)的公式如下：

$ a = {\rm{si}}{{\rm{n}}^{ - 1}}\left[ {{\rm{sin}}\left( {\frac{{\rm{\pi }}}{2} - {Y_g}} \right){\rm{sin}}\left( {{X_g} - {X_t}} \right)/{\rm{sin}}\left( {\frac{s}{R}} \right)} \right] $

(4)

式中R为地球半径。弧长s的计算公式如下：

$ s = R{\rm{co}}{{\rm{s}}^{ - 1}}\left[ \begin{array}{l} {\rm{cos}}\left( {\frac{{\rm{\pi }}}{2} - {Y_t}} \right){\rm{cos}}\left( {\frac{{\rm{\pi }}}{2} - {Y_g}} \right) + \\ {\rm{sin}}\left( {\frac{{\rm{\pi }}}{2} - {Y_t}} \right){\rm{sin}}\left( {\frac{{\rm{\pi }}}{2} - {Y_g}} \right){\rm{cos}}\left( {{X_g} - {X_t}} \right) \end{array} \right] $

(5)

三维格点对应仰角e为：

$ e = {\rm{ta}}{{\rm{n}}^{ - 1}}\left\{ {\frac{{{\rm{cos}}\left( {\frac{s}{n}} \right) - \frac{n}{{n + {Z_g} - {Z_t}}}}}{{{\rm{sin}}\left( {\frac{s}{n}} \right)}}} \right\} $

(6)

式中n为大气折射指数，一般取3/4^[18]。三维格点距离雷达馈源的距离r为：

$ r = \frac{{{\rm{sin}}\left( {\frac{s}{n}} \right)\left( {n + {Z_g} - {Z_t}} \right)}}{{{\rm{cos}}\left( e \right)}} $

(7)

如图 2所示，(r, a, e)为格点在球坐标系中的三维坐标，r、a与e分别为距离、方位角与仰角，三维格点(r, a, e)在垂直方向上的相邻仰角分别为仰角e₁与仰角e₂，当相邻仰角小于15°时，三维格点(r, a, e)处的属性值可以用方位角相同处的相邻仰角格点属性值近似表示^[19]。(r, a, e₁)和(r, a, e₂)分别为三维格点(r, a, e)垂直方向上的最邻近上点与最邻近下点，故三维格点(r, a, e)处的反射率值f(r, a, e)可通过f(r, a, e₁)与f(r, a, e₂)插值得到，即：

$ f\left( {r, a, e} \right) = \frac{{{w_{{e_1}}}f\left( {r, a, {e_1}} \right) + {w_{{e_2}}}f\left( {r, a, {e_2}} \right)}}{{{w_{{e_1}}} + {w_{{e_2}}}}} $

(8)

其中w_e1与w_e2分别为插值时格点(r, a, e₁)与格点(r, a, e₂)的权重值：

$ {w_{{e_1}}} = \frac{{{e_2} - e}}{{{e_2} - {e_1}}} $

(9)

$ {w_{{e_2}}} = \frac{{e - {e_1}}}{{{e_2} - {e_1}}} $

(10)

2 计算机三维可视化

目前利用计算机图形学搭建三维网格模型并对模型填色的方法主要有等值面体绘制法与粒子系统仿真法。由于等值面体绘制法适用于数据量巨大，并且对实况显示速度要求极高的情况，故采用等值面体绘制法对每个体扫仰角的格点按数值填色，位于同一仰角的格点生成等值面体，全部等值面体依次堆叠并通过计算机渲染技术生成全彩立方体^[20]，计算机三维可视化流程如图 3所示。

图 3中griddata插值函数调用自Python3.7的scipy库。双偏振雷达回波数据的三维体绘制主要通过Python3.7的Mayavi库实现，Mayavi高度封装的VTK工具包主要负责调用图形模块，并且VTK利用OpenGL来做底层的渲染，可实现网格平滑、锐化过渡、等值线绘制及点线面渲染等三维展示功能^[21]，mlab库函数调用自Python3.7的Mayavi库。绘图采用的色轴通过Mayavi pipeline中的Launch LUT editor自定义得到，三种双偏振产品的绘制透明度Alpha均设置为0.5，对近场非敏感回波透明度Alpha设置为0.8。

3 三维体绘制个例展示

通过个例全面地展示济南雷达双偏振产品的3D显示效果。受西风槽影响，2020年7月5-6日鲁西北东部、鲁中、鲁东南和半岛的部分地区出现雷雨或阵雨天气，局地出现大雨或暴雨并伴随短时强降水、冰雹和大风，潍坊临朐、临沂沂水和日照莒县出现冰雹天气。选取2020年7月6日13:27潍坊临朐境内冰雹为例，济南双偏振雷达观测到东偏南180 km处出现对流风暴(图 4)，1.5°仰角强度为64.5 dBZ，高度6.8 km，3.3°仰角强度为61.5 dBZ，高度12.7 km，0 ℃层高度(4.2 km)以上存在强的反射率因子悬垂，属于典型冰雹特征。根据双偏振产品Z_DR、K_DP与CC的展示效果与侧重点不同，分别对三种产品的二维PPI产品图与三维体绘制效果图进行展示与比较, 其中三维体绘制效果图均可以通过鼠标灵活转动与放大。

3.1 Z_DR产品个例展示

对2020年7月6日13:27济南双偏振雷达Z_DR产品进行3D建模，Z_DR产品图与Z_DR三维体绘制效果图分别如图 5与图 6所示。

图 5a-d分别为0.5°、1.5°、2.4°与3.3°仰角的济南双偏振雷达Z_DR产品，图中蓝色圆圈与冰雹单体强回波区(≥45 dBZ)相对应。可以看出，0.5°仰角强回波区(3.5 km高度)基本对应大的Z_DR值，1.5°仰角(6.8 km高度)以上强回波区基本对应偏小的Z_DR值，同时强回波区一侧存在明显的旁瓣回波，Z_DR表现为正的大值。图 6a-d分别展示了Z_DR产品3D体绘制的前后斜视图、正俯视图与正侧视图。对比图 5与图 6发现，180 km左右存在高度较高的Z_DR大值区，Z_DR大值区既包含Z_DR柱的信息，还包括强回波旁瓣产生的Z_DR大值区。Z_DR柱内含有湿冰粒子或大的液态雨滴，与强上升气流区相对应。

3.2 K_DP产品个例展示

对2020年7月6日13:27济南双偏振雷达K_DP产品进行3D建模，K_DP产品图与K_DP三维体绘制效果图分别如图 7与图 8所示。

图 7a-d分别为0.5°、1.5°、2.4°与3.3°仰角的济南双偏振雷达K_DP产品，图中白色圆圈与冰雹单体各层强回波区(≥45 dBZ)相对应。可以看出，0.5°和1.5°仰角强回波区存在大的K_DP，2.4°仰角以上K_DP明显减小。图 8a-d分别展示了K_DP产品3D体绘制前后斜视图、正俯视图与正侧视图。对比图 7与图 8发现，图 8可以透视地展现出云体中K_DP大值区，环境0 ℃层高度以上有大的K_DP即为K_DP柱。K_DP柱可作为深厚对流上升气流特征的观测度量^[22-25]，强上升气流可将一定浓度的液态雨滴或湿冰带至较高的高度而导致大的K_DP值，三维K_DP体绘制实现了全角度可视，有助于把握强对流回波内部上升气流的强度。

3.3 CC产品个例展示

对2020年7月6日13:27济南双偏振雷达CC产品3D建模，CC产品图与CC三维体绘制效果图分别如图 9与图 10所示。

图 9a-d分别为0.5°、1.5°、2.4°与3.3°仰角的济南双偏振雷达CC产品。图中白色圆圈与冰雹单体各层强回波区(≥45 dBZ)相对应。可以看出，0.5°和1.5°仰角强回波区CC较小，表现为冰粒子和液态雨滴共存。图 10a-d分别展示了CC产品3D体绘制的前后斜视图、正俯视图与正侧视图。

4 小结

通过对济南雷达双偏振产品的解析与坐标变换，利用Python库函数与计算机图形学插值技术，对Z_DR、K_DP与CC双偏振产品数据进行了三维立体(3D)建模重绘，实现了CC、K_DP与Z_DR双偏振产品的三维可视化透视显示。

双偏振产品的三维可视化透视显示对认识风暴强度及强对流风暴双偏振结构特征提供了直观全面的新视角。全角度地展示了强对流风暴发展旺盛阶段在强上升气流作用下常会出现的典型双偏振参量特征如K_DP柱和Z_DR柱。双偏振产品的三维显示可以清晰识别出Z_DR柱和K_DP柱，间接表明风暴内部上升气流强度较强，将一定浓度的液态雨滴或湿冰带至环境0 ℃层高度以上，进而可判别风暴的强度，有助于强对流天气的预警。

双偏振产品的三维可视化透视显示有助于气象工作者对风暴结构的直观理解，有效地将回波信息的垂直分布特征与水平分布特征相结合，快速识别风暴发展强度，提高预警技术水平。