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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

16、身是無限階的AR模型（如ARMA模型的等效），需要很高的階數(shù)才能逼遠(yuǎn)真正的過程，這時(shí)已給相關(guān)的最大遲后所定出的階數(shù)又可能太小。當(dāng)然，M愈大，用估計(jì)的精度也愈低，所以取很大的階數(shù)未必就好。鑒于極大熵是AR譜，我們可以利用諸如FPE、AIC、BIC等定階準(zhǔn)則進(jìn)行檢驗(yàn)和判定。下面的例子是一組由產(chǎn)生的20個(gè)數(shù)據(jù)其中是白噪聲，它的標(biāo)準(zhǔn)差約為正弦振幅的5%，一個(gè)的紀(jì)錄為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程序計(jì)算可得最佳階數(shù)最小殘差方差為0.0304。圖5-3為極大熵譜曲線，其峰值出現(xiàn)在處。圖5-2 20點(diǎn)數(shù)據(jù)圖圖5-3 20點(diǎn)數(shù)據(jù)的極大熵譜AR（6）這里需要指出的是，用自回歸模型擬合只適用于有譜密度的序列，對(duì)正弦信號(hào)而言其譜密度在給定頻率處為無窮。這個(gè)例子只是說明根據(jù)短的序列樣本以極大熵譜估計(jì)譜線的位置，即正弦振蕩的頻率。 極大熵譜估計(jì)的伯格算法程序見附錄二MEBURG。5.4極大熵譜估計(jì)的LSLUD算法 伯格在極大熵基礎(chǔ)上提出的算法，由于采用了較合理的外推和正反向誤差平方和的

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

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

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

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

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