1 引言

中国是自然灾害频发的国家, 其中气象灾害如洪涝、冰雹、台风等是给中国带来最大损害的自然灾害之一。天气雷达作为一种主动微波遥感技术, 已成为强对流天气的重要观测和预警手段之一。高时空分辨率的降水观测对强对流降水的精确监测、数值模式同化预报至关重要。地基雷达如中国CINRAD SA/SB雷达能提供高时间(约6 min)和空间(约1 km×1 km)分辨率的测量, 但其易受波束阻挡、地物杂波等的影响, 且在远距离受波束展宽的平均作用而降低分辨率(俞小鼎等, 2006)。星载降水雷达可提供近全球范围的三维降水观测, 可探测到地基雷达不易到达的地方如海洋、山区等, 且不受波束阻挡的影响, 但其时空分辨率较低(Kozu, et al, 2001; Hou, et al, 2014)。因此, 如何通过改进降水观测手段或升级观测设备、以及通过观测数据后期处理如插值或超分辨率重建, 获得高分辨率中小尺度强对流降水数据依然是值得研究的课题。

随着越来越多的降水探测工具的出现, 将多源数据融合或联合测量在以后的降水数据处理和应用中将愈加频繁。如2014年2月发射成功的全球降水观测计划(Global Precipitation Mission, GPM)主卫星, GPM计划由主卫星和8颗星座卫星联合观测实现每隔3 h覆盖全球一次(Hou, et al, 2014), 以及中国计划在风云三号02批卫星观测星座中的测雨星上搭载双频测雨雷达(吴琼等, 2013), 这为以后的多源多分辨率降水观测数据的融合和联合处理及应用提供了更多的机会和挑战。在不同平台观测数据联合处理或对比中, 为匹配不同时空分辨率的数据, 插值或重采样处理不可避免(寇蕾蕾等, 2016; Kou, et al, 2018; Wang, et al, 2009)。此外, 在多部地基雷达组网拼图、雨量计校准雷达定量测量降水、雷达资料同化研究等气象应用中, 可能都会用到天气雷达数据精确有效的插值或空间匹配处理(肖艳姣等, 2006; Bocchiola, 2007)。

基于不同的研究目的和应用需求, 有不同的插值方法, 如最近邻点、线性插值、Cressman权重插值、傅里叶插值等(黄云仙等, 2008; 焦鹏程等, 2016; 肖艳姣等, 2006)。最近邻点、线性、Cressman权重插值均属于空间域插值法, 这些插值方法共同的优点是算法简单运算效率高, 但这些滤波函数均为低通滤波器, 对数据中的高频信息有滤除和抑制效应。傅里叶插值通过在频率域高频补零实现时域插值, 改善图像的视在分辨率, 但它将整个序列当做周期信号, 使得待插值点与整个序列包络有关, 这可能造成一些失真的孤立强降水信号。为提高天气雷达时间分辨率, Ruzanski等(2015)提出基于核方法的傅里叶插值法, 它通过加窗来筛选有效的待插值输入数据, 可改善常规傅里叶插值法的周期输入数据局限性。另一些天气雷达图像插值或降尺度方法是从降水数据结构的几何及统计特性分析出发, 目的是为有效保持数据本身的时空变化特点, 如随机分形插值法、基于规整化的高分辨率重建法等(Deidda, 2000; Ebtehaj et al, 2012, 2013)。它们通常通过建立图像的先验模型, 将插值问题转化为约束优化问题, 如Foufoula-Georgiou等(2014)提出的一种基于变分规整化的星载降水雷达数据高分辨率重建方法。

为精确刻画自然图像的高阶边缘统计特性及系数间高阶依赖关系, 人们提出小波域高斯尺度混合(Gaussian Scale Mixtures, GSM)模型(Portilla, et al, 2003; Wainwright, et al, 2000), 通过一个简单局部依赖关系模型和一个控制参数的随机变量来联合建模。由于降水数据尤其强对流降水常呈现出明显的簇拥和聚集性, 出现一些降水单体, 而这些孤立的强降水单体在直方图分布中表现出比常规高斯分布更厚的尾部, 即重尾性。为更好地抓住降水数据结构这些重要的高阶统计特征, Ebtehaj等(2011a)将描绘自然图像的GSM模型引入到天气雷达数据建模中, 并已将此模型应用到了地基和星载雷达降水融合中(Ebtehaj, et al, 2011b)。

