基于SRTM数据地物滤波改进雷达降水估测

赵 文，王澄海，张 强，岳 平，赵 宁，杜莉丽

(1.兰州大学大气科学学院，甘肃省气候资源开发及防灾减灾重点实验室，甘肃 兰州 730000；2.兰州大学地球系统模式研发中心，甘肃 兰州 730000；3.甘肃省气象信息与技术装备保障中心，甘肃 兰州 730020；4.中国气象局兰州干旱气象研究所，甘肃省干旱气候变化与减灾重点实验室，中国气象局干旱气候变化与减灾重点实验室，甘肃 兰州 730020；5.甘肃省天水市气象局，甘肃 天水 741000；6.陕西省气象台，陕西 西安 710014)

引 言

多普勒天气雷达的非气象回波主要有地物回波，超折射回波同波长干扰回波，飞机、船只等回波，海浪回波和由天线辐射特性造成的虚假回波等。在逆温层显著且水汽压随高度迅速减小情况下，地物回波会比通常情况偏多，增加的这部分地物回波被称作超折射回波。对内陆地区而言，地物回波和超折射回波这两项杂波的滤除对于多普勒天气雷达的质量控制以及降水估测[1]具有十分重要的意义。

在WSR-88D雷达中，软件缺省的滤除地物回波方式基于槽口宽度图(notch width map)和杂波过滤旁路图(bypass map)[2]。其中，槽口宽度图给出了一个径向速度范围，对径向速度处于该范围内的回波进行滤除；杂波过滤旁路图说明了不同仰角需要滤除地物回波的范围，其在雷达安装时建立，并在季节变换时更新一次[3]。CINRAD雷达缺省的滤除地物回波方式主要依赖于晴空时的静态杂波图和无限脉冲响应(infinite impulse response, IIR)滤波两种方法[4]。

除了雷达软件自带的地物回波滤除方法外，还有其他一系列方法用于地物回波和超折射回波滤除。如用人工神经网络和模糊逻辑算法(fuzzy logic classifier)判断和识别地物回波和超折射回波[5-9]；引入贝叶斯分类器(Bayes classifier)滤除超折射回波[10]；在地物回波滤波中使用高斯自适应处理模型(Gaussian model adaptive processing, GMAP)[11]，以及在GMAP基础上使用时域高斯自适应处理模型(Gaussian model adaptive processing in time domain,GMAP-TD)，进一步提高对地物回波的滤除水平[12]。

近年来，我国在地物回波的识别和滤波方面也有大量的研究工作。在识别方面，如基于模糊逻辑的分布式超折射地物回波识别研究[13-15]；纯地物回波和地物与降水混合回波两种情况下，利用模糊逻辑法进行地物回波识别[16]；在用反射率因子和差分反射率因子的方差识别地物回波的基础上进一步利用反射率因子垂直廓线订正受地物影响的底层回波[17]；基于同相正交信号(in-phase/quadrature，I/Q)，比较SCI(spectrum clutter identification)算法和杂波消减决策算法(clutter mitigation decision, CMD)对地物回波的识别效果[18]；云雷达的地物回波识别[19]等。而在地物回波滤除方面，有基于偏度方法[20]、综合识别法[21]以及利用地物回波位置、形状和大小基本固定这一特性[22]来进行研究。

以上滤波方式均基于对雷达回波的信号处理。从本质出发，地物回波来源于山脉、丘陵、海(河、湖)岸线、高大建筑物、农作物，甚至高压电源、水塔等形成的干扰回波，其形成与空间地物高程的分布相关。因此，地物高程本身对天气雷达的探测能力和探测质量有着至关重要的影响，国外有研究利用地物高程研究地物回波分布[23-28]。国内利用地物高程的雷达地物回波研究目前多数针对机载雷达[29-31]。对于天气雷达开展的工作不是很多，比如利用地物高程数据评估全国天气雷达网的覆盖能力[32]、天气雷达站选址[33]以及识别和补偿由地物高度造成的雷达波束遮挡[34]。

