基于数字高程模型和反射率因子垂直廓线的天气雷达波束遮挡订正算法

周淑玥, 李 静, 鲜 林, 杜 捷, 王海江

(成都信息工程大学电子工程学院,四川 成都 610225 )

0 引言

由于X 波段双偏振天气雷达的波长相对于S 波段和C 波段更短,所以在穿过强降雨区域和遇到山体阻挡时更容易产生衰减和波束阻挡问题[1]。 受建设条件与选址等方面的影响,很多已部署天气雷达的周围通常存在较多山体结构,雷达波束在传播路径中受到遮挡,会对雷达的探测数据产生极大的影响。 尤其是在雷达的低层仰角,遮挡对于雷达数据的影响更加明显,不利于低空天气观测。 尽管中国天气雷达部署区域正在进一步扩大,雷达的大部分遮挡情况可通过雷达组网的方式得到改善,但受到多种条件的限制,部分区域仍存在波束遮挡。

针对天气雷达的波束遮挡问题,国内外学者提出了遮挡区域的数据订正和填充方法。 Harrold 等[2]提出一种基于雷达反射率因子分布特征计算波束遮挡率从而补偿遮挡区的方案,该方案在波束遮挡率不高的情况下具有一定效果。 Vulpiani 等[3]利用双线偏振雷达在意大利山区的复杂地形条件下,实现了差分传播相位常数对雷达波束遮挡的补偿订正,并将订正结果应用于雷达降雨估计中,结果表明差分传播相位常数能很好地补偿波束被遮挡区域的回波,其降水估测结果优于反射率因子的估计结果。 Ryzhkov 等[4]利用衰减系数实现了在复杂地形条件下的波束遮挡订正和降水量估算,结果表明该类算法能够在雷达系统标定不佳、雨衰严重和雷达波束部分遮挡严重的情况下提供准确而可靠的雨量估计。 张亚萍等[5]针对部分遮挡区域提出基于数字高程模型的平均值距离库填充法,根据相邻仰角层的反射率因子,实现了动态订正回波。肖艳娇等[6]分析了相邻仰角层的反射率因子差与波束阻挡率的关系,证实在标准大气折射条件下使用DEM(digital elevation model)数据计算波束阻挡的合理性,并利用方位伸展对部分遮挡区域进行回波订正。勾亚彬等[7]提出在对雷达组网拼图前,可先将雷达部分遮挡区域的回波去除,以减小不连续效应。

本文基于SRTM(shuttle radar topography mission)的DEM 数据对北京房山区站点的X 双偏振雷达周围75 km进行了波束遮挡率计算,同时利用反射率因子平均垂直廓线对雷达低层中的遮挡区域回波进行修正。 最后采用定量分析对平均垂直廓线法填补遮挡区域的效果和其相关影响因素进行分析。

1 资料

采用的数据包括数字高程地形数据(DEM)、X 波段双偏振雷达体扫数据。 数字高程资料来自由美国航空航天局NASA(National Aeronautics and Space Administration)和国防部国家测绘局NIMA(National Imagery and Mapping Agency)联合测量的中国境内航天飞机雷达地形测绘使命SRTM 数据,选取北京区域分辨率为90 m×90 m 的数据,绘制如图1 的地形图。 图1 中三角形符号代表房山雷达站点位置,由地形图可知该站点的主要遮挡物是西北方向的高大山脉,而在其他方位则地势较为平坦。

图1 北京地区DEM 地形图

实验数据来自北京房山站点的X 波段双偏振雷达探测到的2020 年8 月一天的天气数据,验证数据选取北京房山、昌平、顺义3 个站点X 波段双偏振雷达数据,3 个站点的雷达详细参数如表1 所示。

表1 雷达站点信息

3 部X 波段双偏振天气雷达,观测仰角均为9 层(0.5°,1.5°,2.4°,3.4°,4.3°,6.0°,9.9°,14.6°,19.5°)。

2 基于数字高程模型计算波束遮挡率

2.1 遮挡率计算方法

本文采用数字高程模型模拟雷达探测时的波束遮挡率,在计算时需同时考虑大气折射和雷达波束与周围地形的相互作用。

