1、第五章目錄第五章極大熵譜估計15.1 譜熵和極大熵準則11問題的提出12高斯過程的熵和熵率13功率譜和熵率的關系35.2 極大熵準則的譜估計65.3 極大熵譜估計的伯格算法95.4 極大熵譜估計的LSLUD算法16第五章 極大熵譜估計1967年伯格(JPBurg)剛一發表：極大熵譜分析”的方法就在工程和科技界產生很大影響，成為相當流行的功率譜密度估計方法。伯格在譜估計準則的提出和具體算法上有所創新，由此演變出來的算法有很多種，被統稱為“現代譜分析”。5.1譜熵和極大熵準則1問題的提出從19世紀未舒斯特(Schuster)在利用富氏級數分析信號隱含的周期特性時提出了“周期圖”，到1985年由伯來

2、克曼和杜奇提出了譜估計的“間接法”和1965年FFT算法提出后流行的“直接法”，它們本質上都是把原序列經過開窗截取處理來獲得對序列譜密度的估計。不論對數據加窗還是對自相關函數加窗，其目的都在于使譜估計的方差減小，然而加窗不可避免地產生頻域“泄漏”，使功率譜失真，盡管在窗函數形式的選擇和處理方法上做了很多分析研究，使得以周期圖為基礎的方法達到相當成熟和實用的程度，但是任何抑制旁瓣的方法都是以損失譜分辨力為代價的，這個難題在數據量少的情況下更為突出。問題的實質是：在周期圖估計中，我們對數據或是它的相關函數所做的加窗處理，等于是假定在窗口外數據(或自相關)為零，而窗口內的部分則加上某種形式的修正。這

3、些人為措施使來自觀察的信息受到了一定程度的歪曲。伯格提出的新概念是；和估計的功率譜相對應的自相關和由觀察數據算得的自相關一致，同時對已有的區段之外的自相關值采用外推的辦法求取，而不是一概假定為零，外推的原則是使相應的序列在未知點上取值的可能性具有最大的不確定性，亦即不對結果人為地強添任何增加的信息。數學家申農最早提出“熵”的概念，在統計學中用它作為各種隨機試驗的不肯定性程度的度量。在熱力學和信息論中，“熵”都有其具體的物理背景和應用。后面介紹將會看到，滿足熵極大的譜估計是自回歸模型的譜。1971年凡登包士（Van Den Bos）證明，一維極大熵譜估計和自回歸譜的最小二乘估計是等效的。盡管如此

4、，伯格關于熵譜估計的概念和他對自回歸參數的遞推算法卻獨樹一幟，隨后還有人提出了各種改進算法，但要注意把極大熵概念本身同等法區別開來。2高斯過程的熵和熵率假定我們研究的隨機試驗a只有有限個不相容的結果，它們相應的概率為，且滿足，簡單描述如下：申農找到并證明了可以用這個量來度量的不肯定性的程度：或簡寫成：稱為試驗的熵當隨機變量的可能取值是連續的，則H定義式中的和式用積分代替(5-1-1)其中為隨機變量，對數可以取10或取e為底，在比較熵的大小時并沒有影響，下面為計算方便均以自然對數ln來定義，如x為正態隨機變量，則有 (5-1-2)進一步，如果討論的是時間序列的實現則這一過程的熵用下面N維積分表示

5、： (5-1-3)其中是聯合概率密度函數，若時間序列是高斯的，則 (5-1-4)其中為自協方差陣 (5-1-5)它的i行j列元素為的均值，表均值向量。將式（5-1-4）代入式（5-1-3）求過程x的熵 (5-1-6)式（5-1-6）就是長度為N的正態時間序列的熵。若有正態白噪聲（方差為），則可求得其熵為 (5-1-7)由于H隨N增長而發散，定義熵率h為 (5-1-8)故白噪聲過程的熵率為 (5-1-9)3功率譜和熵率的關系下面給出功率譜和熵率間的一些重要性質和關系。（1）如果隨機向量是隨機向量，則由于 (5-1-10)其中A是N×N非奇異矩陣，X的聯合概率密度為，則由于 (5-1-1

6、1)可得 (5-1-12)（2）若是一個穩定的因果系統的輸入，該系統的傳遞函數為G(B)（這在單位園內無極點），系統單位脈沖響應為。設在時開始輸入，因而系統輸出是平穩的，以 (5-1-13)其中 (5-1-14)證明：若在t=0時開始輸入，則系統輸出為 (5-1-15)。式（5-1-15）是隨機變量通過線性變換Ax成為隨機變量，這里根據式(5-1-12)得除以（t+1）并令則現在只要證明等于式（5-1-13）中的第二項積分就夠了。由于，情況下有這里的線積分是沿單位園進行的，因故式（5-1-13）中的第二項積分等于，所以需要證明由于G(B)在單位園內是解析的，所以上式中的積分路線可以任意小，當。