以往基于地物高程的研究多半着眼于地物回波本身，对地物高程在雷达反演降水方面的应用关注不够。本文从这一角度出发，利用区域(103.5°E—107.5°E，33.0°N —36.5°N)航天飞机雷达地形测绘任务(shuttle radar topography mission, SRTM)资料、2011—2013年甘肃天水CINRAD/CD雷达观测资料，结合陇东南地区国家级自动气象站与区域气象站的小时降水资料进行地物高程在雷达反演降水中的应用研究。首先在一定的假设下确定雷达波束传播的路径，然后结合SRTM资料研究地物回波的分布和滤除，并利用一些较为成熟的方案滤除超折射回波，零度层亮带和非气象斑点回波等，在其基础上，确立适合本区域的Z-I关系参数。

1 资料及强降水过程确定

天水CINRAD/CD雷达(105°38′6″E，34°36′24″N)位于甘肃省天水市秦州区中梁乡西山梁顶，天线馈源海拔为1672.9 m，雷达电磁波波长5.5 cm，水平波束宽度0.95°，垂直波束宽度0.96°，共有9个体扫仰角，探测半径166.5 km，采用RVP7处理器，RVP7对地物回波采用IIR滤波器。该滤波器对强杂波滤除不足，对弱杂波又会过度滤波[35-36]。

选取甘肃天水、平凉2011—2013年10个国家级自动气象站(麦积、甘谷、秦安、武山、天水、清水、张家川、静宁、庄浪、华亭)和240个区域气象观测站的小时降水数据。如果在一次降水过程中，有2个以上国家级自动气象站的最大1 h降水量大于等于8 mm，将其定义为一次强降水，由此提取出6次大范围强降水过程(表1)，只对强降水过程中降水量R≥6 mm的站点数据加以分析。雷达和降水资料及降水分型信息来自天水市气象台。

2 基于SRTM数据计算雷达波束遮挡

为了较准确地确定地物回波分布，需要着眼于两个方面：(1)多普勒天气雷达的波束传播规律；(2)雷达信号覆盖区域内的地物高程。

多普勒天气雷达波束在空气中的传播路径随着大气物理状态的不同而不同。大气中的物理状况错综复杂，通常是时变的非线性系统。为了模拟大气状况从而获知传播路径，需要对大气状态做一定的假设。常用的假设方法有3种：线性折射指数法、指数衰减的折射指数法和射线追踪法[37]。其中射线追踪法精度最高，指数衰减法其次，线性衰减法稍差。若以雷达所在海拔高度为0高度点，线性折射指数法在标准大气折射假定下，3 km高度以下相当精确，而3 km高度以上误差迅速增大。雷达估测降水通常需要选择混合扫描仰角，使得混合扫描高度在1 km左右[32]，而天水雷达站没有探空资料，无法获得折射指数的精确数据，故采用线性折射指数假定。

表1 降水资料信息Tab.1 Information of the precipitation data

2.1 电磁波在大气中的传播

2.1.1 线性折射指数法(等效地球半径法)

假定地球上空折射指数n仅与高度h有关，且折射指数随高度线性减小。公式如下：

(1)

在线性折射指数假设下，可以用等效地球半径法来简化分析。经等效地球半径变换后，雷达波束传播将等效为直线传播，如图1[38]所示。

图1 等效地球半径法模型Fig.1 The effective earth’s radius model

5种地球物理状态(标准大气、临界大气、超折射、零折射和负折射)对应不同的等效地球半径。

(1)标准大气

(2)

图2 计算高度的示意图Fig.2 The diagram of calculating the height

可以推算出：

(3)

(4)

(2)、(4)式联立可以得到沿仰角δ的雷达波束在距离雷达站d的位置波束的海拔高度。由(4)式还可以得到：

(5)

这样，假设距离雷达站某点的海拔高度为H，由(5)式可以得到该点的遮挡仰角，即只有雷达波束仰角高于δ时才能探测到该地上方的大气(这里假设在δ仰角时雷达波束传播到该点之前不会受到阻挡)。

(2)临界大气

对于临界大气折射，可以把地球球面当成一个平面。有R、δ计算公式如下：

(6)

(7)

2.1.2 指数衰减的折射指数

