楊立儒 劉永祥 楊 威 沈親沐

（國防科技大學(xué)電子科學(xué)學(xué)院，湖南長沙 410073）

1 引言

恒虛警檢測（Constant False Alarm Rate，CFAR）［1］和雷達(dá)雜波分類［2］時，需要準(zhǔn)確估計出雜波數(shù)據(jù)分布模型的參數(shù)，較大的估計誤差會導(dǎo)致檢測和分類的性能下降。目前常用來描述和分析雜波幅度分布的模型主要有韋布爾分布、對數(shù)正態(tài)分布、瑞利分布和K 分布等［3］。對于以上分布，參數(shù)估計方法最準(zhǔn)確的是最大似然估計法［4］，但是最大似然估計受限于待估計參數(shù)的解析形式，當(dāng)解析形式復(fù)雜時通常求解困難或計算量大，如K 分布的最大似然估計法［5］。而矩估計方法［6］因其計算量小，克服了基于最大似然估計方法計算量大的缺點。但是，當(dāng)觀測序列的長度受限時，K 分布的高階矩估計法計算會出現(xiàn)偏差，導(dǎo)致計算出的形狀參數(shù)v<0，如二階/四階矩估計法［7］、二階/分?jǐn)?shù)階矩估計法［8］、log-Ⅲ型估計法［9］等都存在這種估計偏差。文獻(xiàn)［10］針對噪聲條件下的K 分布參數(shù)估計問題，分別提出了一種一/二/三階矩估計法和神經(jīng)網(wǎng)絡(luò)估計法，但前者估計有效率并不理想且常出現(xiàn)異常估計值，后者需要大量數(shù)據(jù)進(jìn)行訓(xùn)練。

本文針對該問題，在矩估計方法的基礎(chǔ)上，通過推導(dǎo)原點矩偏導(dǎo)與原點矩之間的關(guān)系，提出了一種新的有效的基于矩估計偏導(dǎo)的參數(shù)估計方法。這種方法對數(shù)據(jù)長度受限造成的偏差較小，能夠有效避免計算出形狀參數(shù)v<0 的情況出現(xiàn)。通過仿真和實測雜波數(shù)據(jù)處理，與常用K 分布的矩估計方法進(jìn)行比較，驗證了本方法的準(zhǔn)確性和有效性。

2 基于原點矩偏導(dǎo)的K 分布雜波參數(shù)估計方法

K分布的概率密度函數(shù)為［11］：

式中，v表示形狀參數(shù)，取值范圍為0

K分布的原點矩（k>0）為：

矩估計方法通常為了求解形狀參數(shù)v和尺度參數(shù)σ，通過選擇兩個不同的原點矩階數(shù)k，得到兩個關(guān)于v、σ和原點矩的關(guān)系式，然后聯(lián)立關(guān)系式求解出v和σ的值。這些方法是在兩個不同階數(shù)k的條件下進(jìn)行求解，通過不同階數(shù)的矩估計值代替理論值以求解形狀參數(shù)v，極可能會放大近似誤差。比如二階/四階矩估計法中，形狀參數(shù)v使用了二階矩和四階矩的近似值和，具體如式（3）所示，其中的近似計算極可能放大形狀參數(shù)v的估計誤差，甚至可能出現(xiàn)v<0的情況。

大量實測數(shù)據(jù)也表明，當(dāng)數(shù)據(jù)長度有限或者存在噪聲的影響，在二階/四階矩估計法中，當(dāng)<2時會出現(xiàn)部分樣本計算出v<0的情況，使用二階/分?jǐn)?shù)階矩估計法、log-Ⅲ型估計法等方法同樣也會出現(xiàn)v<0 的異常值。對于出現(xiàn)v<0 的異常值情況，文獻(xiàn)［12］采用了剔除異常樣本的處理辦法。

針對上述問題，通過分析原點矩表達(dá)式的特點，推導(dǎo)提出一種基于原點矩偏導(dǎo)的K 分布參數(shù)估計方法，該方法僅利用單個k階原點矩，避免了需要多個原點矩的近似值參與運算，不僅能提升形狀參數(shù)的估計有效率，而且有望提高形狀和尺度參數(shù)的估計精度。具體推導(dǎo)過程如下，令

將k看作變量，分別對G(k)和E[zk]求偏導(dǎo)，可得式（6）和式（7）：

又因為伽瑪函數(shù)有以下性質(zhì)

其中，ψ(?)是digamma函數(shù)，即普西函數(shù)（psi function），則式（6）進(jìn)一步求解，可以得到式（9）。

將式（9）代入式（7），可得

因此，聯(lián)立式（10）和式（11）可得式（12）。

分析式（2）發(fā)現(xiàn)，對E[zk]取對數(shù)同樣可以得到ln(2σ)項，即