针对大气折射的影响,在计算雷达每个距离库的高度时,使用引入地球等效半径的大气折射测高公式[8]。

式中,地球曲率半径a为6370 km,ae为地球等效半径,R为径向距离,hs为雷达天线的海拔高度,φ为波束仰角。

雷达波束被地物遮挡的程度可用BBF(beam block fraction) 来描述,它被定义为在波束的有效照射体积内被地表障碍物遮挡而损耗(Pbl)的功率与波束内总功率(Porig)之比[9]。 换句话说,BBF 表示因遮挡引起的功率部分损失比例。 若BBF=1(或100%),则表示理论上电磁功率被完全遮挡。

利用DEM 数据求得雷达站对应的遮蔽角后,通过对地物遮挡区投影到水平和垂直方向的波束内电磁功率分别进行积分来计算波束遮挡率。 由于受地物遮挡的电磁波功率损耗和波束内电磁波总功率均与雷达天线功率方向函数相关,因此需要精确地描述天线的辐射模式[10]。 实际上,气象雷达的天线都是具有高度方向性的定向天线,能将大部分能量集中在一个很窄的波束范围内并朝指定方向发射出去。 考虑到天线辐射能量在波束内分布不均,雷达能流密度在波束中轴线方向上最为集中,且会随轴线方向沿其左右两边迅速下降直到几乎为零。 因此可引入基于高斯分布的天线功率方向函数来近似描述波束内非均匀电磁能量分布模式[11]。在天线功率方向函数中,天气雷达的波束宽度通常采用两个半功率点的夹角表示,其表达式如下:

式中,θ1、φ1分别为雷达波束沿径向在水平和垂直方向的宽度,θ、φ分别为波束内的某一点距离波束中心的方位角和仰角。

根据经典气象雷达方程,在无地物遮挡的情况下,径向距离R处的雷达回波功率Porig为

式中α0为波束中心仰角。 由地物遮挡产生的功率损失与总回波功率之比即为BBF,其表达式:

2.2 遮挡率计算结果与分析

根据DEM 模型模拟的雷达分别在0.5°、1.5°、2.4°仰角的波束遮挡率BBF 计算结果见图2(a)、(c)、(e),对应的雷达反射率PPI(plan position indicator)见图2(b)、(d)、(f)。 通过对比BBF 与雷达反射率PPI 图像,可看出0.5°、1.5°仰角层的BBF 图中深色区域(高遮挡率区域)与雷达反射率图像中的缺测值区域范围基本一致,随着雷达仰角的抬高,因地物而产生的波束遮挡逐渐减小,从2.4°仰角开始完全遮挡区域消失。 0.5°仰角层的径向角约120° ～140°区域与其周围和上层反射率相比明显数值较低,与其BBF 图像中东南角的波束遮挡区域相对应。 由此可见,根据DEM 模型模拟的波束遮挡率基本能代表实际的雷达波束遮挡情况,可根据其对遮挡区域进行反射率因子订正。

图2 房山站点低层仰角遮挡率结果与雷达反射率对比

3 平均垂直廓线法订正

为减少X 波段雷达的衰减特性对于实验结果的影响,在对遮挡区域进行订正之前先对X 波段雷达的各层反射率数据进行衰减订正。 衰减订正算法采用Bringi 提出的AH-KDP算法,衰减订正公式如下:

式中,α为衰减系数,AH为反射率因子衰减率,KDP为差分传播相移率。 衰减订正结果如图3 所示。

图3 房山站点雷达衰减订正后反射率因子PPI

遮挡区域填补前后的对比效果如图4 所示,由图4 可见0.5°与1.5°仰角的因遮挡造成的缺测区域经过处理后明显都得到了订正,且与周围的非遮挡区域的反射率因子无论是在径向还是切向方向都具有较好的连续性。

图4 遮挡区域回波填补前后效果对比

3.1 平均垂直廓线计算

反射率因子垂直廓线生成方法主要有参数法、平均法和识别法3 种,其中平均法计算简单,实用性强,被广泛应用于雷达降水订正[12]。 本文采用平均法生成反射率因子垂直廓线MVPR(mean vertical profile of reflectivity)。