指数衰减的折射指数模型只用到了地面折射指数梯度。假定地球上空折射指数n仅与高度h有关，且折射指数随高度以指数衰减。

n(h)=1+(n0-1)e-ceh

(8)

式中：n0为雷达所在水平面的折射指数；ce为与折射指数梯度相对应的量，计算公式如下：

(9)

2.1.3 射线追踪法

射线追踪法[37]模型需要知道折射指数随高度的分布。假定地球上空折射指数n仅与高度h有关。设雷达波束沿初始仰角θ0传播了R后到达高度为h1的某一点，则

(10)

2.2 数字高程数据

2.2.1 地物高程

地物回波主要是受地形影响而产生的杂波信号，采用雷达探测范围内的地形及高程数据，生成空间高度场，就可以基本确定地物回波的空间分布。

数字高程模型(digital elevation model, DEM)，是在一个规则的格网单元(通常为正方形或经纬度间隔相等的空间单元)上，每一个格网单元都包含一个高程值。常用高程数据[39]如表2所示。SRTM和先进星载热发射和反射辐射仪全球数字高程模型(advanced spaceborne thermal emission and reflection radiometer global digital elevation model, ASTER GDEM)精度最高。SRTM由美国太空总署(NASA)和国防部国家测绘局(NIMA)以及德国与意大利航天机构合作完成联合测量。高精度SRTM数据分两种，即SRTM3和SRTM1，SRTM3分辨率为3″(约90 m)，覆盖全球中纬度地区，SRTM1分辨率为1″，只覆盖美国地区。ASTER由NASA的新一代对地观测卫星Terra的观测制作完成，覆盖范围达到全球陆地表面的99%，在覆盖的所有区域其分辨率都是1″(约30 m)。

