引言

众所周知，VAD技术在水平风场为线性分布的假定条件下，从一阶谐波系数中得到平均风向风速，从二阶谐波系数得到平均水平弹性形变和切变形变，并且在任何水平呈线性和非线性分布的风场，VAD技术都可以提取出零阶谐波项的系数a₁^[1~4]，已经过Srivastava^[5]和Matejka^[6]证明:任何含有小于四阶谐波的风场，零阶谐波项的系数a₁本身包含两项，即平均水平散度和降水质点相对于大气的垂直速度V_f, 如式 (1) 所示:

由上式可见，由a₁不能直接得到水平散度，而且，其中V_f实际上是大气垂直速度W和有效照射体积内不同大小粒子在静止大气中的平均下落末速度W₀之和，即使可准确得到V_f, 但仍不能得到W。为此，许多人都想从a₁中提取出水平散度，通过连续方程估计出大气的垂直速度W，其中Srivastava^[5]提出了提取水平散度的方法，即EVAD技术，他假定在某一水平层内水平散度和垂直速度V_f不变，用最小二乘方法，从VAD分析中得到的零阶谐波项a₁中分离出水平散度和降水质点相对于大气的垂直速度V f, 但由于此方法假设条件有问题，且违背了回归变量的方差在独立变量的范围内为一常数的最小二乘分析的基本假定，故分离水平散度的效果不好。为此，Matejka^[6]在Srivastava的基础上提出了用带有权重的最小二乘法分离水平散度的EVAD技术，分离的效果有所改善，但在信噪比较小、速度谱宽较宽的条件下，还存在明显的误差，为此，本文提出了另一个权重函数，应用模拟信号技术，对模拟的多普勒速度场进行数值实验，数值实验表明，本方法在信噪比较小、速度谱宽较宽的条件下，对水平散度有明显的改善。

1.   基本原理及计算方法

1.1   EVAD算法

EVAD分析分两步:①进行EVAD分析的第一步是在含有0，1，2，3阶谐波的非线性风场的假设下，应用单多普勒雷达采集的速度资料，在每一个距离圈上对V_r进行谐波分析，可得不同仰角上每个距离圈的零阶谐波项a₁及a₁的方差var (a₁)。② 3种算法的最小二乘回归求平均水平散度。EVAD分析假设在水平薄层内，水平散度和垂直速度为常数，故EVAD分析是在一定厚度的水平层上进行的，薄层的厚度至少应该与沿着波束的数据间隔一样厚，通常Δh取在150~1000 m (本文Δh为300 m)，将速度资料按Δh垂直分层。根据公式 (1)，

其中g₁=1，g₂=-2sinα/(rcosα)。利用在每一水平层内不同仰角和半径的资料以及第一步中得到的每一水平层内不同仰角上每一个距离圈上的零阶谐波项a₁及a₁的方差var (a₁)，对式 (2) 进行最小二乘回归得到的系数p_i的解为:

用式 (3) 便能计算出每一水平层的p₁，即每一水平层的平均水平散度

Srivastava等人提出的EVAD算法中，系数p_i的解式 (3) 是对式 (2) 进行不带权重的最小二乘回归得到的，式 (3) 中:

Matejka和本文建议的EVAD算法中，系数p_i的解式 (3) 都是对式 (2) 进行带权重的最小二乘回归得到的，式 (3) 中:

M为这一层内距离圈的总数，γ_k为这一层内第k个距离圈的权重。但这两种算法的权重函数不同。Matejka的权重函数为:

其中，r_k=R_kcosα_k, R_k为这一层内第k个距离圈的径向距离，α_k为这一层内第k个距离圈所在的高度仰角，n_k为这一水平层内第k个距离圈所在仰角上的距离圈总数。

本文在Matejka的权重函数的基础上，对其中的n_k和var (a₁) 进行了两点改善，即，用Δz_k代替n_k, 且对Thomas的权重中的var (a₁) 设定一个阈值σ_k。改善后的权重函数为:

其中，Δz_k=|h_k-h′|/ h₀，h_k为这一水平层内第k距离圈所在的垂直高度，h′为这一层的中心所在的高度，h₀为薄层厚度的一半，即h₀=Δh/2，当Δ_zk=0(即h_k=h′) 时，这一层内第k距离圈的权重就等于第k圈所在仰角上第k-1圈和第k +1圈的权重的算术平均，即γ_k=(γ_k-1+γ_{k +1})/2。这里有两个原因表明为什么要对Matejka的权重函数中的n_k和var (a₁) 进行改善。第一，Srivastava的EVAD分析是在某一水平层内进行的，很明显，在有一定厚度的水平层内，低仰角探测的距离圈数比高仰角的多，所以在进行最小二乘时，低仰角的影响大。为此，Matejka将n_k考虑到权重中去，但1/n_k只是对一层内同一仰角上的距离圈的算术平均，而不是对一层内所有距离圈的加权平均，所以对其进行了改善。而且，由于Srivastava假定在某一水平层内，水平散度保持不变，然而在实际中，水平层内的水平散度是变化的，故为了补偿假定条件的影响，在水平层内，取某一高度的水平散度做为标准，本文以这一层的中心所在高度上的水平散度为标准，由距离权重1/Δ_z作为每一个距离圈的权重，使靠近水平层中心的距离圈上的权重大，而远离中心的距离圈上的权重小，从而补偿了假定的影响。第二，虽然Matejka把方差var (a₁) 考虑到权重中去，但是当某距离圈上的方差很大时，即表明此距离圈上的数据质量太差，则这个距离圈上的权重非常小，在这种情况下，完全可以将这个距离圈从此层中去掉。故将数据质量差的距离圈从此层中去掉，然后再进行有权重的最小二乘回归。所以，对Matejka权重中的var (a₁) 设定一个阈值σ_k, 即:

其中，var (a₁)k为某一水平层内第k个距离圈上a₁的方差，为这个水平层内所有距离圈上a₁的方差的算术平均。对这一层内每一个距离圈都求出Δ_{var (a1)k}=将这一层内的距离圈保留，而把Δ_{var (a1)k}> σ_k的距离圈都删除后，再对这一层内剩下的距离圈，用本文建议的权重函数进行最小二乘回归，得到这一层的平均水平散度。本文将这种方法称为汤-陶法，简称T-T法。

本文通过模拟降水回波信号的方法来验证此方法，故下面根据文献[7]介绍一下模拟降水回波信号的基本理论。

1.2   回波信号模拟的基本理论

①模拟信号的采样

降水回波的功率谱密度通常呈高斯谱分布，而谱分析的信号是离散的谱信号，必须对原谱进行离散化，用采样公式可得到离散谱序列{S_n}。以V_r=11 m/s, 即径向速度V_r所对应的多普勒频率f₀=120 Hz (f₀=2 V_r/ λ，波长λ=10 cm)，频率谱宽σ_f=40(或速度谱宽σ_v=2 m/s)，Nyquist频率F_n=256 Hz为例，图 1是在原谱的基础上取信噪比SNR=40 dB时得到的。

②模拟信号加噪

为了使模拟出的时域信号尽可能符合实际情况，需要对采样后的离散谱附加一定信噪比的白噪声，用加噪公式可得到加噪后的离散功率谱{P_n}，如图 1所示 (其中P为加噪声后的离散功率谱值，P_max为最大的功率谱值)。

图  1  径向速度V_r为11 m/s (即对应的多普勒频率f₀为120 Hz, 速度谱宽σ_v为2 m/ s) 时, 附加40 dB信噪比 (SNR=40 dB) 白噪声后的离散功率谱密度图 (其中Nyquist频率F_n=256 Hz)

下载: 全尺寸图片 幻灯片

③附加相位谱

对实谱序列附加上相位谱，考虑到噪声的随机性，附加的是[0，2π]区间上概率均匀分布的随机相位谱，得到的复谱序列，再进行反傅立叶变换，得到回波的时域信号。

①~③的具体做法和其中用的公式见文献[7]。

④用加噪后的时域信号，根据文献[7]中的脉冲对处理法 (PPP法) 得到加噪后的径向速度

模拟风场的每一个径向速度所对应的多普勒频率f_d(f_d=2 V_r / λ) 做为模拟回波信号的平均多普勒频率f₀，按照上述①~③步骤对模拟风场的每一个径向速度加不同信噪比和不同速度谱宽的噪声，模拟出每一个径向速度的带噪声的回波信号，最后按④将模拟的回波信号还原成径向速度，从而模拟得到带有噪声的径向速度。