降水回波的小尺度变化对强对流天气监测和预报至关重要。本文基于雷达降水数据空间域统计特征, 采用小波域GSM模型来构造雷达回波强度数据先验模型, 并将其应用到天气雷达图像插值中, 以得到能有效反映强降水回波及小尺度变化细节的高分辨率雷达回波强度图像。基于个例数据分析, 全面讨论了天气雷达回波强度数据小波域的非高斯重尾分布特性、邻域间相关性及尺度间的依赖性, 并将这些特性与GSM模型进行了匹配; 进而基于GSM先验模型和尺度间的关系模型, 利用贝叶斯最小均方差估计出小尺度高频子带小波系数; 最后, 进行小波逆变换, 实现天气雷达图像高分辨率插值。该方法通过学习降水回波小波域尺度内和尺度间的统计关系, 建立相应模型估计更小尺度高频系数, 其充分考虑了数据本身的特点, 能较好重建雷达回波强度图像高分辨率细节且不会引入虚假高频信息, 可进一步改善高分辨率中小尺度降水数据结构精细探测及中尺度预报中的雷达数据同化等应用。

2 雷达回波强度数据小波域GSM建模 2.1 雷达回波强度数据小波域统计特性

对于天气雷达回波强度数据, 强降水对流单体常常群簇出现在一片弱降水中, 从而表现出高度聚集性及稀疏相关性。从图像的特征来看, 即表现出无降水与有降水、强降水与弱降水间的方向边缘性。小波变换通过对信号进行不同尺度不同方向的带通滤波, 从而有效获取信号的局部不连续性或起伏变化等细节特征。如雷达降水数据的群聚性在小波频率域, 常表现出明显的小波系数分布的尖峰和重尾特性(Ebtehaj et al, 2010, 2011a; Harris, et al, 2001)。

图 1是2008年5月27日南京地区一次雷暴过程个例(以下简称个例1), 其中图 1a是原始体扫数据经3次线性插值重采样到笛卡尔坐标系后的3 km等高平面显示(CAPPI), 网格化分辨率为1 km×1 km, 图 1b—d分别是图 1a中数据小波分解后的水平向、垂直向和对角向子带的小波系数, 分别对应图 1a中横向细节、垂直细节和对角细节信息, 图 1e—g分别对应图 1b—d中小波系数的概率分布统计。从图 1e—g可看出, 不管是水平、垂直或对角方向子带, 雷达图像的小波系数的概率分布均主要呈现为大量小值和少量大值的构成形式, 数值方差较大, 即非高斯“重尾”特性。这种重尾性反映了降水回波中的小尺度变化及一些强回波奇异值, 这也是对雷达回波强度数据建模时需要注意的一个重要特点。

需要注意的是, 在进行图 1的小波分解时, 采用一种非抽样平稳离散小波变换(Nason, et al, 1995), 它在分解时去除“下采样操作”, 消除了“频率混跌项”, 从而不会因为重构处理带来虚假回波。此外, 在进行小波系数的统计分析时, 去除了无用背景零点(无降水回波)的影响, 而保留了因分解过程中与滤波器卷积作用而生成的有用零点(Kuma, et al, 1994)。

图 2是另一个降水个例(南京地区, 2011年7月25日, 以下简称个例2)及其小波系数的统计概率分布, 其中图 2a是雷达回波, 图 2b中实线是回波数据小波分解后水平子带的概率统计, 图 2b中虚线显示了常规的高斯分布。与图 1中的小波系数概率分布类似, 图 2b中的概率分布也表现出明显的重尾性, 即离峰值较远的少量大值也有相当的出现概率。但个例2中回波面积相对较少, 而强回波占大部分, 图 2b中概率分布的方差略小于图 1e中方差。

