基于EGARCH過程的電磁頻譜占用狀態(tài)波動特性分析

王 磊 蘇東林 謝樹果 王國玉

(北京航空航天大學(xué)電子信息工程學(xué)院 北京 100191)

1 引言

無線電頻譜的占用(空閑)狀態(tài)是衡量電磁頻譜利用程度和電磁環(huán)境變化特性的重要指標(biāo)。根據(jù)大量電磁頻譜感知數(shù)據(jù)，建立反映頻譜占用狀態(tài)變化規(guī)律的數(shù)學(xué)模型，進(jìn)而實現(xiàn)頻譜占用度的精確擬合與預(yù)測，對于提高無線電管理中頻率分配的準(zhǔn)確性，深入掌握空間電磁環(huán)境變化規(guī)律，以及提升認(rèn)知無線電中動態(tài)頻譜接入效率具有重要價值。受電磁輻射源頻譜參數(shù)變化、電波傳播的衰落特性以及無線電監(jiān)測設(shè)備的性能差異等因素影響，頻譜占用狀態(tài)表現(xiàn)為與時間、空間和頻率等相關(guān)聯(lián)的非平穩(wěn)隨機過程，采用確定性方法難以準(zhǔn)確刻畫頻譜行為和電磁環(huán)境變化的內(nèi)在特征。

當(dāng)前國內(nèi)外對電磁頻譜狀態(tài)的建模和變化研究基于兩種隨機時序分析方法展開。一是將無線電信道隨時間變化描述為在占用和空閑兩個離散狀態(tài)之間轉(zhuǎn)換的馬爾科夫過程，并利用信道狀態(tài)轉(zhuǎn)移概率反映信道狀態(tài)在時域上的相關(guān)性[1-3]。這種方法在解決信道狀態(tài)隨時間長度增加呈指數(shù)增長的問題時遇到很大困難，而且單純提高馬爾科夫階數(shù)并不能帶來頻譜狀態(tài)預(yù)測準(zhǔn)確度的顯著改善[4]。二是基于無線電業(yè)務(wù)頻段采用頻譜占用度等統(tǒng)計量來描述頻譜狀態(tài)的時域變化，形成以時間先后順序排列的頻譜占用狀態(tài)時間序列。這種方法在不需要先驗知識的前提下，能夠根據(jù)頻譜狀態(tài)的樣本時序建立相關(guān)時間序列模型，定量刻畫頻譜狀態(tài)的動態(tài)變化規(guī)律，進(jìn)而預(yù)測和控制頻譜狀態(tài)的未來行為。基于這種研究方法，文獻(xiàn)[5-7]根據(jù) GSM900/1800頻段和電視頻段頻譜測試數(shù)據(jù)集和頻譜占用時間序列，建立描述其平穩(wěn)部分變化特性和自回歸(Auto Regressive,AR)時間序列模型，文獻(xiàn)[8,9]首先對 GSM900下行頻譜占用度數(shù)據(jù)集進(jìn)行預(yù)處理，而后根據(jù)數(shù)據(jù)的周期性特征建立自回歸移動平均(Auto Regressive Moving Average, ARMA)模型，首次反映該頻段頻譜占用狀態(tài)的非平穩(wěn)特性。然而上述研究均未對反映頻譜狀態(tài)時變特性的條件二階矩等非線性特征做進(jìn)一步討論。

本文在現(xiàn)有頻譜占用時間序列研究的基礎(chǔ)上，進(jìn)一步對頻譜占用狀態(tài)變化的時域相關(guān)特性進(jìn)行研究。首先通過對頻譜占用狀態(tài)ARMA模型的剩余殘差進(jìn)行異方差性檢驗，表明頻譜占用狀態(tài)具有隨時間變化的“波動集聚性”，即頻譜占用狀態(tài)轉(zhuǎn)換具有時變性和簇集性；其次基于指數(shù)廣義自回歸條件異方差(Exponential Generalized Auto Regressive Conditional Heteroskedasticity, EGARCH)過程構(gòu)建頻譜占用狀態(tài)變化等效模型，并通過對基于實測數(shù)據(jù)的頻譜占用度序列條件方差估計，證明該模型的擬合和預(yù)測精度均優(yōu)于ARMA模型；最后通過對模型參數(shù)進(jìn)行適應(yīng)性檢驗以及對“杠桿效應(yīng)”項的分析，表明頻譜占用狀態(tài)變化對電磁環(huán)境波動的影響具有非對稱性。研究結(jié)果證明采用EGARCH模型能夠量化反映頻譜占用狀態(tài)的復(fù)雜非線性時變過程。