在1°×1°经纬网格内，SRTM数据有1201×1201个，ASTER GDTM数据有3601×3601个，雷达波束在径向的分辨率为250 m，在径向的法线方向分辨率更低，SRTM数据精度已经足够，且数据处理量小，对计算资源要求低，故本文采用SRTM数据，数据来源于中国科学院计算机网络信息中心国际科学数据镜像网站(http://www.gscloud.cn)。

对中国地区SRTM数据高程误差的研究[40]表明，SRTM高程数据的平均误差为-0.35 m，90%误差为7.4 m，达到绝对高程误差小于16 m的设计要求。对中国西北部高山草原地区SRTM数据的研究也表明该数据在西北高山草原地区同样具有极高的可信度[41]。

表2 常用高程数据Tab.2 The common DEM data

经过计算，天水雷达覆盖区域基本处于103.5°E—107.5°E、33°N—36.5°N网格范围内，故对SRTM数据做上述范围内的拼合，得到4801×4201组网格高程数据。

2.2.2 高程插值

雷达有360个径向束，为了获得地物遮挡情况，需要求得360个径向束所对应的海拔高度，为了确定其海拔高度，首先需要知道每个点的具体经纬度坐标。

用λ和φ分别表示经度和纬度。下面求沿雷达站方位角θ(指从某点的指北方向线起，依顺时针方向到目标方向线之间的水平夹角)前进距离d到达点的经纬度，同样用R来表示斜距，设雷达站所在位置为A(λ1,φ1)，所求点为B(λ2,φ2)。公式如下：

δ=d/R

(11)

(12)

(13)

△λ=δ·sinθ/q

(14)

φ2=φ1·cosθ

(15)

λ2=λ1+△λ

(16)

根据公式(11)—(16)，就可以求出从雷达站出发的波束沿特定方位角沿地球面行进特定距离后到达的位置，由此可以得出360个方位角上每行进250 m到达点的经纬度坐标，从250 m直到166.5 km，每条径向上有666个点，共360×666个点的经纬度。考虑到雷达波束的方位角并不恒定，为了避免每次都计算一次方位角，并且方便后面计算波束遮挡率，这里加大分辨率，以0.1°为间隔，从0°一直到359.9°共3600个方位角每行进250 m所到达的位置，共3600×666个点，然后将SRTM数据插值到这3600×666个点上以获取地物高程分布。

高程空间插值有很多方法。比如ArcGIS软件中常常用到反距离加权法(inverse distance to a power, IDW )、克里金法(Kriging)、最近距离法(nearest neighbor, NN)、样条函数法(spline)等。本文选取双三次Bezier曲面和B样条曲面这两种样条函数进行插值试验，按照徐良等[42]提出的检验方法，选择SRTM数据无缺失的16个点，分别用两种方法进行插值，并与真实值进行对比。双三次Bezier曲面插值16个点的平均相对误差为1.20%，而B样条曲面插值16个点的平均相对误差仅为0.71%，表明B样条曲面效果更好。故采用B样条曲面进行高程空间插值。

2.3 地物回波和波束遮挡率

2.3.1 地物回波分布

由B样条曲面在3600个径向上插值后，就可以计算地物遮挡情况，每一条径向，从雷达站开始，由近及远， 每250 m取一点，根据该点高程和前一点遮挡角来计算该点遮挡角，由于缺乏探空资料，无法确定超折射回波参数，这里仅计算标准大气折射和临界大气折射2种情况。天水雷达站周边海拔分布如图3所示。

图3 天水雷达站周边海拔(单位：m)Fig.3 The elevation of the region around Tianshui radar station (Unit: m)

结合2.1节，分别计算标准大气折射和临界大气折射情况下天水雷达站周边区域的雷达波束遮挡角，如图4所示。标准大气折射时整个雷达覆盖范围内最大遮挡角为1.21°，临界大气折射时最大遮挡角为1.56°。

2.3.2 波束遮挡率

实际中，雷达波束不是一条直线，而是一个椭圆锥体，如图5所示。天水雷达垂直波束宽度为0.96°，水平波束宽度为0.94°，雷达波束就在这两个波束宽度所定义的椭圆椎体内。故应该考虑不同方位雷达波束遮挡率，其定义为有效照射体积(波束)内因地形等障碍物遮挡而损耗的功率比[34]。

图4 天水雷达站周边遮挡角(a)标准大气折射，(b)临界大气折射Fig.4 The blocking angle of the region around Tianshui radar station(a) standard refraction, (b) critical refraction

图5 雷达波束遮挡示意图 (a)径向束方向 ，(b)径向束法线方向Fig.5 The diagram of blockage of radar beam(a) radial direction， (b) normal direction

计算波束遮挡率的方法是以0.1°为基本计算单元，以雷达波束仰角为中心轴，向上、向下、向左、向右各扩展n=15条，一共是1.5°×1.5°大小的区域(31条×31条)。对每个单元分别检验其遮挡情况，则遮挡率O的计算公式如下：

(17)

其中：

(18)

(19)

式中：θ1为水平波束分辨率；φ1为垂直波束分辨率；m为第n个方位上的遮挡标号，取值从-16到15；B(n)为第n个方位的遮挡程度；W(|n|)为第n个方位在总波束功率中所占的功率比。

由此得到标准大气折射下雷达的波束遮挡率，发现在2.4°仰角，已经完全没有遮挡，0.5°和1.5°仰角的遮挡率如图6所示。根据计算结果，剔除掉波束遮挡率大于0的资料。

图6 标准大气折射下雷达波束遮挡率 (单位：%)(a)0.5°仰角，(b)1.5°仰角Fig.6 The blockage rate of radar beam for standard atmospheric refraction (Unit: %)(a) 0.5°elevation angle, (b) 1.5°elevation angle

3 Z-I关系及效果评估

除了地物杂波外，雷达主要的杂波信号还有超折射回波、零度层亮带和非气象斑点回波等，这些信号也会显著影响雷达的探测质量。对于这些杂波信号，已经有很多成熟的滤波方案，如刘黎平等[13]发展了由KESSINGER等[8-9]提出的基于模糊逻辑的超折射地物回波滤波方法；江源等[14]统计了该方案中用于检测地物和超折射回波的几种识别因子在C波段、S波段的正常降水回波以及地物和超折射回波中的分布差异；李丰等[15]进一步利用上述统计结果，针对C波段多普勒天气雷达对该方案中的参数进行了优化。庄薇等[43]参考ZHANG等[44]和张乐坚等[45]的研究，提出了高原地带零度层亮带的识别和订正方法。肖艳姣[34]探讨了非气象斑点回波的滤波方法。本研究除了利用SRTM资料滤除地物杂波外，还基于以上这些方案对这3种杂波信号进行了滤除。

3.1 数据预处理及算法

3.1.1 时次及仰角选取

利用已经滤除上述杂波的雷达反射率因子对Z-I关系的本地化进行研究。通常认为雷达反射率因子和降水强度之间存在一定的定量关系，即Z-I关系：

Z=AIb

(20)

只要获得了上式中的A，b两个参数，就可以由反射率因子对降水进行估测。不同雷达对应的这两个参数有所差异，通常的Z-I关系主要基于S波段雷达给出，对于C波段雷达的适用性不是特别好，所以上述关系的本地化研究对提高雷达预测水平有很重要的意义。不同类型降水的物理机制和大气条件有明显差异，因此分别对对流性降水、混合性降水和非对流性降水3种类型降水得出其各自本地化的Z-I关系。由于不同降水系统的降水效率等原因，猜测雷达回波对降水有一定超前，为此，对小时降水和1 h内雷达反射率因子的滑动平均做相关分析。而雷达反射率因子和小时降水量不是简单的线性关系，由反射率因子dbZ=10lgZ可以推出：dbZ=10lgA+10blgI，故lgZ和lgI成线性关系。这里研究1 h内lgZ的滑动平均值与小时lgI之间的相关关系。

天水雷达1 h体扫次数(11或12次)不固定。故计算6次降水过程中lgZ的滑动平均值与降水时次lgI的相关系数(图7)。可以看出，在时间超前上，对流性降水情况较复杂，这里取均值5 min；非对流性降水，反射率因子超前降水5 min；混合性降水，反射率因子和降水同步。采用1 km混合体扫反射率进行雷达降水反演，在距离雷达站0～16 km采用3.4°仰角，16～23 km采用2.4°仰角，23 km以上采用1.5°仰角。对于受遮挡部分，用不受遮挡(即遮挡率为0)的最低仰角数据代替。在每个时次均用双线性插值将反射率因子插值到各站点上空。

图7 6次降水过程lg Z滑动平均值与小时lg I的相关系数(数字12代表参与滑动平均反射率因子均在降水时次内；11则代表去掉最后反射率因子并在最前面相应加入一组值进行滑动平均，依次类推)Fig.7 Correlation coefficient between moving average value of lg Z and 1 hour lg Iduring 6 precipitation processes(The number 12 represents that the reflectivity factors participating in the moving average are with in the precipitation hour; and number 11 represents that the last reflectivity factor is removed and the one before the first values are added to the front for moving average, and so on)

3.1.2 最优化拟合法

雷达反演降水中，最常用的方法是最优化拟合法，即事先假定一个Z-I关系，将所有Z值换算成雨强I，对I进行时间积分获得小时降水量的雷达估算值Hi，选定判别函数CTF(Hi与气象站实测降水量Gi的函数)，将雷达数据和气象站数据代入，如CTF过大，就不断调整Z-I关系中的参数A和b，直到CTF最小。

目前比较理想的一个判别函数形式为

CTF=min{∑i[(Hi-Gi)2+|Hi-Gi|]}

(21)

采用此判别函数，对6次3类降水过程分类型进行最优化拟合，在降雨过程中A的变化范围通常为16～1200，b值为1～2.87，考虑到该结论基于S波段雷达得出，未必切合C波段雷达，本文取A值范围为1～1200，步长为1；b值范围为0.5～3，步长为0.1。针对每一种类型降水，在上述范围内，确定CTF最小时A和b值，所得公式为最优化拟合公式。

为了检验拟合效果，与实践中用两种经验得出的效果较好的Z-I关系进行比较。对暴雨来讲，常用的一种Z-I关系为

Z=486I1.37

(22)

我国新一代天气雷达默认的Z-I关系式为

Z=300I1.4

(23)

通过比率(RATIO)、平均相对误差(ARE)、均方根差(RMSE)和相关系数(COR)4个参数[46]比较这两个经验公式和最优化拟合法结果，其中，RATIO越接近1，ARE和RMSE越小，COR越大，拟合效果越好。具体公式如下：

(24)

(25)

(26)

(27)

(28)

定义一个新的参数反演效果系数λ21来表征其他两个经验公式相对于最优化拟合公式的反演效果，反演效果系数越小，证明其他公式的反演效果相对于最优化拟合公式来说越差，λ21公式如下：

(29)

式中：τ1、ARE1、RMSE1和COR1为最优化拟合公式对应的统计量。τ2、ARE2、RMSE2和COR2为其他两种经验Z-I关系对应的统计量。

3.2 拟合结果及效果检验

表3列出最优化拟合结果。可以看出，天水雷达的Z-I关系中，A值相对S波段雷达的结果偏小，b则偏大。为进一步验证本研究中拟合所得结果的合理性，对3种类型降水的最优化拟合公式与其他两种经验公式的结果进行比较，结果如图8所示。可以看出，所有指标最优化拟合公式均优于两种经验公式。两种经验公式相对于最优化拟合公式的反演效果系数如图9所示，同样表明最优化拟合公式相对于经验公式有明显改进。

表3 最优化拟合结果Tab.3 The results of optimal fitting

天水雷达Z-I关系中的A参数显著偏小，这似乎不能仅仅用雷达波段以及地理和气象条件的不同来解释，为了探究A参数小的原因，利用临近的庆阳西峰新一代天气雷达与天水雷达进行了比较。庆阳雷达与天水雷达的型号(CINRAD/CD)与探测半径(166.5 km)均一致；二者的距离库长有一定差异，庆阳雷达为300 m，天水雷达则为250 m，二者之间不存在严重的波束遮挡，因此有很好的可比性。两部雷达探测范围如图10所示，可以看出，二者探测范围的重合区域大致在106°E—107°E、35°N—36°N范围内，且这一范围和两部雷达的距离差距不大，因此选择该范围为比较区域。对前面给出的3种降水类型，各自选取一个时次，使得该重合区域在该时次内有显著降水。由于0.5°仰角有明显的遮挡，所以选择1.5°仰角进行比较，结果如图11所示。可以看出，无论是哪一种类型降水，天水雷达回波的反射率因子均小于同一时次庆阳雷达，因此天水新一代天气雷达可能存在系统性的回波强度偏低。

图8 不同参数Z-I关系比较(a)RATIO ，(b)ARE，(c)RMSE，(d)CORFig.8 The comparation of Z-I function with different parameters(a) RATIO, (b) ARE, (c) RMSE, (d) COR

图9 两种经验公式相对于最优化拟合公式的反演效果系数Fig.9 The inversion-effect coefficient of two empirical formulas with respect to the optimal fitting formula

图10 两部雷达的探测范围Fig.10 The detection ranges of two radars

图11 不同类型降水时天水(a、c、e)、西峰(b、d、f)两部雷达在重叠区的反射率因子比较(a、b)对流性降水，(c、d)混合性降水，(e、f)非对流性降水Fig.11 Comparison of reflectivity factors in overlapping area of Tianshui (a, c, e) and Xifeng (b, d, f) radars under different types of precipitation(a, b) convective precipitation, (c, d) mixed precipitation, (e, f) non-convective precipitation

4 结论与展望

(1)对反射率因子和降水的同步性研究表明，混合性降水中反射率因子和降水同步；非对流性降水中反射率因子超前降水一个体扫；对流性降水情况较为复杂，不可一概而论。

(2)天水雷达的最优化Z-I关系和其他两种经验公式不同，A值偏小而b值偏大，而国内常用的Z-I经验公式对天水雷达的估测水平不高。通过与邻近的庆阳西峰新一代天气雷达重合范围内的反射率因子进行比较，发现天水雷达存在系统性的回波强度偏低，其原因有待进一步研究。

由于天水雷达站没有探空数据，故在研究雷达波束传播时无法采用相对更准确的射线追踪法和指数模型。下一步应该着手在天水雷达站进行探空工作，以获得探空数据从而更精确地模拟大气物理状态以及雷达传播情况。