除了小波系数边缘概率分布的重尾性, 由于降水数据的多尺度特征和群聚特性以及小波分解尺度间的相似性, 雷达回波强度在小波域还表现出“尺度内相关性”和“尺度间依赖性”(Ebtehaj et al, 2011a, 2011b)。个例1水平、垂直和对角子带内5×5邻域小波系数的协方差矩阵如图 3所示(个例2图像效果与此类似, 图略), 反映了降水回波强度数据空间分布邻域间的相关性。图中越白的像素表示相关性越高(向1靠近), 越暗的像素表示相关性越低(向0靠近), 非对角线上的非零元素值反映了系数局部的相关性。从图 3可看出, 虽然进行了近似正交的小波分解, 但较近的邻域范围内数据相关性依然较大, 较远地方的相关性较弱。这说明经过小波分解后的降水回波只与周边邻域范围的像素有关, 在进行插值估计更高分辨率回波系数时, 仅考虑待插值点周边局部降水回波数据协方差矩阵即可。

图 4是图 1中雷达反射率因子数据进行二次小波分解后, 尺度1(对应较精细尺度)和尺度2(对应较粗糙尺度)小波系数的二维联合直方图, 它反映了降水数据小波系数尺度间的依赖性。由图 4可以看出, 水平向和垂直向子带系数表现为明显倾斜的非零协方差矩阵元素, 说明降水回波在水平和垂直变化上有高阶相关性, 尺度间的依赖性较强, 而对角向尺度间数据近似不相关, 依赖性较弱。

2.2 雷达降水回波强度数据小波域GSM模型

为了精确刻画雷达回波强度数据的非高斯重尾特性及尺度内邻域系数间相关性, 采用GSM模型来构造降水回波强度数据小波系数的先验模型。

假设天气雷达图像分解后的子带内小波系数每一局部邻域系数X符合GSM模型, 则向量X为高斯混合尺度, 即它可表示为0均值高斯向量U和独立正尺度随机因子的乘积(Wainwright, et al, 2000)

(1)

式中, 表示在分布上相等, 因子z为权系数, 则向量X的概率密度由U的协方差矩阵C_U和随机数z的密度p_z(z)决定

(2)

式中, N_n为向量X和U的维数, 即所选择的邻域的大小。由式(1)可得

(3)

式中, E(z)表示对变量z取总体平均。式(3)说明向量X的协方差C_X主要由C_U和变量z的平均得到。一般情况下, 设置随机数z的平均值为1, 则C_X与C_U一致。

由于邻域相互重叠, 每个系数同属于多个领域, 这种局部模型暗含全局特征。由式(2)可看出, GSM模型是一种广义的混合高斯模型, 它相当于对高斯分布进行加权混合, 加权系数为正随机数z。随机数z的方差将控制图 2中非高斯分布的尖峰和重尾大小, 而模型中邻域的选择即反映了图 5中系数尺度内相关性。GSM模型既可描述小波系数的边缘形状又可描述邻域系数模值间的相关性, 利用这种模型可抓住天气雷达图像中回波强度的小尺度变化特征, 同时可考虑到降水回波数据间局部关联性。

在已知变量z的基础上, GSM是一种相对高斯分布模型。为表征图 4中不同尺度小波系数间的依赖性, 这里引入多尺度高斯马尔科夫随机过程以建立不同尺度间小波系数的联系(Willsky, 2002)

(4)

式中, X₁表示尺度1的状态向量, X₂表示尺度2的状态向量, A为转移矩阵, ω是均值为0的高斯白噪声。式(4)即可反映不同尺度间系数的相关性, 且由转移矩阵A控制, A值越大, 则不同尺度间系数的相关性越强, 反之越小。

3 雷达降水回波强度数据小波域GSM插值

从雷达降水回波强度数据的小波域特征可看出, GSM模型联合高斯马尔科夫过程能很好地刻画降水数据小波域的重尾性、尺度内邻域相关性及尺度间联系性。由低分辨率雷达回波强度图像分解得到的小波系数与所建立的模型匹配, 以得到小波系数的GSM模型中各参数, 然后再由建立的先验分布模型估计更小尺度的小波系数, 最后进行小波逆变换, 则可得到更高分辨率的雷达降水回波图像, 实现小波域内的高分辨率插值过程。

更小尺度的小波系数可由式(4)中尺度间状态变量依赖关系得出, 若已知观测向量X₂, 且雷达降水回波强度数据已假定为式(1)中GSM模型, 则每个邻域精细尺度状态向量X₁的估计可近似于局部维纳线性估计