7、故上式右邊等于，證明完成。（3）若是正態過程，其功率譜滿足 (5-1-15)則有 (5-1-16)注：這一結論是不難看出的，因為非白正態過程的功率譜密度可以看作是方差為1的白噪聲通過頻率響應模的平方等的線性系統所產生的過程的譜，因此利用式（5-1-13）和（5-1-9）就可導出式（5-1-16）。式（5-1-16）給出了過程的功率譜密度和它的熵率之間的關系式，由于右邊第一項是常數，比較的大小等價于比較第二項積分的大小，因此稱的譜熵，并以它作為推導極大熵譜估計的出發點。例如已知過程的方差為，即 (5-1-17)要導出能使為最大的功率譜。這個問題可以通過求泛函極值來解決。以表拉格倫日乘子作泛函其變

8、為達到極值的條件為,故應有代回式（5-1-17）可得，即必須為常數，因此只有當過程為白噪聲時才能使熵率達到最大，這里約束條件是方差為。5.2 極大熵準則的譜估計根據伯格所提出的概念，功紡譜密度估計的準則應當是：設表示估計的譜，則它在滿足約束條件 (5-2-1)的同時，應使譜熵達到極大，其中是的正、實、偶函數，這樣對應的R(r)自然也是r的偶函數。下面論證滿足以上要求的所應具有的形式。設已知自相關函數R(r)在內的2M1個值，以表拉格倫日乘子作泛函由得(5-2-2)這里應有(5-2-3)以保證是實的。將式（5-2-2）代入式（5-2-1）得(5-2-4)用后移算子代替變量，上式所寫成(5-2-5

9、)該積分沿B平面單位園進行，基于式（5-2-3）有(5-2-6)其中：不難看出故多項式的全部零點均在B平面單位圓外，而的全部零點均在單位圓內，兩部分零點是互為倒數分布的。將式（5-2-6）代入式（5-2-5）得構造和式由于在單位圓內沒有零點，上式被積函數除了r=0時有極點在原點外，r1時在單位圓內是解析的，根據柯西留數定理可得詳細寫，就是若令(5-2-7)則有(5-2-8)利用式（5-2-7）可將式（5-2-6）寫成(5-2-9)其中將式（5-2-9）代入式（5-2-2），并考慮到的偶函數性質，得(5-2-10)由結果可以看出，在已知頭M+1個相關函數值時，以它作為約束條件推出的極大熵譜是AR

10、模型的功率譜。一般情況下，我們直接觀察到的是過程的數據，并不知道相關函數的準確值。因此通常根據為最小來求，這將得到一組n個方程，其形式和式（5-2-8）一樣，只是用估計的自相關代替R(r)，所以當采用這種估計值時，極大熵法和最小二乘法估計的結果是相同的。5.3 極大熵譜估計的伯格算法上面已經指出，由于不確切知道頭M+1個自相關的真值一般只能用它的估計值代替，在第三章中提到的兩種估計算式為和是非負定的，方差較小，但估計偏度隨r增加而增大。R是漸近無偏的，但不能保證非負定性質。因此以上兩種算法各有不足之處。伯格采用的算法是從一階模型開始逐步增加階數的遞推算法，每步遞推都能保證相應的自相關序列是非負

11、定的，而且得到的模型也是平穩的，不僅如此，從后面介紹還可看出，由于采用了雙向（正、反）預報誤差平方和為最小，提高了數據的利用率，充分挖取數據內含的信息，因此特別有利于短數據的分析和建模。先看一階模型這里上標是用來標明它屬于一階模型，以下都用它作為遞推次數的記號，對平方和作極小化來選擇時可能會出現，從而使模型不平穩，例如序列值為1，2時由由于在均方意義下最優的模型參數只取決于序我的自相關而非序列值本身，而序列按反向的時間順序排列并不改變自相關函數（見2.1節），伯格提出以作極小化來求，上式右邊第二個括弧內的可以看作是由“預報”的誤差，不難看出，由此得到的滿足。因此這種將正反雙方預報誤差的平方和作

