最大熵法在风廓线雷达谱分析中的应用研究

第３　２卷，第４期 光谱学与光谱分析Ｖｏｌ．３２，Ｎｏ．４，ｐｐ１０８５－１０８９２　０　１　２年４月 Ｓｐｅｃｔｒｏｓｃｏｐｙ　ａｎｄ　Ｓｐｅｃｔｒａｌ　Ａｎａｌｙｓｉｓ　Ａｐｒｉｌ，２０１２　最大熵法在风廓线雷达谱分析中的应用研究胡明宝１，２，郑国光１，３，张培昌１１．南京信息工程大学，江苏南京　２１００４４ ２．解放军理工大学气象学院，江苏南京　２１１１０１３．中国气象局，北京　１０００８１摘　要　为了解现代谱分析技术在风廓线雷达中应用的可行性，通过采用实测的和模拟的风廓线雷达回波信号，对比研究了ＦＦＴ法与最大熵法的谱分析效果。

结果表明：（１）当回波信号比较强时，两者都可以得到较好的谱分析效果；但是当回波信号较弱时，最大熵法分析效果优于ＦＦＴ法，最大熵法对地杂波具有较好的抑制能力。

（２）最大熵谱比较光滑，表明最大熵法对随机白噪声也有一定的抑制作用。

（３）最大熵法的递推阶数对谱分析结果有一定影响，最终预测误差准则确定的递推阶数一般偏小，采用１５阶的递推阶数进行最大熵法分析取得了较好结果。

由于风廓线雷达回波通常都较弱，因此研究结果可望用于改善信号处理效果。

关键词　风廓线雷达；ＦＦＴ法；最大熵法；谱分析中图分类号：ＴＮ９５９．４ 文献标识码：Ａ ＤＯＩ：１０．３９６４／ｊ．ｉｓｓｎ．１０００－０５９３（２０１２）０４－１０８５－０５　收稿日期：２０１１－０７－１７，修订日期：２０１１－１０－０８　基金项目：国家自然科学基金项目（４１００５０１８）资助　作者简介：胡明宝，１９６４年生，解放军理工大学教授 ｅ－ｍａｉｌ：ｈｕｍｂ２００５＠ｓｉｎａ．ｃｏｍ引　言 风廓线雷达采用多普勒收发体制，通过接受湍流微团对雷达波的散射回波，来测量风随高度的分布。

由于风廓线雷达的波束比较宽，在水平和垂直方向都有较大的伸展尺度，因此在雷达的每一个采样体积内，都包含有很多个湍流微团，雷达回波信号是由它们对雷达波共同作用的结果，因此需要对回波信号进行多普勒频谱分析。

目前风廓线雷达信号处理中采用的谱分析方法是周期图法［１］，这种方法不仅计算量大，而且去杂的效果很有限［２］。

因为周期图法主要利用快速傅氏变换（ｆａｓｔ　Ｆｏｕｒｉｅｒ　ｔｒａｎｓｌａｔｅ，ＦＦＴ）将信号由时域变换到频域，一旦噪声、杂波与雷达回波信号混杂在一起时，就会严重污染谱分析结果，难以得到有用的雷达回波信息。

１９９７年小波变换技术开始应用于去除风廓线雷达信号中的地杂波和间歇性杂波［３－６］，取得了一定的效果。

应用小波分析的主要困难在于需要找到一组合适的小波基，要对杂波环境具有自适应性，因为小波基及其阈值的针对性太强，限制了小波分析在风廓线雷达实际探测中的应用效果。

最大熵谱分析法具有一定的自适应性，且有更高的谱分辨率和估计精度［７］。

最大熵法已在地震波信号分析中得到应用［８］，并在合成孔径雷达、跟踪雷达等的信号分析领域开展了研究工作［９，１０］，而在大气科学中主要应用于气候变化和天气分析［１１，１２］，在气象雷达信号处理中尚未见到相关的研究报道。

气象雷达由于探测的是分布式软体目标，回波信号有其自身特点，因此有研究的必要。

本研究采用实测的和模拟的风廓线雷达时域回波复数信号（即（Ｉ，Ｑ）信号），对比分析了ＦＦＴ法、最大熵法的谱分析效果，探索提高风廓线雷达谱分析质量的方法。

１　最大熵法１．１　基本原理最大熵法用于谱分析法，相当于一个具有最小相位的预测误差滤波器［１３］，该滤波器的各阶系数ａ（１，Ｍ），ａ（２，Ｍ），…，ａ（Ｍ，Ｍ）和滤波器的输出功率ＰＭ，共同确定了输入数据序列｛Ｘ（ｎ）｝的功率谱密度函数Ｓ（ｆ），也称为最大熵谱。

Ｓ（ｆ）＝ＰＭΔｔ１－∑Ｍｍ＝１ａ（ｍ，Ｍ）ｅｘｐ（－ｊ２πｆｍΔｔ）２（１） 所谓“最小相位滤波器”是指该滤波器传递函数的所有零点都落在单位圆内，在单位圆外和圆上无极点和零点。