(5)

由于小波分解是不断增加幅值的过程, 小波系数尺度间将呈指数退化的形式, 式(4)并不是一个完全平稳随机过程, 因此利用式(5)估计精细尺度状态变量X₁前须先计算不同尺度系数间的近似退化关系。对于q阶矩, 高频子带内的小波系数近似呈如下关系(Deidda, 2000)

(6)

式中, k表示分解级数(如k=1—4, 表示尺度2—16 km, 图像原始分辨率为1 km), c_q表示常数因子, α_q表示尺度因子。α_q反映了小波系数随尺度的大致变化规律。图 5是个例1和个例2数据水平子带小波系数一阶矩(对应均值)和二阶矩(对应方差)随尺度变化, 其中*和о符号表示直接计算出的结果, 直线表示式(6)函数形式的拟合结果。经计算, 这两个个例的α₁分别为1.49、1.51, α₂分别为2.98、3.01。

对尺度间系数进行指数退化关系的补偿后, 可利用式(5)估计精细尺度的状态向量。由于邻域间相互重叠, 对于整幅图像, 需要从系数中估计位于邻域中心的参考系数

(7)

式中, (X)_C表示取向量X的中心值, X_1C是邻域向量X₁的中心系数值, z_MAP是概率p(z|X₂)的最大后验估计值。式(7)说明, 更高分辨率的降水回波系数的估计考虑到了相邻尺度内系数局部协方差矩阵C_U、不同尺度间的回波系数相关矩阵A以及控制降水回波强度数据分布形状和重尾性的因子z的影响, 重建的高分辨率降水回波将有效保持数据中的那些重要空间分布统计特征。

式(7)中X₂是原始低分辨率图像的分解小波系数向量, C_U与C_X2等效, 转移矩阵A可通过计算不同尺度小波系数的局部邻域相关系数矩阵得到, 为得到精细尺度邻域中心系数的估计值, 还须计算z_MAP。对于z_MAP, 可通过最大后验估计方法得到

(8)

由式(2)的GSM模型可知, p(X₂|z)呈条件高斯分布。而z的概率分布p(z)在描述降水数据重尾分布特性时通常假定为一种对数高斯分布(lgz~N(μ_z, σ_z²))(Ebtehaj, et al, 2011a)

(9)

由于E(z)=1, 则z的概率密度函数隐含有如下关系

(10)

进而可得到E(z²)=exp(σ_z²)。再由式(3)中GSM模型的关系E(X⁴)=E(z²)E(U⁴)可计算得到σ_z²。至此, 式(8)中只剩下一个未知数z, 可通过基于多维二分的数值法(Portilla, et al, 2003)求解最大后验方程组, 即可得到变量z_MAP值。最后, 将这些估计出的参数代入式(7), 可计算得到更小尺度每个邻域中心的小波系数, 再结合原始低通尺度系数, 进行逆小波变换, 则可得到更高分辨率雷达回波强度图像。

总的来说, 基于小波域GSM模型的天气雷达图像高分辨率插值主要是基于先验GSM模型的更小尺度小波系数的估计, 这又涉及原始雷达回波强度数据与先验GSM模型的匹配和参数估计, 如z、α_q、A的估计。上述算法归纳总结如下:(1)原始低分辨率雷达回波强度图像进行2级以上的小波分解, 得到不同尺度不同方向子带的小波系数; (2)基于分解得到的水平、垂直、对角方向子带的小波系数, 计算它们的概率分布及邻域协方差矩阵, 构建不同方向子带的如式(1)所示的GSM模型, 并基于式(8)利用最大后验概率(MAP)法估计GSM模型中加权系数z值, 其中z假定为式(9)所示的对数高斯分布; (3)根据不同尺度间小波系数关系计算邻域转移矩阵A, 并根据一阶矩和二阶矩关系, 拟合得到不同尺度系数的指数关系(式(6)); (4)基于步骤(2)和(3)得到的各参数, 根据式(7)利用贝叶斯估计计算得到水平、垂直、对角方向子带更小尺度的小波系数; (5)根据已有的尺度系数和小波系数以及估计得到的各方向子带的更小尺度小波系数, 如式(11)所示, 进行小波逆变换, 则可得到更高分辨率的雷达降水图像f。

