基于粒子群優化的偏斜α穩定分布參數估計

李靜威 全厚德 崔佩璋

（軍械工程學院光學與電子工程系，河北 石家莊050003）

引 言

短波（3～30MHz）是進行遠距離無線通信的重要波段，而大氣噪聲是影響短波天波通信的重要因素之一.長期以來，該波段的信號處理研究往往假設環境噪聲為高斯噪聲，但在實際應用和研究中都發現，由于大氣放電，環境噪聲并不符合高斯噪聲的統計特性，而具有較強的沖擊特性.特別在15MHz以下，大氣噪聲沖擊特性非常明顯［1-5］，給短波通信造成了很大影響.沖擊噪聲的統計特性具有非高斯和重尾的特性，而這些特性恰恰符合α穩定（αS）分布［6］.環境噪聲特性研究是一項基礎性工作，因此其參數估計受到研究者的廣泛關注.

國內外學者對αS分布模型的參數估計提出了多種方法，主要有特征函數法（CFM）［7-8］、分位數法（QM）［9-10］、極大似然估計法（MLE）［11-13］、極大值估計法（EVM）［14-15］、矩估計法（MM）［15-16］、聯合估計（JE）［17］以及貝葉斯參數估計（BE）［18-19］等.其中，CFM、QM、MLE是通用的參數估計方法，但是αS分布獨有的非線性、無限方差、無閉式的概率密度函數（PDF）等特性給估計帶來了很多困難.EVM是基于αS分布重尾性的估計方法，但是EVM需要查表確定樣本的分組大小.JE是對以上幾種方法的綜合應用.在很多情況下，可采用整數階矩進行統計參數估計，但是αS分布的二階及高階矩是無限的，因此并不適用.分數低階矩（FLOM）的提出為αS分布的參數估計提供了特有的方法.以上方法除BE外待估參數都必須按照特定順序逐一估出，無法同時得到.然而，BE雖能同時進行多參數估計，但需提前給出3種先驗分布中的6個參數，參數設置的復雜性及隨機性對后續的估計影響較大.

粒子群優化（PSO）［20］算法是 Kennedy和Eberhart于1995年提出的一種隨機全局優化算法.該算法具有在較高維空間發現全局最優解的能力，可解決大量非線性、不可微和多峰值等復雜優化問題，實現簡單，收斂速度較快，穩定性好［21］.本文針對αS分布的非線性、無限方差及無閉式PDF的實際選取PSO算法，結合αS分布的FLOM 特性，提出了PSO-FLOM方法，實現了αS分布多個參數的同時估計.蒙特卡羅（MC）仿真實驗證明了本方法的有效性和準確性.

1 α穩定分布

α穩定分布是目前惟一符合廣義中心極限定理的分布類型，即無限多個可能方差無限大的獨立同分布的隨機變量和的極限分布是αS分布.因其沒有閉式PDF，只能由特征函數進行表征.目前，αS分布具有三種不同的參數系表征方法，分別稱為S標準參數系、S1參數系和S2參數系［22-23］.S參數系在信號處理中使用普遍，也是文獻中最為常用的［24］，因此僅討論S參數系下的參數估計.

定義1 稱隨機變量X具有一個穩定分布，如果存在四個參數 （α，β，σ，μ），其中0＜α≤2，-1≤β≤1，0＜σ，μ為實數，其特征函數具有式（1）形式：

α稱為特征指數，決定該分布脈沖特性的程度.α值越小，所對應分布的拖尾越厚；相反，隨著α值變大，所對應的拖尾變薄，且脈沖特性減弱.β稱為對稱參數，用于確定分布的偏斜程度.β＝0對應于對稱分布，稱為對稱α穩定（SαS）分布；-1≤β＜0對應左偏斜分布；0＜β≤1對應右偏斜分布.σ為尺度參數，又稱分散系數，它是關于樣本相對于均值的分散程度的度量.μ為位置參數，用于確定PDF的位置.