12、極小化的辦法在一階的情況下是可行的如果自然地推廣到二階，則有對它作極小化求然而伯格指出，這樣求得的模型參數并不總能滿足平穩性條件，但他注意到利文森（Levinson）提出的遞推算法可以做到這一點，這種算法由AR(1)推出AR（2）是按以下關系：伯格決定不由，作極小化，而是把看作是的函數，然后對極小化，為此寫出其中表示一階模型的正向預報誤差，表示相應的反向預報誤差。由可得因恒大于零，可以證明。后面將會看到這對模型的平穩性和自相關函數的非負定性都是很關鍵的。求后，AR(2)的參數是從AR(2)向AR(3)遞推的方式和前面類似，令于是其中分別為二階模型的正摻向預報誤差同樣，由求得且有這樣的梯推繼續下

13、去，到M階時的一般算式有(5-3-1)現在來看和自相關函數值的遞推式，令階數從1起和的遞推計算將利用自相關函數陣的對稱托普里茨（Toeplitz）性質，以二階向三階遞推為例，由式（5-2-8）寫出(5-3-2)如果將該矩陣所對應的議程組及變量的順序反過來，則有(5-3-3)把式（5-3-2）的相關陣放大到4×4；可得(5-3-4)其中(5-3-5)對式（5-3-3）也作類似擴大，然后將兩個4×4陣按下面方式組合(5-3-6)根據式（5-3-1）中關于自回歸參數的遞推關系可見上式左邊有關參數的矩陣乃是三階模型的參數，因此上式右邊等于，于是有(5-3-7)解之得利用式（5-3-

14、5）和式（5-3-7）得推廣到一般，可綜合如下遞推算式：(5-3-8)(5-3-9)(5-3-10)以上三式和（5-3-1）諸式組成了一套極大熵譜估計的伯格遞推算法（程序見附錄二MEBURG）。（1） 遞推所得的參數滿足平穩性條件，即的根全部在B平面單位圓以外，或者等效地說中任一都滿足。證：令并將表為，則B應當是一個半徑小于1的圓，或者說是圓心在(1,0)但不包圍原點的圓（見圖5-1）。這必等價于當f由-1/2時變到圖5-1 的根滿足平穩條件時的曲線1/2時的幅角增量為零。當全部的模均小于1時，總的幅角增量也是零，或曲線不包圍原點。而今(5-3-11)這里模故當是不包圍原點的，即式（5-3-1

15、1）右邊·部分的零點都在B平面單位圓外，如果前一步遞推得到的已滿足平穩條件，則也將滿足平穩條件。由于從一階開始遞推時已有，且以后每步遞推均有，因此每步遞推所得的參數必然均能滿足平穩條件。（2）遞推所得的自相關序列滿足非負定條件。證：由于，根據式（5-3-9）必有。再以表示由，構成的相關函數陣，則式（5-2-8）可寫成由克萊姆法則知由于按歸納法可知每次遞推所得的，因此上述遞推算法得出的序列構成正定列。關于由極大熵譜獲得的模型階數問題，由式（5-2-10）可見，其階數M是已給自相關估值的最大遲后，當數據個數為N時最大可能的遲后值為N-1，這可能并非是過程的真正階數。而另一方面，如果序列本

16、身是無限階的AR模型（如ARMA模型的等效），需要很高的階數才能逼遠真正的過程，這時已給相關的最大遲后所定出的階數又可能太小。當然，M愈大，用估計的精度也愈低，所以取很大的階數未必就好。鑒于極大熵是AR譜，我們可以利用諸如FPE、AIC、BIC等定階準則進行檢驗和判定。下面的例子是一組由產生的20個數據其中是白噪聲，它的標準差約為正弦振幅的5%，一個的紀錄為0.1410，1.0509，1.7826。2.6804，3.0536, 2.9605, 2.7524, 2.1767, 1.6413, 1.371, 0.1217, -0.9359, -1.8501, -2.5495, -2.5454, -

17、2.9358, -3.0448, -2.2961, -1.7726, -0.9091（見圖5-2），利用MEBURG程序計算可得最佳階數最小殘差方差為0.0304。圖5-3為極大熵譜曲線，其峰值出現在處。圖5-2 20點數據圖圖5-3 20點數據的極大熵譜AR（6）這里需要指出的是，用自回歸模型擬合只適用于有譜密度的序列，對正弦信號而言其譜密度在給定頻率處為無窮。這個例子只是說明根據短的序列樣本以極大熵譜估計譜線的位置，即正弦振蕩的頻率。 極大熵譜估計的伯格算法程序見附錄二MEBURG。5.4極大熵譜估計的LSLUD算法 伯格在極大熵基礎上提出的算法，由于采用了較合理的外推和正反向誤差平方和的