2 基于EGARCH過程的時間序列建模方法

針對 ARMA不能反映時間序列方差隨時間變化的缺陷，Nelson在自回歸條件異方差 ARCH(Auto Regressive Conditional Heteroskedasticity)模型和廣義自回歸條件異方差 GARCH(Generalized ARCH)模型的基礎(chǔ)上，提出用于擬合時間序列中的方差變化特性和杠桿效應(yīng)的EGARCH模型[10]。它采用自然對數(shù)形式將ARMA時間序列模型中的剩余殘差分解為自回歸項、殘差項和杠桿效應(yīng)項，即在前t-1期的信息集It-1給定的條件下，對于平穩(wěn)隨機序列{Xt,t=0,±1, ±2,…}，EGARCH模型結(jié)構(gòu)可以用公式(1)表述。

約束條件：

式中，R和M分別為{Xt}的ARMA模型中自回歸階數(shù)和移動平均階數(shù)，φi和θi分別為自回歸和移動平均多項式系數(shù)，{εt} 是殘差序列，服從獨立同標(biāo)準(zhǔn)正態(tài)分布。P和W分別為擬合異方差函數(shù)值GARCH模型中的自回歸和移動平均階數(shù)。Gi,Ai和Lj分別為殘差序列的自回歸項、誤差項和杠桿效應(yīng)項的系數(shù)，K為常數(shù)。

EGARCH模型的建立過程分為ARCH效應(yīng)檢驗、模型定階與參數(shù)估計以及模型的適應(yīng)性檢驗 3個步驟。

(1)ARCH效應(yīng)檢驗 相比ARMA時間序列模型假定序列的方差為常數(shù)，EGARCH模型假定序列的方差隨時間的變化而變化(條件異方差性)，即存在ARCH效應(yīng)。條件異方差特性檢驗主要包括兩種方法，即拉格朗日乘子法(LM檢驗)和Ljung-Box Q檢驗。LM檢驗針對殘差時間序列ε(t)的平方序列建立L階自回歸模型為

檢驗統(tǒng)計量為拉格朗日乘子N×R2(其中N為樣本個數(shù)，R2為可決系數(shù))。

Ljung-Box Q檢驗的統(tǒng)計量為

其中ρ(k)2為殘差平方序列的k階自相關(guān)函數(shù)。在給定的顯著性水平α下，若序列中不存在L階ARCH效應(yīng)的零假設(shè)成立，則上述統(tǒng)計量均服從自由度為L的χ2分布。

(2)模型定階與參數(shù)估計 EGARCH 模型的階數(shù)和參數(shù)分別采用 AIC(Akaike Information Criterion)準(zhǔn)則和最大似然方法進(jìn)行估計。對于殘差時間序列ε(t)，條件對數(shù)似然函數(shù)以及用于確定最佳模型階數(shù)的AIC值分別表示為

其中l(wèi)gL與m分別為模型的對數(shù)似然函數(shù)與獨立參數(shù)。對AIC求最小值可以確定EGARCH模型的階數(shù)，求lgL的最大值可以得到模型的參數(shù)估計。

(3)模型的適應(yīng)性檢驗 利用估計出的EGARCH模型可以得到條件方差序列，進(jìn)而得到EGARCH過程的殘差序列ηt=εt/σt。殘差序列平方ηt2的Ljung-Box Q統(tǒng)計量為

通過對殘差序列進(jìn)行相關(guān)性檢驗可以確定模型參數(shù)的適應(yīng)性，對殘差序列平方進(jìn)行上述檢驗可以確認(rèn)殘差序列中是否還存在條件異方差性，從而判定EGARCH模型對頻譜占用度時間序列的適用性。