2.   数值模拟试验及结果分析

2.1   数值模拟试验

为比较Srivastava、Matejka和本文建议的T-T法得到的水平散度的精度，本文进行了数值模拟试验。

首先，考虑到EVAD算法中径向速度展开的傅氏级数中含有方位角θ的0，1，2，3阶谐波项，本文设计模拟风场含有0，1，2，3阶谐波，此风场的水平风场呈非线性分布。令u, v分量符合含有θ的0，1，2，3阶谐波的非线性分布条件，将它们分别对雷达中心进行泰勒展开，则有

式中的θ、β分别为方位角、高度仰角，R为径向距离，u₀，v₀以及u, v的各阶偏导数均在扫描中心取值。模拟风场的特征参量如表 1所示。将上述公式带入下式:

表  1  模拟含0, 1, 2, 3阶谐波的非线性风场的特征参量

下载: 导出CSV 

| 显示表格

再由表 1的特征参量可得到全方位角上径向速度模拟场。此模拟风场是按美国标准探测模式即在14个高度仰角β上进行模拟，方位角按采样间隔为1°的均匀采样进行模拟，并且在每条径向上每隔150m模拟一个径向速度资料库，共100个资料库。

其次，本文设计的非线性风场，是通过对某一固定距离圈上对V_r(θ) 作谐波分析的基础上得到的，所以，此模拟风场上的每一个径向速度都在VAD和EVAD显示的速度-方位曲线上，此径向速度不含有雷达的机内噪声。为了使模拟风场更加接近实际风场，按照1.2节的步骤对上述模拟风场的每一个径向速度加不同信噪比和不同速度谱宽的噪声，模拟得到了带有噪声的径向速度场。从而可以通过数学模拟的方式更好地检验本文建议的EVAD方法求得的水平散度的精度。由于计算量太大，本次模拟仅对第一水平层的径向速度加噪声，模拟中用到的假设:①脉冲重复频率PRF为512 Hz, 即Nyquist频率F_n为256 Hz; ②雷达发射波的波长为10 cm; ③信噪比SNR的取值范围为0~40 dB; ④PPP法中取64个脉冲对。

再次，对模拟的第一水平层内带有噪声的径向速度按1.1节中的步骤分别用Srivastava、Matejka和本文建议的T-T法求出平均水平散度。通过按1.1节的步骤对式 (2) 求解，得到的最小二乘回归的系数即可得到第一水平层内的平均水平散度。最后将3种方法得到的平均水平散度，与表 1中模拟风场所设定的平均水平散度 (div V_h=2×10^-5s^-2) 相比较，通过平均水平散度的相对误差ε_D(ε_D=((D′-D)/D)×100 %，D为表 1中模拟风场所设定的平均水平散度，D′为对模拟的非线性风场加噪声后得到径向速度场进行EVAD分析得到的平均水平散度) 来检验本文建议的EVAD算法求得的水平散度的精度。

2.2   不同信噪比SNR、不同速度谱宽σ_v条件下，3种EVAD算法的结果分析

模拟回波信号由平均多普勒频率、频率谱宽σ_f(或速度谱宽σ_v)、信噪比SNR这3个参数决定，故速度谱宽、信噪比的变化将引起径向速度以及用EVAD方法得到的水平散度的变化。

(1) 速度谱宽σ_v一定 (σ_v=2 m/s)，信噪比SNR变化

本文分别对模拟的带有SNR分别为35，20，15 dB噪声的3个径向速度资料进行EVAD分析，结果见表 2:当SNR减小时，3种EVAD算法得到的第一层内的平均水平散度以及它的相对误差ε_D都随之增加，Srivastava方法的ε_D都大于10 %，而本文建议的方法 (T-T法) 的ε_D都控制在10 %以内，则说明T-T法受信噪比的影响最小; 当SNR相同时，水平散度的相对误差由大到小依次为:Srivastava方法最大，Matejka方法次之，T-T法最小; 当SNR=15 dB时，Matejka方法的ε_D比Srivastava方法的ε_D减小了6.5 %，但Matejka方法的ε_D仍大于10 %，而T-T法的ε_D不但比Srivastava方法的ε_D减小了19.6 %，且ε_D小于10 %，这说明当SNR越小，即噪声大时，T-T法提高水平散度的精度的幅度越大，显然，信噪比小时，本文建议的方法 (T-T法) 更能有效的提高水平散度的精度。