(11)

式中, WT^-1表示小波逆变换, a_k表示第k级的尺度系数, d_Hi、d_Vi、d_Di分别表示第i尺度(如i=1表示尺度1、分辨率2 km, 尺度2表示分辨率4 km)的水平H、垂直V、对角向D子带小波系数。或直接将原始低分辨率图像(记为a₁)作为低通的尺度系数, 联合估计出尺度1的小波系数, 进行小波逆变换

(12)

4 个例试验

以2008年5月27日和2011年7月25日两次强降水为个例试验说明基于小波域GSM模型插值算法的有效性, 并将此方法结果与常规的双线性插值结果进行了对比。

首先以水平向子带系数估计为例, 说明GSM模型匹配中参数估计结果的有效性。图 6a中实线是基于MAP法估计的个例1中水平子带GSM模型中的加权系数z值的概率分布, 其中横坐标为lgz。由于模型中假定z为对数高斯分布, 则lgz为高斯分布。为说明估计参数z的分布, 图 6a中虚线显示了均值为μ_z、均方差为σ_z的理想高斯分布, 其中μ_z、σ_z根据输入数据计算得到。由图 6a可看出, 基于MAP法估计出的参数z近似为对数高斯分布, 与理论模型基本吻合。图 6b是利用得到的个例1水平子带先验GSM模型, 根据图 6c中原始低分辨率图像分解得到的小波系数估计得到的更高分辨率的小波系数。对比图 6b、c可看出, 更小尺度的小波系数比粗尺度的小波系数能体现出更丰富更精细的信息。

图 7是个例1插值试验结果, 其中图 7a是原始高分辨率图像(分辨率1 km×1 km), 作为结果对比中的真值, 图 7b是原始高分辨率图像进行4×4格点平均后的降分辨率图像(分辨率4 km×4 km), 作为待插值的原始低分辨率输入图像。图 7c是对图 7b中图像进行双线性插值后的结果(分辨率1 km×1 km)。从图 7可以看出, 双线性插值后图像整体变得更平滑, 但图中的大多数细节也被平滑掉, 强回波中心变少或少数消失。这与双线性插值的低通滤波效应有关。图 7d是基于小波域GSM模型高分辨率插值后的结果(分辨率1 km×1 km)。由于插值后分辨率从4 km变成1 km, 图 7中的小尺度小波系数须经过两次估计。对比图 7a、d可看出, 小波域GSM插值可恢复大多数高频细节信息, 且强回波位置与范围也与原始高分辨图像比较一致, 它的低频平滑部分与图 7b信息基本一致。由于小波域GSM插值主要是针对雷达降水图像自身特点建立模型估计出更精细尺度的小波系数, 此方法能很好地重建图中易出现变化的高频细节, 使回波结构更加精细, 凸显强回波位置, 这对强对流降水回波的识别和监测非常有用。

表 1是个例1插值试验整体像素(>0 dBz)和强回波(>40 dBz)双线性插值结果(GR_in)、小波域GSM插值结果(GR_GSM)与原始高分辨率图像(GR_ori)的统计比较, 其中熵差的计算式(13)(Gonzalez, et al, 2008)为

(13)

式中, EN_o表示原始高分辨率图像的熵值, EN_i表示插值图像的熵值, p(Z_k)表示降水回波反射率因子Z_k位于第k个直方图区间的概率。熵值反映了雷达强度回波空间分布变化或规律性, 熵值越大变化规律越不确定则越难估计或预测。而熵差反映了插值图像与原始真值图像回波强度大小空间分布特征的差异, 熵差越大, 则插值后的图像与原始真值回波分布规律差别越大, 反之越小。