3 電磁頻譜監(jiān)測和數(shù)據(jù)預(yù)處理

本論文監(jiān)測數(shù)據(jù)來自對北京航空航天大學(xué)新主樓(北緯 39°54′ 15′′ ，東經(jīng) 116°24′ 27′′)周邊電磁環(huán)境連續(xù)7天的實際監(jiān)測(2011年10月23日2:30至2011年10月29日2:30)，監(jiān)測設(shè)備由Agilent N9340B頻譜分析儀、數(shù)據(jù)采集軟件和CS-AOS30-3000V有源全向天線組成，監(jiān)測頻段為30-3000 MHz，頻率掃描分辨率為100 kHz，測量時間樣本個數(shù)共1482個。

設(shè)M(fi,tj)(i=1,2,…,n;j=1,2,…,m)表示監(jiān)測設(shè)備在fi頻點，tj時刻測得的無線電信號場強值(dBμ V/m)，用“1”表示該信號被占用，“0”表示該信號空閑，即

其中η為頻譜占用統(tǒng)計門限值。

定義頻譜占用度SOR(Spectrum Occupancy Rate)為頻譜監(jiān)測時段內(nèi)特定無線電業(yè)務(wù)頻段占用信道個數(shù)與該頻段包括的總信道個數(shù)之比。即

其中sk(k=1,2,…,K)表示第K個無線電業(yè)務(wù)，N表示該業(yè)務(wù)所包含的信道個數(shù)。SOR隨時間先后變化取值構(gòu)成頻譜占用時間序列，且SOR∈[0,1]。

為進(jìn)一步描述頻譜占用度隨時間變化的波動特性，在對頻譜占用度時間序列去周期性的基礎(chǔ)上，再做一階差分處理，來反映頻譜占用度隨時間變化的“時變增量”，即定義頻段占用度時變序列

其中T表示頻譜占用度時間序列的周期長度，j表示差分階數(shù)。顯然，若SOBT值越大，則表示某頻段頻譜隨時間變化強度越大，因而SOBT序列可反映以無線電業(yè)務(wù)頻段為基準(zhǔn)的頻譜狀態(tài)和電磁環(huán)境波動程度。

根據(jù)實際監(jiān)測數(shù)據(jù)以及式(6)和式(7)處理方法，得到調(diào)頻廣播頻段(88-108 MHz)，電視 7-12頻道(167-223 MHz)，移動業(yè)務(wù)頻段(406-470 MHz)和GSM900上行頻段(885-915 MHz)的頻譜占用度時間序列如圖1所示。

圖1 4個無線電業(yè)務(wù)頻段頻譜占用度時間序列

4 EGARCH建模

傳統(tǒng)上對頻譜占用度時間序列建模采用的ARMA模型擬合，未考慮檢驗多項式的剩余殘差中是否還存在序列相關(guān)信息。本文通過對ARMA模型的殘差序列進(jìn)行異方差檢驗，構(gòu)建基于 EGARCH過程的頻譜占用度時間序列模型。主要流程包括序列平穩(wěn)化處理、條件異方差檢驗、EGARCH模型定階與參數(shù)估計以及條件方差預(yù)測等步驟，如圖2所示。

4.1 頻譜占用度時間序列平穩(wěn)化處理

由圖1可知，各無線電業(yè)務(wù)頻段頻譜占用度時間序列在監(jiān)測時段內(nèi)以天為單位呈現(xiàn)明顯的周期性特征，故先對時間序列進(jìn)行差分處理，差分階數(shù)為一天的采樣數(shù)，得到調(diào)整后的頻譜占用度周期性調(diào)整時間序列，如圖3所示。由該組時間序列可以看出，各業(yè)務(wù)頻段頻譜占用度時間序列的周期特征得到明顯消除，序列呈現(xiàn)平穩(wěn)變化特征，但仍具有強烈的時變特征。通過對序列進(jìn)行相關(guān)性分析，發(fā)現(xiàn)其具有較長階的自相關(guān)性，為消除相關(guān)性影響，降低ARMA模型階數(shù)，再次對該序列進(jìn)行一階差分，得到SOBT序列，如圖4所示(以GSM900上行業(yè)務(wù)為例)。由該圖可以看出，SOBT序列存在較低階數(shù)的自相關(guān)和偏相關(guān)性，可以對該組序列構(gòu)建ARMA模型。