因此，聯(lián)立式（12）和式（13）可得

當(dāng)k的值給定后，式（14）就是關(guān)于形狀參數(shù)v的方程，其是由單個k階原點矩推導(dǎo)得到，避免了不同階原點矩的求解。進(jìn)一步，式（14）還可表示為

公式（15）左側(cè)部分是待求變量v的函數(shù)，令

當(dāng)v>0時，f(v)是單調(diào)減函數(shù)，證明如下。

證明

普西函數(shù)的級數(shù)表達(dá)式

其中，γ為歐拉常數(shù)。

將普西函數(shù)的級數(shù)表達(dá)式代入式（18）并整理得

因此，式（15）通過數(shù)值計算求解方程可以得到形狀參數(shù)v的唯一解。利用式（2），不難推導(dǎo)得出尺度參數(shù)估計式，如式（20）所示。將估計的形狀參數(shù)帶入式（20）即可得尺度參數(shù)σ的估計值。

由式（15）和式（20）可以看出，本文方法形狀參數(shù)v和尺度參數(shù)σ是通過單個k階原點矩求解得到，相當(dāng)于在同一標(biāo)尺下計算形狀參數(shù)v和尺度參數(shù)σ，避免了需要兩個不同階原點矩的近似值參與運算，進(jìn)而可以避免不同階原點矩聯(lián)合計算造成估計誤差的放大，因此能提高估計的性能。

3 參數(shù)估計的有效率和估計精度分析

本部分對所提參數(shù)估計方法的有效率和估計精度進(jìn)行了實驗分析，然后對比了不同估計方法對雜波數(shù)據(jù)的尾部擬合情況，最后分析了噪聲對參數(shù)估計方法的影響。

3.1 估計有效率和估計精度分析

首先使用仿真數(shù)據(jù)驗證本文方法的估計有效率和估計精度。因為服從K 分布的雜波形狀參數(shù)通常在0.1～10之間變化［16］，本文實驗定義仿真參數(shù)v取值為［0.5∶0.5∶10］，σ取值服從均勻分布U～（0，5），利用零記憶非線性（Zero-memory Non-linear，ZMNL）方法［17］分別進(jìn)行1000次仿真，隨機(jī)生成長度256 點的K 分布數(shù)據(jù)（長度受限）。采用常用的矩估計方法，如二階/分?jǐn)?shù)階矩估計法、log-Ⅲ型估計法、二階/四階矩估計法、log（z）的一二階矩估計法［13］、基于zrlog（z）期望估計法［18］、基于zlog（z）期望估計法［19］和本文方法進(jìn)行參數(shù)估計，當(dāng)仿真參數(shù)v=3，原點矩階數(shù)k=0.5 時，7 種方法的估計有效率ρ和估計精度MSD的統(tǒng)計情況如表1所示。

表1 可以看出，當(dāng)仿真參數(shù)v=3 時，文獻(xiàn)［13］、文獻(xiàn)［18］和本文方法估計有效率為100%，但是估計精度比本文方法要差。表1 中MSD 的統(tǒng)計結(jié)果是經(jīng)過剔除掉無效估計樣本后進(jìn)行的統(tǒng)計。

表1 v=3時估計有效率和估計精度統(tǒng)計情況Tab.1 Statistics of estimation efficiency and estimation accuracy when v=3

對100%估計有效率的方法做1000個樣本的估計精度統(tǒng)計，表2所示。

盡管三種方法的估計有效率都為100%，但是表2表明了本文方法的估計精度更高。七種方法的無效估計次數(shù)與仿真參數(shù)v之間的關(guān)系如圖1所示。

表2 v=3時1000個樣本的估計精度統(tǒng)計情況Tab.2 Statistics of estimation accuracy of 1000 samples when v=3

圖1 可以看出，當(dāng)仿真參數(shù)v小于等于2 時，所有方法的無效估計次數(shù)都為0。隨著v的增加，二階/分?jǐn)?shù)階估計法、log-Ⅲ型估計法、二階/四階矩估計法、文獻(xiàn)［19］方法的無效估計次數(shù)也在不斷增多，文獻(xiàn)［13］、文獻(xiàn)［18］方法和本文方法的無效估計次數(shù)依然為0。

圖1 各方法估計有效率與仿真參數(shù)v之間的關(guān)系Fig.1 Relationship between estimation efficiency and simulation parameter v of each methods

仿真參數(shù)v=3時，文獻(xiàn)［18］和本文方法中都存在階數(shù)k，當(dāng)k的取值為［0.1∶0.1∶5］時，估計有效率ρ都能達(dá)到100%，但是本文方法的估計精度更高，k與估計精度MSD之間的關(guān)系如圖2所示。