特別的，α＝2時，該分布為高斯分布；α＝1，β＝0時為柯西分布；α＝1/2，β＝-1時為泊松分布.除這三種特殊情況外，αS分布的PDF沒有封閉的表達式.四個指數的靈活取值使該分布具有廣泛的適應性.

2 Kuru-FLOM估計

若隨機變量X滿足

令γ＝σα，則X的分布亦可表示為

對于α≠1，偏斜αS分布的絕對FLOM理論值為

其符號FLOM理論值為

一個具有N個樣本點的隨機變量X，其絕對FLOM和符號FLOM的估計為

通過FLOM的乘積與比值，利用伽馬函數式（8）的性質，可得到一些簡單的α、β和γ的閉合表達式.

1）α估計式

2）α估計出來后，利用估值α＾和β的估計式（10）進行β的估計.

式中θ可由Sp/Ap、ApS-p、ApA-p、SpS-p等得到.但后兩種無法確定β的符號.

3）估計出α和θ后，利用下式估計γ為

并根據γ＝σα，得到σ的估值.

3 PSO-FLOM估計

3.1 PSO算法原理

PSO算法首先初始化一群隨機粒子（隨機解），每個粒子具有位置和速度兩個特征.粒子位置坐標對應的目標函數值為該粒子的適應度，算法通過適應度來衡量粒子位置的優劣.假設D維搜索空間中第i個粒子的位置和速度分別為Xi＝（xi1xi2… xiD）、Vi＝（vi1vi2… viD），在每次迭代中，粒子通過跟蹤兩個最優解來更新自己，第一個是該粒子本身最優解，即個體極值，Pi＝（pi1pi2… piD）；另一個是整個種群目前找到的最優解，即全局最優解Pg.找到兩個最優值后，第i個粒子的第j維向量根據下面兩式來更新自己的速度和位置.

式中：w為慣性權因子，它使粒子保持運動的慣性，使其有能力探索新的區域；c1和c2為正的學習因子，分別調節全局最好粒子和個體最好粒子飛行的最大步長；r1和r2為0到1之間均勻分布的隨機數.

更新過程中，粒子每維的坐標限制在 ［xmin，xmax］.根據αS分布參數定義域確定取值范圍，提出位置映射進行限制.以α∈（0，2］為例，當粒子該維的更新值為α′，α的取值為

按照同樣方法，根據β和σ值域分別進行相應的映射，保證迭代值落在各自定義域中.

3.2 目標函數確定

3.2.1 目標函數選取

根據FLOM理論式和估計式，選用絕對FLOM.偏斜αS分布模型參數估計就是使式（15）最小化的參數取值，即

式中p為FLOM的階次.

3.2.2 階數p的選擇

p≥α時，FLOM值為無窮，因此需p＜α.FLOM的協方差在p＝α/2時也是無限的［15］.根據式（15），分別計算了在β＝0.5、σ＝1、μ＝0、α＝［0.5，0.9，1.3，1.7］四種情況下，p＜α時，目標函數隨p值變化的情況.目標函數值為1時，p值分別為 ［0.17，0.32，0.50，0.74］，分 別 對 應 于 α 值 的0.34、0.36、0.38和0.44倍.以α＝1.7為例，圖1給出了兩者的關系曲線.從圖1可以看出：α/2＜p＜α時，誤差較為平穩，基本保持在1值附近；p＝α/2附近時誤差急劇增大；p＜α/2時誤差具有較明顯下降趨勢，且越接近0值下降趨勢越強.因此，進行估計時，p值至少應小于α/2.依據這一原則，p值的大小并不會影響估計算法的實質.

3.2.3 β符號的確定

圖1 目標函數隨p值變化曲線

借鑒簡便的樣本分位數法，首先確定β的符號以縮小β值的搜索范圍，提高搜索效率.若Xmd表示隨機序列的中位數，Xmin表示序列中的最小值，Xmax表示序列中的最大值，則Xmin和Xmax亦分別代表了分布左右兩端的分位數.

綜上分析，PSO-FLOM流程如圖2所示.

圖2 PSO-FLOM流程圖