圖2 頻譜占用度時間序列EGARCH建模流程

圖3 4個業(yè)務(wù)頻段頻譜占用度周期性調(diào)整時間序列

4.2 頻譜占用度更新時間序列條件異方差性檢驗

通過對各業(yè)務(wù)頻段頻譜占用度時間序列進(jìn)行ARMA建模，發(fā)現(xiàn)殘差序列仍然存在一定的相關(guān)性，綜合圖3和圖4可以得出SOBT序列存在如下特征：

(1)各業(yè)務(wù)頻段SOBT序列的方差隨時間變化較為劇烈，出現(xiàn)“波動集聚”特征，即方差在一定時段內(nèi)比較小，而在另一時段內(nèi)比較大。

(2)SOBT序列不符合正態(tài)分布，且樣本自相關(guān)性較小(接近于隨機分布)，樣本平方存在顯著的相關(guān)性。

上述特征直觀表明SOBT序列可能存在條件異方差性。進(jìn)一步利用式(3)和式(4)對各業(yè)務(wù)頻段時間序列進(jìn)行24階LM 檢驗和Ljung-Box Q檢驗，結(jié)果如表1所示，由于各階LM統(tǒng)計量和Q統(tǒng)計量均大于 0.05顯著水平的臨界值，即拒絕序列不存在ARCH效應(yīng)的零假設(shè)，因而可確認(rèn)各無線電業(yè)務(wù)頻段SOBT序列存在條件異方差特征。

表1 4個業(yè)務(wù)頻段ARCH效應(yīng)檢驗結(jié)果(α=0.05)

4.3 EGARCH模型定階與參數(shù)估計

EGARCH過程的階數(shù)由SOBT序列ARMA模型剩余殘差序列存在的條件異方差效應(yīng)決定。如果檢驗表明時間序列存在較高階的條件異方差效應(yīng)，則可以考慮擬合低階 EGARCH模型。同時可以結(jié)合所擬合的各階模型的對數(shù)似然函數(shù)與AIC值進(jìn)行模型定階與參數(shù)估計。

根據(jù)各業(yè)務(wù)頻段SOBT序列的自相關(guān)和偏相關(guān)系數(shù)，首先建立調(diào)頻廣播和電視頻道、移動業(yè)務(wù)(406-470 MHz)和GSM900上行業(yè)務(wù)的ARMA(2,1)模型，而后根據(jù)式(5)和式(6)建立剩余殘差序列的一階條件異方差模型，其中各業(yè)務(wù)頻段 EGARCH模型參數(shù)估計結(jié)果見表2。根據(jù)表2并結(jié)合式(1)即可確定各業(yè)務(wù)頻譜占用度時變序列的 EGARCH模型具體表達(dá)式。由表2知，各模型L1值均不等于零，表明各業(yè)務(wù)頻譜占用波動均存在不同程度的“杠桿效應(yīng)”，具體來講，對于調(diào)頻廣播頻段、電視 7-12頻道和移動業(yè)務(wù)頻段，由于其L1值均大于零，故這些頻段頻譜占用度增加更易引起電磁環(huán)境波動；而對于GSM900上行頻段，由于其L1值小于零，故該頻段頻譜占用度減小更易引起電磁環(huán)境波動。

4.4 EGARCH模型適應(yīng)性檢驗

圖4 GSM900上行頻段SOBT正態(tài)檢驗結(jié)果和相關(guān)性系數(shù)

表2 頻譜占用度EGARCH模型參數(shù)估計結(jié)果