本文方法中，當(dāng)k=0.2時MSD最小，為1.8955×10-3。當(dāng)k>0.9時，隨著k的增加MSD也會增加，當(dāng)k=0.7時，部分仿真樣本使用文獻(xiàn)［18］方法開始出現(xiàn)非正常運算產(chǎn)生的奇異值NaN，當(dāng)k越大出現(xiàn)奇異值的樣本越多。使用文獻(xiàn)［18］方法出現(xiàn)NaN的數(shù)量比使用本文方法的數(shù)量多，如k=5 時，1000 個樣本中使用文獻(xiàn)［18］方法出現(xiàn)奇異值69個，本文方法25個。出現(xiàn)NaN是因為隨著k的增大，部分樣本估計出的v會出現(xiàn)比較大的情況，甚至大于100，而gamma（100）是個非常大的數(shù)，導(dǎo)致容易出現(xiàn)Inf/Inf 的情況。圖2中，k>0.6時的MSD 值是剔除掉出現(xiàn)奇異值樣本的統(tǒng)計結(jié)果，k越大需要剔除的樣本越多。

圖2 仿真數(shù)據(jù)實驗的階數(shù)k與MSD之間的關(guān)系Fig.2 Relationship between order v and MSD in simulation data experiment

3.2 尾部擬合情況分析

累積分布函數(shù)（Cumulative Distribution Function，CDF）［20］能描述隨機(jī)變量的統(tǒng)計規(guī)律，可以對K 分布的尾部情況進(jìn)行描述。對仿真參數(shù)v=3，σ=1，樣本長度256 點的1000 組仿真數(shù)據(jù)使用CDF 進(jìn)一步驗證各方法對K 分布拖尾部分的擬合情況。圖3是各方法與仿真數(shù)據(jù)的尾部擬合情況。CDF 擬合結(jié)果來看，除了文獻(xiàn)［13］方法在尾部擬合效果較差外，其他方法在尾部擬合的都很好?？傮w而言，文獻(xiàn)［13］方法性能最差，其次為二階/四階估計法、文獻(xiàn)［18］方法、log-Ⅲ型估計法和文獻(xiàn)［19］方法、二階/分?jǐn)?shù)階矩估計法，本文方法性能最佳。

圖3 K分布的尾部擬合情況Fig.3 The tail fitting in the K-distribution simulated data

3.3 噪聲對模型影響分析

當(dāng)K 分布在有噪聲影響的情況下，傳統(tǒng)矩估計方法模型已不適用，需要對矩估計模型進(jìn)行修正。

這里用估計量的相對偏差|-v|/v表示尺度參數(shù)估計的精度。對仿真參數(shù)v=2，σ=2.2，雜噪比分別為CNR=30 和CNR=10，樣本長度為256 點的1000組仿真數(shù)據(jù)進(jìn)行分析，各方法的平均相對偏差如表3和表4所示。

表3 有30 dB雜噪比的平均相對偏差統(tǒng)計情況Tab.3 Statistics of average relative deviation with 30 dB CNR

表4 有10 dB雜噪比的平均相對偏差統(tǒng)計情況Tab.4 Statistics of average relative deviation with 10 dB CNR

在有噪聲的情況下，文獻(xiàn)［18］方法和本文方法的估計有效率都為100%，但是本文方法的平均相對偏差和估計方差最小。文獻(xiàn)［10］方法對數(shù)據(jù)敏感，部分樣本的估計值特別大甚至達(dá)到106級別，導(dǎo)致平均相對偏差和估計方差都很大。將文獻(xiàn)［10］方法估計出v>100 的樣本剔除，文獻(xiàn)［10］方法的估計精度有所提升，證明了文獻(xiàn)［10］方法對數(shù)據(jù)的敏感性。文獻(xiàn)［10］中的矩估計方法因為存在高階矩的原因，矩估計結(jié)果會出現(xiàn)特別大的v值，對數(shù)據(jù)要求較高?；谠c矩偏導(dǎo)的估計方法不會出現(xiàn)v<0的情況，并且偏差比傳統(tǒng)矩估計方法的估計偏差小，估計方差也是最小的。

在沒有噪聲影響的情況下，使用仿真參數(shù)v=2，樣本長度256點的1000組仿真數(shù)據(jù)進(jìn)行分析，各方法的平均相對偏差如表5所示。在理想無噪聲條件下，本文方法性能仍然優(yōu)于所比較的現(xiàn)有方法。

表5 無噪聲的平均相對偏差統(tǒng)計情況Tab.5 Statistics of average relative deviation without noise