因此，最小相位滤波器是一种使输出信号的相位延迟达到最小的线性滤波器，这种滤波器只改变信号的相位谱，不改变信号的振幅谱，也就不改变信号分析的功率谱密度函数。

Ｂｕｒｇ在使滤波器的前向与后向预测误差能量之和为最小的约束下，给出了各阶滤波器系数的递推算法。

该算法直接利用观测数据序列｛Ｘ（ｎ）｝来计算滤波器系数，所需计算时间少，因此获得了广泛应用。

递推公式如下ｆＫ（ｎ）＝ｆＫ－１（ｎ＋１）－ａ（Ｋ，Ｋ）ｂＫ－１（ｎ）（２）ｂＫ（ｎ）＝ｂＫ－１（ｎ）－ａ＊（Ｋ，Ｋ）ｆＫ－１（ｎ＋１）（３）ａ（ｋ，Ｋ）＝ａ（ｋ，Ｋ－１）－ａ（Ｋ，Ｋ）ａ＊（Ｋ－１，Ｋ－１）（４）ａ（Ｋ，Ｋ）＝２∑Ｎ－１－Ｋｎ＝０ｂ＊Ｋ－１（ｎ）ｆＫ－１（ｎ＋１）∑Ｎ－１－Ｋｎ＝０（｜ｂＫ－１（ｎ）｜２＋｜ｆＫ－１（ｎ＋１）｜２）（５）其中Ｋ＝１，２，…，Ｍ。

ｋ＝１，２，…，Ｋ－１。

ｎ＝０，１，２，…，Ｎ－１。

且令ｆ０（ｎ）＝ｂ０（ｎ）＝Ｘ（ｎ），ａ（０，０）＝１。

Ｋ阶预测误差滤波器输出功率的递推公式为Ｐ０＝１Ｎ∑Ｎ－１ｎ＝０｜Ｘ（ｎ）｜２（６）ＰＫ＝ＰＫ－１（１－｜ａ（Ｋ，Ｋ）｜２）（７） 式（２）—式（７）组成了Ｂｕｒｇ递推算法的计算公式，随着Ｋ从１逐渐递推到某个合适的最高阶Ｍ，用上述公式便可以计算出所有的预测误差滤波器系数，然后用式（１）计算出最大熵谱。

１．２　递推阶数的确定方法在上述递推过程中，何时中止递推是一个很重要的问题，它不仅影响计算时间，还关系到谱估计的质量，中止递推时的Ｋ被称为递推阶数（也称预测误差滤波器阶数）。

阶数小，计算时间短，但是谱会被严重平滑，降低了谱的分辨率；而阶数过大，谱分辨率会提高，但是计算时间长，所以有个最佳取值问题。

目前一般采用最终预测误差（ｆｉｎａｌ　ｐｒｅ－ｄｉｃｔｉｏｎ　ｅｒｒｏｒ，ＦＰＥ）准则来确定递推阶数。

ＦＰＥ准则的计算式为ＦＰＥ（Ｋ）＝Ｎ＋Ｋ＋１Ｎ－Ｋ－１ＰＫ（８）式中ＰＫ为预测误差滤波器的输出功率。

由式（７）可知ＰＫ随Ｋ的增加而减少，而式（８）中右边前半部分随Ｋ的增加而增加，则ＦＰＥ（Ｋ）必在某个Ｋ值处取得最小值，此时的Ｋ值即为所要确定的阶数。

２　实测数据分析 在雷达厂家的协助下，从对流层风廓线雷达的信号处理分系统中，采集了Ｉ和Ｑ信号。

图１显示的是雷达中模式指向时，探测的第一个距离库的（Ｉ，Ｑ）采样序列。

图１中上面一排为Ｉ信号，下面一排为Ｑ信号，横坐标为采样点数，纵坐标表示信号的强度，对信号强度进行了归一化处理，显示的是相对于整个采样序列中信号强度最大值的百分数，Ｉ和Ｑ的最大值分别为８８和１０３。

利用图１所示的时域信号，分别进行了ＦＦＴ法与最大熵法谱分析，结果见图２。

其中最大熵法采用ＦＰＥ准则来确定滤波器阶数，然后将各阶滤波器系数代入式（１）求得最大熵谱。

图２中横坐标表示多普勒速度，最大不模糊速度为２０．１ｍ·ｓ－１，纵坐标表示探测高度（单位ｍ）。