表  2  不同SNR条件下, 3种EVAD算法的结果对比

下载: 导出CSV 

| 显示表格

(2) 信噪比SNR一定，速度谱宽σ_v变化

在SNR=20 dB，σ_v=6，5，4，3，2 m/s的条件下，模拟的5个带有噪声的径向速度场，分别进行EVAD分析，结果如表 3所示:当速度谱宽增大时，3种方法得到的第一层内的平均水平散度以及它的相对误差ε_D都随之增加，但本文建议的方法 (T-T法) 得到的ε_D基本上控制在10 %以内，说明本文建议的方法受谱宽的影响最小; 当谱宽相同时，水平散度的相对误差由大到小依次为:Srivastava方法最大，Matejka方法次之，T-T法最小; 谱宽越宽，本文建议的方法 (T-T法) 提高水平散度精度的幅度越大，显然，谱宽大时，本文建议的方法更能有效的提高水平散度的精度。

表  3  不同谱宽条件下, 3种EVAD算法的结果对比

下载: 导出CSV 

| 显示表格

2.3   非全方位均匀采样和非全方位不均匀采样的条件下，3种EVAD算法的结果分析

在实际境况下，测站周围的降水回波分布不均匀，同一距离圈上存在缺测的情况; 另外在资料使用前必须做质量检验、剔除奇异点，这样也会出现一些空隙、缺测区域。为此，模拟了两种缺测区域:随机取的非连续性缺口 (代表积累的V_r(θ) 缺测区) 和随机取的连续性缺口 (代表无回波区缺口)，即非全方位不均匀采样和非全方位均匀采样。从而使模拟风场的速度资料更加接近真实情况。

(1) 随机取非连续性的缺口即非全方位不均匀采样

对模拟的带有SNR为35 dB噪声的径向速度资料，随机取非连续缺口分别为10°、20°、…、60°后再进行EVAD分析，结果如表 4所示:当缺口增大时，3种EVAD算法得到的平均水平散度都随之增加，其中Matejka方法和本文建议的方法的平均水平散度缓慢增加; 当缺口相同时，本文建议的方法 (T-T法) 的水平散度最靠近模拟风场设定值 (div V_h=2 ×0^-5s^-2)，Matejka次之，Srivastava偏离模拟风场设定值的程度最大。显然，缺测区对T-T法得到的水平散度的影响小，且缺口越大，T-T法提高水平散度精度的幅度越大，即本文建议的方法 (T-T法) 更能有效的提高水平散度的精度。

表  4  不同非连续性缺口的条件下, 3种EVAD算法的结果对比

下载: 导出CSV 

| 显示表格

(2) 随机取连续性的缺口即非全方位均匀采样

对模拟的带有SNR为35 dB噪声的径向速度资料，随机取4个位置不同的20°连续性缺口，分别对带有不同位置的缺口的速度资料进行EVAD分析，结果如表 5所示:当缺口位置相同时，本文建议的方法 (T-T法) 的水平散度最靠近模拟风场设定值，Matejka次之，Srivastava偏离模拟风场设定值的程度最大。显然，无回波区对本文建议的方法 (T-T法) 得到的水平散度的影响小。

表  5  缺口范围相同但位置不同的条件下, 3种EVAD算法的结果对比

下载: 导出CSV 

| 显示表格

3.   结论

通过理论分析和对模拟资料的计算分析，结果表明:

(1) 速度谱宽σ_v一定，与Srivastava和Matejka的EVAD方法相比，本文建议的方法得到的水平散度受信噪比的影响最小，并且信噪比越小，本文建议的方法提高水平散度的精度的幅度越大，即信噪比小时，本文建议的方法更能有效的提高水平散度的精度。

(2) 信噪比SNR一定，与Srivastava和Matejka的EVAD方法相比，本文建议的方法得到的水平散度受速度谱宽的影响最小，并且谱宽越宽，本文建议的方法提高水平散度的精度的幅度越大，即谱宽大时，本文建议的方法更能有效的提高水平散度的精度。

(3) 与Srivastava和Matejka的EVAD方法相比，缺测区或无回波区对本文建议的方法得到的水平散度的影响最小，且本文建议的方法得到的水平散度精度高。