實驗結(jié)果表明，當(dāng)K 分布雜波中含有噪聲時，本文方法的性能超過其他方法的性能，但是估計結(jié)果的平均相對偏差仍然較大，下一步需要對雜波模型進(jìn)行修正；當(dāng)K 分布雜波中不含有噪聲時，雜波的噪聲修正模型不適用且對數(shù)據(jù)敏感，本文方法的性能同樣超過其他方法的性能。以上數(shù)據(jù)分析表明了，相較于其他矩估計方法，本文方法的魯棒性相對更好。

4 IPIX海雜波實測數(shù)據(jù)分析

為進(jìn)一步驗證本文方法的有效性，實測數(shù)據(jù)采用IPIX 雷達(dá)1998 年測量的數(shù)據(jù)［21］，數(shù)據(jù)形式為復(fù)數(shù)形式，包含I 通道和Q 通道，脈沖重復(fù)頻率為1 kHz。IPIX 數(shù)據(jù)在加拿大格里姆斯比市安大略湖岸邊采集得到，包含不同日期不同時間不同氣象條件下的湖面回波。對160 組不同距離單元，長度為256 點數(shù)據(jù)取模進(jìn)行分析，分別進(jìn)行不同方法的參數(shù)估計，估計有效率ρ如表6所示。

表6 可以看出，所有方法的估計有效率都為100%，本文方法（k=0.5 階）的估計精度是最高的。文獻(xiàn)［13］方法出現(xiàn)6個樣本計算出異常均方差（v的估計值大于100），文獻(xiàn)［18］方法出現(xiàn)4 個樣本計算出異常均方差。表6 中MSD 的統(tǒng)計結(jié)果是經(jīng)過剔除掉異常均方差后進(jìn)行的統(tǒng)計。

表6 實測數(shù)據(jù)估計有效率和估計精度統(tǒng)計情況Tab.6 Statistics of estimation efficiency and estimation accuracy in measured data

實測數(shù)據(jù)中，當(dāng)k的取值為［0.1∶0.1∶5］時，文獻(xiàn)［18］和本文方法估計有效率ρ依然能達(dá)到100%，但是本文方法的估計精度更高，k與MSD 之間的關(guān)系如圖4所示。

本文方法中，當(dāng)k=0.3時MSD最小，為6.7876×10-3；當(dāng)k>0.9 時，隨著k的增加MSD 也會增加；當(dāng)k=0.7時，部分仿真樣本使用文獻(xiàn)［18］方法開始出現(xiàn)非正常運算產(chǎn)生的奇異值NaN，當(dāng)k越大出現(xiàn)奇異值的樣本越多。使用文獻(xiàn)［18］方法出現(xiàn)NaN 的數(shù)量比使用本文方法的數(shù)量多，如k=5 時，1000 個樣本中使用文獻(xiàn)［18］方法出現(xiàn)奇異值62個，本文方法23 個。圖4 是剔除掉出現(xiàn)奇異值樣本的MSD 統(tǒng)計結(jié)果。

圖4 實測數(shù)據(jù)實驗的階數(shù)k與MSD之間的關(guān)系Fig.4 Relationship between experimental order v and MSD in measured data

通過實驗分析，當(dāng)使用原點矩進(jìn)行K 分布的參數(shù)估計時，階數(shù)k的取值較小時能得到較為理想的估計結(jié)果和精度，這一結(jié)論與文獻(xiàn)［22］相同。文獻(xiàn)［22］指出，K 分布的參數(shù)估計中高階矩對數(shù)據(jù)較為敏感，應(yīng)盡量選取低階矩。二階/四階矩估計法、二階/分?jǐn)?shù)階矩估計法、log-Ⅲ型估計法、log（z）的一二階矩估計法、基于zrlog（z）期望估計法和基于zlog（z）期望估計法等方法都有高階矩（k≥1）的使用，并且以上方法為了求解方便，使用了不同的階數(shù)k進(jìn)行聯(lián)立方程求解，造成計算誤差的擴(kuò)大，本文方法的推導(dǎo)是基于單個階數(shù)k，所以計算誤差最小。

5 結(jié)論

本文提出基于原點矩偏導(dǎo)的K 分布雜波參數(shù)估計方法，避免了常規(guī)矩估計方法在處理雜波數(shù)據(jù)時會導(dǎo)致出現(xiàn)部分錯誤參數(shù)估計的問題。利用該方法對仿真和實測的K 分布數(shù)據(jù)進(jìn)行了參數(shù)估計，結(jié)果表明該方法具有100%的估計有效率和最小的均方誤差。通過實驗分析，當(dāng)使用原點矩進(jìn)行K 分布的參數(shù)估計時，階數(shù)k的取值較小時能得到較為理想的估計結(jié)果和精度。通過選取合適的階數(shù)k，可以提升觀測序列長度受限的K 分布雜波參數(shù)估計結(jié)果和精度，進(jìn)而提升檢測器和分類器的性能。