利用估計的各業(yè)務(wù)頻段頻譜占用度條件標(biāo)準(zhǔn)差序列，可以得到相應(yīng)的殘差序列。通過計算調(diào)頻廣播和電視頻道、移動業(yè)務(wù)(406-470 MHz)和GSM900上行業(yè)務(wù)頻譜占用度 EGARCH模型剩余殘差的自相關(guān)系數(shù)和殘差平方自相關(guān)系數(shù)，兩者均沒有顯著的自相關(guān)性，說明EGARCH(1,1)模型殘差序列是相互獨立的，已不存在 ARCH效應(yīng)。同時，按式(7)對各業(yè)務(wù)頻段殘差序列平方進(jìn)行Ljung-Box Q統(tǒng)計量檢驗，接受不存在ARCH效應(yīng)的零假設(shè)，即各業(yè)務(wù)頻段SOBT序列均通過 EGARCH模型適應(yīng)性檢驗。

4.5 頻譜占用度時間序列條件方差估計

根據(jù)式(1)模型表達(dá)式和表2所列參數(shù)估計值，可以得到各無線電業(yè)務(wù)頻段頻譜占用度時變序列的條件方差表達(dá)式，如對于調(diào)頻廣播業(yè)務(wù)該式為

由式(11)和式(12)通過遞推可得到該序列的條件方差估計值。通過對各業(yè)務(wù)頻段頻譜占用度時間序列分別建立ARMA模型和EGARCH模型可以進(jìn)行相應(yīng)頻譜占用度時間序列預(yù)測。圖5為對4個業(yè)務(wù)頻段頻譜占用度分別建立 ARMA(2, 1)模型和EGARCH(1, 1)模型得到的條件方差估計結(jié)果。

通過對4個無線電業(yè)務(wù)頻段頻譜占用時間序列條件方差大小進(jìn)行對比分析得出，采用 EGARCH模型得到的序列方差隨時間發(fā)生不均勻變化，而采用ARMA模型得到的序列方差為常數(shù)，說明無線電業(yè)務(wù)頻段頻譜占用度具有顯著的波動集聚特征，即大的波動后面傾向于跟著一個大的波動，小的波動后面傾向于跟著一個小的波動，即頻譜占用度狀態(tài)的前后之間存在著某種相關(guān)關(guān)系；通過對無線電業(yè)務(wù)頻段頻譜占用度的波動幅度與占用度絕對值大小以及頻率高低等因素觀察，88-108 MHz頻段的頻譜占用度波動幅度最大，其次為167-223 MHz頻段和885-915 MHz頻段，406-470 MHz頻段頻譜占用度幅度變化最小。這說明，頻譜占用度越高，占用度的波動幅度越大，頻率越低，占用度的波動幅度一般越小。

圖5 4個業(yè)務(wù)頻段頻譜占用度時間序列條件方差估計(周變化)

圖6 4個業(yè)務(wù)頻段頻譜占用度時間序列條件方差估計(日變化)

為了進(jìn)一步分析無線電業(yè)務(wù)頻譜占用度在一天中的波動情況，我們在上述以周為時段進(jìn)行研究的基礎(chǔ)上，分析了88-108 MHz頻段、167-223 MHz頻段、406-470 MHz 以及885-915 MHz頻段的頻譜占用度在一天中以小時為單位的波動狀態(tài)分布，如圖6所示。由圖可以看出，4個無線電業(yè)務(wù)頻譜占用度在一天中的波動狀況基本相似，每天會分別在8:30-10:30時段和13:30-14:30時段出現(xiàn)兩個頻譜占用度變化峰值，這兩個時段正好與監(jiān)測地域(校園)人為用頻活動最為密集的時段相吻合。同時，GSM900 上行業(yè)務(wù)頻段(885-915 MHz)頻譜占用度從8:00-24:00時段均保持高位波動，而在其余時段波動幅度相對平坦，這種變化樣式與人們使用移動通信頻率和作息時間正好吻合。上述分析結(jié)論說明EGARCH模型能夠較好地擬合和解釋頻譜占用度的波動現(xiàn)象。

5 結(jié)論

基于電磁環(huán)境實測數(shù)據(jù)，對無線電業(yè)務(wù)頻譜占用時間序列進(jìn)行建模與分析，是量化描述無線電頻譜狀態(tài)波動特性、揭示電磁環(huán)境復(fù)雜變化機理的可行方法。通過本文的研究工作，得出以下幾點結(jié)論：