18、極小化，比傳統的方法提高了分辨力。或者說在同樣的分辨力下只需要較少的數據量。在伯格方法引起人們廣泛重視和應用的同時，也發現其不足之處，這就是“譜峰偏移”和“譜線分裂”，前者是指峰值頻率估值和真值之間的偏離度，后者是指本來只有一個譜峰，但在估計譜中卻出現兩個或多個相距很近的譜峰。這類現象在耶爾一瓦克爾的AR譜估計中也是存在的(凱伊(Kay)和馬波(Marple)已曾指出過)，而伯格算法依然存在這兩個問題。泛杰爾(Fougere)首先指出；當數據中信噪比高以及所取階數較高的情況下，伯格算法容易產生譜線分裂，在用周期為了的正弦信號疊加白噪聲作樣本進行分析中發現，當數據長度為T4的奇數倍，以及正弦的初

19、始相位為的奇數倍時也容易出現分裂。譜峰偏移也和初始相位有關。而數據量加大，譜線分裂和譜峰偏移量都會減弱。為了解決上述問題，人們已經并且還在做出努力，通常在研究中都采用正反向預報誤差，而實踐結果表明，AR模型的最小二乘(LS)解的頻率偏移小，在短數據下，不受托普里茨(Toeplitz)矩陣形式限制的LS方法可以得到更好的結果。不過LS方法計算量大模，型結果可能非平穩等問題則需要妥善解決。這方面近幾年來已取得一些肯定的結果，本節介紹的LS-LUD算法就是在最小二乘方法基礎上采用上下三角陣分解(LUD)的算法求解AR譜。我們知道。若給定時間序列為，當用n個既往觀察數據作一步正向預報時。其預報誤差為(

20、5-4-1)其中(5-4-2)根據LS準則，。類似地，當用n個后來觀察數據作反向預報時。其預報誤差記作(5-4-3)其中(5-4-4)而可由的準則確定利用矩陣方程表示時。式（5-4-1）可寫成(5-4-5)其中式（5-4-5）是含有n個未知參數的N-n個方程組。殘差平方和記作根據可得n個方程式正規方程：(5-4-6)或(5-4-7)其中 類似地，由反向預報誤差可以得出正規方程(5-4-8)或(5-4-9)其中如果同時考慮正反向預報誤差，并對兩種預報采用同樣的自回歸系數，則可寫出含有這n個參數的2(N-n)個方程式(5-4-10)其中而實現的正規方程為(5-4-11)其中實際確定參數的過程包含兩

21、大步驟，一是正規方程的列寫，二是正規方程的求解。一般說，對于含n個未知參數的M個超定方程組，其正規議程的建立需要次運算，用巧列斯基（Cholesky）方法求解需要次運算（一次運算是指一個乘法或除法加上一個加法或減法，在比較運算量時只考慮M和n的最高次冪）。對應于式（5-4-10），M=2(N-n)，但在求解AR模型參數一特定情況下，由n階建立n+1階正規方程可以采用遞推的方法，下面將表明：對于單向（正向或反向）預報來說，建立正規議程的計算量為，而對于雙向預報的式（5-4-10），其計算量為。由于階數往往要從1開始搜索，在合適階數未知的情況下，從1到某個階數n2逐個建立正規方程的總計算量為：以正

22、向預報中由n=2推出n=3的情況為例，考慮到正規方程矩陣的對稱性，只寫它的下三角部分注意到中由構成的2×2子矩陣的各元素等于中的元素減去元素所定義的內積中的第一項，即最后一行中除了第一個元素外，都可由的最后一行推得，即此外，除了最后一個元素外，亦可由推得可見，除了要算一個內積中的其余元素只要一次運算便可求得。一般情況下，由需要多于一次的運算。關于三種方式下正規方程的建立過程用編程格式列出如下：（1）正向預報方式算法（2）反向預報方式算法（3）雙向預報方式算法關于n×n正規方程的求解，第一步是將R（它是對稱正定陣）分解成（L為下三角陣），所需運算量為，第二步是解，然后由，需要的運算量為，兩步共需的運算，這樣從到某一階數n獨立求解所需的總計算量為。由于計算中采用了三角陣的分解（巧列斯基分解），故簡稱LUD算法（Lower and Upper Triangular Matrix Decomposition）。總的說來，在情況下，全部計算中占主要的計算是建立正規方程，雖然LS-LUD算法的計算量比伯格算法大，但相差不是數量級的。