Ｆｉｇ．１　Ｔｈｅ　ｔｉｍｅ　ｓｉｇｎａｌ　ｏｆ（Ｉ，Ｑ）Ｆｉｇ．２　Ｔｈｅ　ｃｏｍｐａｒｉｓｏｎ　ｏｆ　ＦＦＴ　ａｎｄ　ＭＥＭＦｉｇ．３　Ｔｈｅ　ｄｅｐｅｎｄｅｎｃｅ　ｏｆ　ＦＰＥ（Ｋ）ｖｓ　Ｋ６８０１光谱学与光谱分析 第３２卷 对比图２中两幅图可以看出，除了在３　１５０～４　３５０ｍ高度之间外，两种谱分析方法的结果都能够显示出谱峰，且上下高度间具有连续性，但是最大熵法在各高度的谱平滑很严重，这说明ＦＰＥ准则确定的滤波器阶数偏低。

对ＦＰＥ准则计算结果随Ｋ的变化情况进行了分析，结果见图３。

图３中横坐标表示Ｋ值，纵坐标为ＦＰＥ（Ｋ）的值。

可以看出在Ｋ＝４时，ＦＰＥ（Ｋ）取得最小值。

对风廓线雷达探测的４０个高度层的时域数据都分别采用ＦＰＥ准则，统计分析了所确定的滤波器阶数分布情况（图略），发现ＦＰＥ准则确定的阶数主要在２～５之间。

同时还发现，当回波信号强时，ＦＰＥ准则确定的阶数就大，最大熵法得到的功率谱峰就明显，即ＦＰＥ准则确定的阶数大小与信噪比有关。

在弱信号时，ＦＰＥ准则确定的阶数太低，使得计算出的功率谱被严重平滑，不利于信号谱峰显示出来。

为此，不采用该准则来确定阶数，而是对设定的不同滤波器阶数时，最大熵法分析结果进行了运行试验，分析结果见图４。

图４中坐标与图２相同。

Ｆｉｇ．４　Ｔｈｅ　ｃｏｍｐａｒｉｓｏｎ　ｏｆ　ｄｉｆｆｅｒｅｎｔ　ｉｔｅｒａｔｉｖｅ　ｓｔｅｐｓ　ｉｎ　ＭＥＭ 图４表明，达到１０阶以上时，最大熵法取得了比较好的效果，谱峰明显，且表现出了上下高度间的连续性。

图２中在３　１５０～４　３５０ｍ高度之间信号谱峰不明显的情况，在图４（ｃ）中有明显改善，这表明最大熵法有优越性。

风廓线雷达的信号通常都很微弱，ＦＰＥ准则确定的阶数常常偏低，经过反复试验，根据谱分析效果，认为采用指定１５阶的滤波器阶数进行最大熵法分析是一个合理的可接收值。

３　数值模拟检验 以上在讨论ＦＦＴ法、最大熵法谱分析效果时，采用的是实测时域数据，因为无法得知谱分布的真实情况，因此只能根据分析出的谱峰的明显程度，以及上下高度层之间谱峰位置的连续性来判断某种方法的优劣。

为了有一个客观对比标准，进行了数值模拟试验，介绍如下。

３．１　时域信号的模拟采用Ｓｉｒｍａｎｓ等［１４］介绍的方法来产生模拟的功率谱及其相应的时域信号，即首先给出已知谱宽的高斯谱分布函数，对其进行量化后，加上具有一定信噪比的随机白噪声，得到比较符合实际情况的离散谱，称为实谱。

对实谱序列附加上在［０，２π］区间上概率均匀分布的随机相位谱，得到复谱序列，最后对复谱序列进行ＦＦＴ逆变换，得到模拟的时域信号。

考虑到风廓线雷达探测时，不仅有随机白噪声的影响，还常常受到地杂波的影响，地杂波常出现在零频附近，容易和湍流信号混在一起，为了与实际情况尽可能相似，还要在零频附近迭加杂波谱。

文献［１５］指出大气湍流和各类降雨回波的多普勒速度谱宽值一般大于１ｍ·ｓ－１，有树木的山脉和城市建筑等杂波频谱的谱宽一般小于０．５ｍ·ｓ－１，而气象目标回波谱和杂波谱都可以采用高斯分布来表示。

据此，进行了湍流回波信号谱和地杂谱的混合模拟。

具体做法是：首先产生一个谱峰位于１ｍ·ｓ－１处谱宽也为１ｍ·ｓ－１的湍流信号谱，再产生一个峰值位于０．２ｍ·ｓ－１处谱宽为０．４ｍ·ｓ－１的地杂波谱。

对这两个谱采用一样的离散点数一样的谱线间隔，进行离散采样。

由于杂波和信号的强度一般不一样的，因此在谱离散之后，设置杂波的强度为１０ｄＢ，湍流信号谱的强度分别为－５，０和５ｄＢ，进行谱幅度倍增，然后对两个谱在相同的频点处的值进行迭加，得到实谱，结果如图５所示。

在图５中，横坐标表示多普勒速度，最大不模糊速度为２５．６ｍ·ｓ－１。

图５中第一排为加白噪声以后的谱，第二排为原高斯谱。

由于在不同的频点处，所加噪声是随机的，因此谱上脉动的特点很明显，与风廓线雷达实测谱很相似。