对一个仰角为θ的PPI 划分高度层如图5 所示,Δh为层间高度,方位角夹角1°,Δh高度对应的径向上增量为Δr的单元称为bΔh,它与沿着雷达径向夹角为1°范围的1 km增量b不一定相等,而与θ,h有关。

图5 MVPR 计算示意图

计算MVPR 的方法如下:

式中θ为雷达扫描仰角,α为方位角,r(h+Δh)和r(h)表示竖直高度的距离,N代表雷达体扫个数。 式中分子是高度的反射率因子之和,分母是参考高度的径向增量和。

针对MVPR 的生成过程主要有选取的生成区域和分层厚度两个影响因素:选取的MVPR 生成区域如果距离雷达较近,可获取更多底层数据,但容易受到地物和雷达旁瓣干扰,距离雷达太远,波束垂直间隔过大,雷达低层仰角数据量较少。 综合考虑下,选择距离雷达1 ～50 km区间的观测数据,使得到的MVPR 相对光滑且有足够的垂直分辨率。

除了距离,选取生成MVPR 的区域时,应该选取波束遮挡率在10%以下的非遮挡区域,才能够为后续的遮挡区域反射率因子填充生成可靠的垂直廓线模型,选择径向上210° ～240°作为生成MVPR 区域。 分层厚度的选择上,分层厚度若选择过大,生成的反射率因子垂直廓线的形状细节特征容易丢失;若分层厚度选择过小,廓线细节体现充分但计算量将大幅增加且廓线形状不稳定,因此分层厚度选择为250 m。

根据上述MVPR 的区域选取和计算方法,对于房山地区的一次降水体扫数据进行廓线计算,其结果如图6 所示。 通过廓线可以看出该体扫数据的较强回波集中在6 km高度以下,且低空的反射率主要分布在25 ～30 dBZ,与图3 中各层仰角的PPI 图像回波情况基本一致。

图6 MVPR 廓线图

3.2 根据廓线填补遮挡区域

仰角上每个距离库的高度H可根据考虑大气折射和地球曲率的测高公式算得,然后根据垂直分层距离确定波束遮挡点位于MVPR 上的区间位置,再通过插值得到MVPR 对应高度上的反射率因子值Z1(单位dBZ)和该距离库上层无遮挡区域仰角在MVPR 上的反射率因子Z2。Zm为该距离库对应的上层无遮挡仰角的实际反射率因子,则通过计算,可得到完全遮挡区的反射率因子值Z:

3.3 算法效果分析

为分析算法的填补效果,将一个遮挡率为10%以下的无遮挡区的去除回波,并将其假设为完全遮挡区。然后利用MVPR 和该区域上层仰角的反射率因子填补该区域,将该区域填补值与实际反射率因子值进行对比,分析两者差异。

采用房山、昌平、顺义3 个雷达站点的0.5°仰角填补后反射率与实际观测反射率数据进行相关性分析。 表2 列出3 个站点不同时间在20°和170°两个方位角的填补数据与对应位置“真值”的相关系数,结果显示相关系数集中分布在0.9 左右,证明填补后的值与实际反射率因子一致性较高。 表2 中方位角20°的填补值和真值的相关性大部分都比170°方位角的高,这是由于各方位角的反射率因子垂直反射率分布情况与生成的平均反射率因子垂直廓线差异大小不同,因此在各方位角表现的填补效果也有差异。

表2 填补值与“真值”的相关系数

4 结束语

经过对比DEM 模型模拟的波束遮挡率与实际反射率因子PPI 图片可知,波束遮挡率较大处与反射率因子缺测处具有较好的对应关系,因此使用DEM 模型模拟的雷达波束在传递过程中的遮挡率BBF 能够很好的代表雷达实际的遮挡情况。 在利用经过衰减订正的反射率因子数据生成平均反射率因子垂直廓线后,根据波束遮挡率使用MVPR 进行遮挡填充。 经过对该方法的定量分析后,证明该方法在波束遮挡区域填充的反射率因子与“真值”一致性较好,具有很好的适用性。