表1 β的符號判斷

4 仿真實驗與分析

根據αS分布參數特點，結合PSO算法的經驗和準則，PSO-FLOM參數設置為：粒子群規模30，維數3，α、β、σ值域分別為（0，2］、［0，1］、［0，2］，粒子初始速度為0.5＊rand，學習因子c1、c1分別設為2.8和1.3，迭代次數為500.不失一般性，FLOM階次p＝0.3.PSO算法本質是一種隨機算法，即使用同樣的參數，每次求解的結果也可能不同，故進行200次MC仿真后取其均值.

表2給出了Kuru-FLOM和PSO-FLOM兩種方法的估計結果，同時給出了估計值誤差相對真值的誤差比（ER）.實驗1顯示，PSO-FLOM的估計參數精度分別提高了2.4%、4.2%、2.2%；實驗2顯示，PSO-FLOM 的估計精度分別提高了2.8%、4.4%、2.4%；實驗3顯示，PSO-FLOM 的估計精度分別提高了-0.8%、5.6%、0.2%.實驗1和實驗2在估計精度上都有了明顯的提高，而實驗3除α估計精度略有下降外，其余2個參數估計精度也有所提高.其原因為實驗3的α真值有所降低，根據FLOM階數p值的選取規律，要達到與實驗1、2相應的估計精度，p值需適當減小.總之，PSO-FOLM方法對比Kuru-FLOM方法的估計精度有所提高或相當.在多參數估計方式上，前者能同時給出多個參數估值，而后者必須依次給出參數估值.

表2 試驗結果對比

根據仿真實驗估計結果，圖3分別給出3個實驗估值與真值情況下得到的αS分布隨機序列的PDF曲線.可以看出：實驗1和實驗2中由PSOFLOM方法的估計結果產生的PDF曲線更加接近真值曲線，說明其是有效的，且較為精確；實驗3中兩種估計方法所得的PDF曲線相當.因為αS分布由4個參數綜合決定，每個參數存在偏差都會對估計曲線產生影響，即使估計精度較高，估值曲線與真值曲線也并不會完全重合，在α值較小時，綜合影響的程度更加明顯.

圖3 仿真實驗結果PDF曲線對比

5 結 論

提出了基于粒子群優化的偏斜α穩定分布參數估計，并在對稱參數β符號確定和FLOM階次p選取上給出了相應分析，進而確定了算法的目標函數和PSO-FLOM估計方法的流程，最后通過仿真實驗驗證了該方法的有效性.仿真實驗顯示該方法PSO參數設置簡單，克服了Kuru-FLOM估計方法必須逐步給出α、β、σ估計值的不足，能夠同時得到3個參數的估計值，估計精度有所提高或與之相當.雖然偏斜αS分布的FLOM參數估計研究還限于位置參數μ＝0的情況，但無疑PSO-FLOM給偏斜α穩定分布參數估計提供了一種簡單有效的新選擇.

［1］KRAMER J D R，WILLIAMS R T.High frequency atmospheric noise mitigation［C］//IEEE International Conference on Acoustics，Speech，and Signal Processing.Adelaide，April 19-22，1994：101-106.

［2］LAUBER W R，BERTRAND J M.HF atmospheric noise levels in the Canadian arctic［J］.IEEE Trans on Electromagnetic Compatiblity，1994，36（2）：104-109.

［3］SHINDE M P，GUPTA S N.A model of HF impulsive atmospheric noise［J］.IEEE Trans on Electromagnetic Compatiblity，1974，16（2）：71-75.

［4］AGARWALA R A.Frequency characteristics of atmopheric noise［J］.IEEE Trans on Electromagnetic Compatiblity，1983，25（2）：119-120.

［5］黃德耀，王聚杰，楊維富，等.福建地區大氣無線電噪聲的某些特征［J］.電波科學學報，2002，17（6）：650-655.HUANG Deyao，WANG Jujie，YANG Weifu，et al..Atmospheric radio noise characteristics of Fujian regions［J］.Chinese Journal of Radio Science，2002，17（6）：650-655.（in Chinese）