(1)無線電業(yè)務(wù)頻譜占用時間序列方差具有顯著的條件異方差特性，采用 EGARCH模型相較傳統(tǒng)的 ARMA模型能夠客觀反映頻譜占用序列條件二階矩的時變特性，且模型擬合和預(yù)測精度更高。

(2)無線電業(yè)務(wù)頻譜占用狀態(tài)在時域上表現(xiàn)出顯著的周期性、共動性和簇集性變化特征，即頻譜占用狀態(tài)在時域上表現(xiàn)出以天為周期交替變化，各業(yè)務(wù)頻段頻譜占用度在某些時段變化幅度產(chǎn)生強烈波動，在其余時段波動幅度變化較小，且電磁環(huán)境波動時段與人為用頻活動具有強關(guān)聯(lián)性。

(3)無線電業(yè)務(wù)頻譜占用狀態(tài)在頻域上表現(xiàn)出顯著的波動差異和“杠桿效應(yīng)”，即頻譜占用狀態(tài)變化幅度隨著頻率的升高呈現(xiàn)減小趨勢，同時頻譜占用狀態(tài)的變化對整體電磁環(huán)境波動的影響具有非對稱性。在某些頻段，如調(diào)頻廣播頻段、電視7-12頻道以及移動業(yè)務(wù)頻段，頻譜占用度增加會帶來較大的電磁環(huán)境波動，而在其它頻段，如GSM900上行頻段，頻譜占用度減小會帶來較大的電磁環(huán)境波動。

上述結(jié)論主要通過對日常電磁環(huán)境進(jìn)行監(jiān)測和數(shù)據(jù)分析得出的，遇有重大活動或特殊用頻情況下的電磁環(huán)境異常波動規(guī)律還需進(jìn)一步研究。

[1]López-Benítez M and Casadevall F. Discrete-time spectrum occupancy model based on Markov chain and duty cycle models[C]. IEEE International Symposium on Dynamic Spectrum Access Networks, Aachen, Germany, 2011: 90-99.

[2]Li Yang, Dong Yu-ning, Zhang Hui,et al.. Spectrum usage prediction based on high-order Markov model for cognitive radio networks[C]. 10th IEEE International Conference on Computer and Information Technology, Bradford, UK, 2010:2784-2788.

[3]Datla D, Wyglinski A M, and Minden G J. A spectrum surveying framework for dynamic spectrum access networks[J].IEEE Transactions on Vehicular Technology,2009, 58(8): 4158-4168.

[4]Song Cheng-qi, Chen Da-wei, and Zhang Qian. Understand the predictability of wireless spectrum: a large-scale empirical study[C]. IEEE International Conference on Communications,Cape Town, South Africa, 2010: 1-5.

[5]Yin Si-xing, Chen Da-wei, Zhang Qian,et al.. Mining spectrum usage data: a large-scale spectrum measurement study[J].IEEE Transactions on Mobile Computing, 2011,11(6): 1033-1046.

[6]Gorcin A, Celebi H, Khalid A Q,et al.. An autoregressive approach for spectrum occupancy modeling and prediction based on synchronous measurements[C]. IEEE 22th International Symposium on Personal, Indoor and Mobile Radio Commission, Toronto, Canada, 2011: 705-709.

[7]Pagadarai S and Wyglinski A M. A quantitative assessment of wireless spectrum measurements for dynamic spectrum access[C]. Proceedings of the 4th International Conference on CROW NCOM, Hanover, Germany, 2009: 1-5.

[8]Yarkan S and Arslan H. Binary time series approach to spectrum prediction for cognitive radio[C]. IEEE Vehicular Technology Conference, Baltimore, MD, USA, 2007:1563-1567.

[9]Wang Zhe and Salous S. Time series ARIMA model of spectrum occupancy for cognitive radio[C]. IET Seminar on Cognitive Radio and Software Defined Radios, 2008: 1-4.

[10]Wei W S. Time Series Analysis: Univariate and Multivariate Methods[M]. Second Edition, New Jersey, USA, Pearson-Addison Wesley, 2006: 345-359.