从表 1中整体像素对比来看, GR_in与GR_ori和GR_GSM与GR_ori的均方根误差基本接近, 但GR_in与GR_ori的平均偏差和熵差均明显大于GR_GSM与GR_ori的结果, 说明GR_GSM与真值GR_ori的反射率因子大小更接近, 且整体分布规律一致, 而GR_in相对于GR_GSM破坏了原始降水回波GR_ori的空间分布结构特征, 如GR_in弱化了一些强回波到弱回波之间的变化梯度。GR_GSM与GR_ori的熵差接近于0, 反映了两者反射率因子大小变化的一致, 信息量的一致, GR_GSM很好地保持了原有降水回波的空间分布特点, GR_GSM相对于真值GR_ori, 具有更小的不确定性。

为更好说明小波域GSM插值能很好恢复强回波细节的优势, 表 1也给出了插值图像和真值图像大于40 dBz回波的统计比较结果。由于双线性插值的平滑效应, 其强回波像素个数明显少于真值中强回波像素个数, 而基于小波域GSM模型的插值能有效重建高分辨率细节信息, GR_GSM像素个数明显多于GR_in, 与真值GR_ori的结果比较接近。且在强回波时, GR_in与GR_ori的平均偏差更大, 而GR_GSM强回波平均值与真值的平均值几乎一致, GR_GSM与GR_ori的均方根误差也更小。在凸显强回波区这方面, 小波域GSM插值明显优于双线性插值, GR_GSM能恢复出低分辨率图像中丢失的大部分高频变化区域和强回波值, 且小波域GSM插值考虑了数据本身的特点, 重建的强回波范围、大小以及回波分布特征与原始图像更接近。

图 8及表 2是个例2的插值试验图像结果和统计比较结果, 其设置与图 7和表 1一样。个例2的性能与个例1的性能基本一致, 但由于个例2中强回波的比例更大, GR_GSM与GR_ori相对于GR_in与GR_ori的平均偏差和均方根误差优势更明显。双线性插值平滑掉了大多数强回波, 而基于小波域GSM模型的插值能恢复大多数强回波中心及变化细节, 且分布规律与原始高分辨率图像基本一致。

5 结论

天气雷达回波强度数据的小波系数反映的是降水变化、边缘等细节信息, 其具有尺度间统计自相似性, 表现出稳定的分形特征。除了尺度间的依赖性, 降水回波的局部不连续性和一些强回波奇异值等特征外, 在小波域还表现出明显的非高斯重尾分布特性。文中在分析雷达回波强度数据小波域统计特点基础上, 提出一种基于小波域GSM模型的天气雷达图像高分辨率插值方法。此方法通过建立先验模型以刻画雷达回波强度数据的重要统计特征, 然后利用此模型估计出更高分辨率的高频细节信息, 不仅考虑了原有数据邻域间相关性及不同分辨率之间的关联性, 还捕捉到了降水回波空间分布的非高斯统计特性, 插值出的回波数据能有效反映此降水回波数据的空间变化统计特征, 可突出雷达图像中强回波以及回波突变区域和回波细节结构信息。需要注意的是, 此方法估计模型参数时利用了降水回波数据的先验数学统计特征, 若原有低分辨率图像已粗糙到几乎看不到纹理、边缘等变化信息, 则很难正确地恢复出高分辨率降水回波细节结构; 且重构出的高分辨率细节特征主要反映的是局部邻域的高频统计特征, 可能不能重现完全一样的回波几何形状。

此方法不同于常规的基于空间低通滤波的天气雷达图像插值方法, 它利用了雷达回波强度数据的高通小波域系数, 关注和强调了降水回波高频变化细节信息的重建和恢复。另外, 它利用GSM模型作为降水回波高通小波系数的先验模型, 表征降水回波的空间分布统计特征, 既考虑到了降水回波的局部关联性, 又能有效反映降水回波的非高斯重尾特性, 从而能有效抓住降水分布中的强回波值及回波强度变化细节。将此插值结果应用到中小尺度对流天气过程精细探测或中尺度预报的数据同化中, 可有效改善降水回波局部结构特征的检测和预报。GSM模型为条件高斯分布, 这也为参数估计和计算带来便利, 此模型和方法也可考虑拓展应用到多传感器数据融合或数据同化系统中, 改进传统的基于线性的融合或基于高斯假定的同化方法如卡尔曼滤波法所带来的平均效应, 从而有效保持降水回波分布中具有非高斯重尾性能的极大值和奇异点, 这将十分有利于强对流灾害性天气的监测和预报。