［6］SHAO M，NIKIAS C L.Signal processing with fractional lower order moments：stable processes and their application［J］.Proceedings of the IEEE，1993，81（7）：986-1010.

［7］KOUTROUVELIS I A.Regression-type estimation of the parameters of stable laws［J］.Journal of the American Statistical Associstion，1980，75（4）：918-928.

［8］WILLIAMS S K.Characteristic function based estimation of stable distribution parameters［M］.Berlin：Birkhauser，1998.

［9］ROLL F F E.Parameter estimations for symmetric stable distribution［J］.Journal of the American Statistical Associstion，1971，66（2）：331-338.

［10］MCCULLOCHJ H.Simple consistent estimators of stable distributions［J］.Communications in Statistics-Simulation and Computation，1986，15 （4）：1109-1136.

［11］MITTNICS，RAVHEV S T，DOGANOGLU T.Maximum likelihood estimation of stable paretian models［J］.Mathematical and Computer Modelling，1999，29（1）：275-293.

［12］BODENSCHATZ J S，NIKIAS C L.Maximum-likelihood symmetric alpha-stable parameter estimation［J］.IEEE Trans on Signal Processing，1999，47（5）：1382-1384.

［13］NOLAN J P.Maximum likelihood estimation and diagnostics for stable distributions［D］.Washington：American Univeristy，2002.

［14］TSIHRINTZIS G A，NIKIAS C L.Fast estimation of the parameters of alpha-stable impulsive interference［J］.IEEE Trans on Signal Processing，1996，44（6）：1492-1503.

［15］KURUOGLU E E.Dendity parameter estimation of skewedα-stable distributions［J］.IEEE Trans on Signal Processing，2001，49（10）：2192-2201.

［16］MA X Y，NIKIAS C L.Parameter estimation and blind channel identification in impulsive signal environments［J］.IEEE Trans on Signal Processing，1995，43（12）：2884-2897.

［17］LI X J，JIN L W，LI X T.Joint parametrs estimation for general alpha stable random noise［C］// 2008 Congress on Image and Signal Processing.Sanya，May 27-30，2008：465-469.

［18］LOMBARDI M J.Bayesian inference forα-stable distributiongs：A random walk MCMC approach［J］.Computational Statistics ＆ Data Analysis，2007，51（5）：2688-2700.

［19］郝燕玲，單志明，沈 鋒.基于自適應Metropolis算法的α穩定分布參數估計［J］.系統工程與電子技術，2012，34（2）：236-242.HAO Yanling，SHAN Zhiming，SHEN Feng.Parameter estimation ofα-stable distributions based on adaptive Metropolis algorithm［J］.Systems Engineering and Electronics，2012，34（2）：236-242.

［20］KENNEDY J，EBERHART R.Particle swarm optimization［C］//Proc of the IEEE International Conference on Neural Networks.Perth，27Nov 1995-01 Dec 1995：1942-1948.

［21］曹運華，吳振森，齊利華，等.粒子群算法在BRDF模型參數優化中的應用［J］.電波科學學報，2008，23（4）：765-768.CAO Yunhua，WU Zhensen，QI Lihua，et al..Application of particle swarm optimization algorithms to parameter optimization of BRDF model［J］.2008，23（4）：765-768.（in Chinese）

［22］NOLAN J P.Parameterizations and modes of stable distributions［J］.Statistics ＆ Probability Letter，1998，38（2）：187-195.

［23］RAFAL W.On the chambers-mallows-stuck method for simulating skewed stable random vairables［J］.Statistics ＆ Probability Letter，1996，28（2）：165-171.

［24］李旭濤，朱光喜，王首勇，等.Alpha穩定分布的參數表征及仿真［J］.信號處理，2007，23（6）：814-817.LI Xutao，ZHU Guangxi，WANG Shouyong，et al..Parameterizations and simulation of alpha stable distribution［J］.Signal Processing，2007，23（6）：814-817.（in Chinese）