基于閃爍譜擬合的不均勻體參數最大似然估計

陳濤 劉鈍 於曉 郝藝 郭珊 朱慶林

（1.中國人民解放軍91198 部隊, 北海 536100；2.中國電波傳播研究所, 青島 266107）

0 引 言

電離層閃爍是影響GNSS 及其應用的重要空間環境因素，對其進行監測預測是應對閃爍影響的有效措施.閃爍影響不僅與造成閃爍的小尺度電離層電子密度不均勻體特征有關(包括擾動強度、譜指數、內外尺度等)，還與用戶的觀測幾何有關[1-2].直接利用用戶測量結果如閃爍指數S4、電子總含量(total electron content, TEC)變化率指數(rate of TEC index,ROTI)等進行閃爍預測建模，會使模型變得復雜且預測性能有限[3-4].因此，已有的閃爍模型均采用對不均勻體特征進行建模，并通過信號傳播模型獲得接收機處信號變化特征.如WBMOD 模型中對不均勻體強度、譜指數、漂移速度等參數進行建模；全球電離層閃爍模型(Global Ionospheric Scintillation Model，GISM)通過對電子密度起伏變化建模[5-6]，并通過相位屏傳播模型模擬產生閃爍信號[7].

利用閃爍觀測數據對擾動強度、譜指數、漂移速度等不均勻體參數進行有效估計是建立閃爍預測模型的基礎.通常利用衛星信號對閃爍進行觀測，隨著GNSS 應用的發展，利用GNSS 信號進行電離層閃爍監測已成為一種重要手段.對衛星信號進行高頻采樣觀測并做譜分析，可以獲得閃爍信號頻譜.理論分析表明，造成閃爍的電離層不均勻體空間譜特征，與接收的衛星閃爍信號的頻譜特征具有相關性[8-9].因此可以通過對衛星閃爍信號頻譜特征的深入分析，獲得電離層不均勻體參數的估計值.

利用譜擬合方法對空間環境參數進行估計是行星際閃爍、X 射線等天文學研究中的一種常用方法[10-12].但GNSS 信號電離層閃爍測量須要建立合適的衛星閃爍信號功率譜模型，才能利用衛星信號譜分析結果進行閃爍參數估計.

閃爍理論表明，弱散射情況下衛星信號的相位譜和功率譜在高頻端表現出冪率譜下降特征.C.Rino基于相位屏理論，建立了弱散射情況下的信號譜模型，給出衛星信號譜的表達式，并指出在高頻端存在冪率譜近似Tscin/f p， 其中Tscin為1 Hz 處的相位譜強度，p為譜指數[13].基于此模型，一些學者通過對GPS 衛星閃爍信號功率譜高頻端一定頻率范圍內的譜系數進行擬合，實現對譜強度和譜指數的估計[9,14].但該方法在實際應用中受到很多限制.為保證信號功率譜在高頻端呈現出冪率譜變化特征，須選擇一定強度的閃爍信號，并須要合理選擇擬合的頻率范圍.實際觀測表明，衛星信號閃爍功率譜在一定強度閃爍下會出現擴展或增強[15]，在信號高頻端很難看到明顯的冪率譜變化特征.同時，不均勻體的漂移速度變化、衛星自身運動，以及接收機觀測噪聲影響，難以確定一個固定的頻率范圍進行擬合[9,14].為此，C.S.Carrano 等人[15-16]進一步利用C.Rino 建立的強擾動下閃爍信號功率譜模型[17]進行直接擬合，對譜強度、譜指數等參數進行最優估計.該方法的優勢是選擇最小截止頻率和接收機高頻觀測噪聲之間一段相對確定的頻譜范圍進行擬合，提高了譜擬合進行參數估計的可用性，但該方法缺乏對估計結果的合理評價.后面分析可以看出，準確的譜模型能夠保證參數估計結果符合特定分布，從而實現對結果置信度的估計.近年來C.Rino 等人針對GNSS 應用提出了雙指數閃爍譜模型，并通過數值進行模擬發現，雙指數譜模型可以表現出不同閃爍情況下的信號功率譜變化特征[15,18]；進一步又利用雙指數譜模型對天基測量數據進行了分析，驗證了模型的有效性[19].

已有的長期閃爍觀測數據主要為單站測量的衛星信標數據，如GPS 信號、通信衛星信號等，因此本文利用閃爍功率譜擬合方法對閃爍相關參數進行最優估計，并利用仿真數據和實測數據對該方法的性能進行分析評估，統計分析結果表明了利用譜擬合方法估計閃爍參數的有效性.

1 電離層閃爍譜模型

C.Rino 等人提出的雙指數電離層閃爍譜模型已廣泛應用到GNSS 閃爍影響分析中[15,20-21]，雙指數譜模型中定義了如下不均勻體相位譜密度函數模型[18]：

利用菲涅爾半徑 ρF對空間波數進行歸一化處理，可以獲得更為統一的相位譜密度形式[18]

式中：譜強度參數為

利用歸一化的相位譜密度函數，接收信號強度的功率譜可以表示為[15,17-18]

式中： η =r/ρF為 利用菲涅爾半徑 ρF對 空間間隔r進行歸一化的結果； γ為結構交互函數，定義為[18]

χ為積分變量.

實際觀測中獲得的是閃爍影響下信號起伏變化的時間序列，而時間序列變化由衛星信號穿越電離層不均勻體的相對運動引起.因此，利用有效掃描速度veff可以進行空域譜和時域譜參數之間的轉換，即μ =2πf/fF， 其中菲涅爾頻率fF=veff/ρF.進而獲得閃爍影響下信號強度變化的時域譜形式為[21-22]

由閃爍理論進一步可知，在已知信號強度功率譜情況下，存在下述關系式[15,18]：

2 閃爍譜參數的最大似然估計(MLE)

利用閃爍影響下的信號強度時域譜模型，采用最小二乘估計和最大似然估計(maximum likelihood estimation, MLE)對實測信號強度譜進行擬合[16,22]，實現對信號強度譜模型中不均勻體參數的估計.在已知信號概率分布情況下，利用MLE 可以獲得對信號參數的無偏估計.

定義 θ = (U,p1,p2,μ0,fF)為信號強度譜中的參數.對于時域頻譜中的第i個頻率分量fi，假設實測數據的譜分析結果為=Im()，信號強度譜模型的預測結果為Ii=I(;θ) ， 則二者的比值Ri=Im(fi)/I(fi;θ)為一個隨機變量.從譜分析理論可知，傅里葉變換后某一特定頻率處的頻譜分量估計符合正態分布[23].若信號強度譜模型是準確的，則Ri的分布為Ri～d，其中參數d= 2M，M為Im(fi)估計中的數據分段數.對上式進行變換，dRi符 合自由度為d的卡方分布，即dRi～χ12,23-24].

對于功率譜的第i個頻率分量fi，給定功率譜模型估計結果Ii情況下，實際信號功率譜測量分量的概率密度函數為

假定功率譜各頻率分量上的測量是統計獨立的，則在給定功率譜模型預測結果時，實測信號功率譜的聯合密度函數為

因此，信號強度譜參數的似然函數為L(θ|Im)=p(Im|I(θ)).式(8)可進一步表示為對數似然函數形式：

3 仿真分析與試驗驗證

Monte Carlo 方法是一種統計試驗方法，通過構造和系統性能近似的概率模型，對系統的隨機性能進行模擬和統計分析[26].由于譜強度、譜指數等參數難以通過實際觀測獲得準確值，無法對譜擬合最優估計結果進行評價，因此利用Monte Carlo 方法對譜參數MLE 方法的性能進行分析.設定合適的參數(包括譜強度、譜指數、漂移速度和中斷頻率)，采用第2 節中的雙參數譜模型，利用相位屏方法可以對閃爍影響下的信號變化進行仿真[20-21,27]；利用仿真獲得的閃爍信號進行譜分析，得到閃爍譜；對閃爍譜進行擬合，并利用MLE 獲得對閃爍譜參數的最優估計；對上述過程進行多次實現(本文為100 次)；對多次MLE 結果與閃爍譜參數預設值(作為準確值)進行比較，分析估計方法的性能.

理論上可對譜模型中的所有參數θ =(U,p1,p2,μ0,fF)進行估計.由于漂移速度可以利用多站組網測量實現[3-4]，因此菲涅爾頻率fF可以很容易獲得.對漂移速度估計的有效性分析將會單獨著文論述，這里僅對譜強度U和譜指數p1,p2估計結果的有效性進行分析，分單指數譜和雙指數譜兩種情況進行.其中，單指數譜可通過設置p1=p2來實現，此時 μ0可取任意值.

3.1 單指數譜估計結果分析

單指數譜分析中，設置譜強度U=0.6，譜指數p1=p2=3， 菲涅爾頻率fF=2 Hz，利用相位屏方法對閃爍信號進行仿真，并對仿真的閃爍信號進行譜參數MLE.

圖1 給出了一次利用相位屏方法仿真獲得的GNSS 閃爍影響下的信號時間變化序列，長度為5 min.由該信號變化序列計算獲得的閃爍指數為S4=0.59.

圖1 仿真的單指數譜GNSS 閃爍信號時間序列Fig.1 Simulated GNSS signal time series under scintillation with one-component power law spectrum

圖2 給出了對圖1 中閃爍信號時間序列的譜分析和譜擬合結果.其中，藍色曲線為仿真數據的譜分析結果.譜分析采用Welch 方法實現，將仿真數據分成5 個數據段，每段長度1 min，對每段數據分別進行譜分析并進行平均獲得最終的譜擬合結果.

圖2 仿真的GNSS 閃爍信號單指數功率譜及MLE 理論功率譜曲線Fig.2 One-component power spectrum of simulated GNSS scintillation signals and that of MLE theory

在譜擬合中進一步考慮噪聲模型可以提高MLE 的穩定性[10-11,22].因此，在理論閃爍譜基礎上增加了簡單的白噪聲模型來進行參數估計.圖2 中的粉色曲線為增加噪聲模型后獲得的功率譜最優估計結果，紅色曲線為不考慮噪聲時的理論譜模型結果.可以看出，通過MLE 獲得的理論譜模型可較好地對實測功率譜進行擬合.

利用MLE 方法對該仿真數據的功率譜進行擬合獲得的最優參數估計結果為U=0.55，p2=3.01，與理論值相近.由式(6)可知，對信號功率譜進行積分可以得到閃爍指數.利用最優譜擬合曲線計算閃爍指數，得到S4=0.55，與實際仿真數據計算結果接近.上述結果說明了MLE 方法在單指數譜形式下估計的有效性.

實測功率譜與理論譜模型的比值dRi符合卡方分布[12,22-25].圖3 給出了利用圖1 數據獲得的功率譜和估計理論譜比值的統計分布.可以看出，dRi統計分布與卡方分布變化一致.統計學中利用Kolmogorov-Smirnov 檢驗(K-S 檢驗)確定觀測數據的潛在概率分布是否與假設分布存在顯著差異，K-S 檢驗門限值為0.05，當觀測數據的概率分布符合已知概率分布時，K-S 檢驗統計量小于預定門限值[28].對dRi分布進一步進行Kolmogorov-Smirnov 測試，得到K-S 統計檢驗量為0.025，說明dRi分布分析結果的正確性.

圖3 單指數譜仿真數據功率譜與MLE 理論譜比值的統計分布Fig.3 Statistic distribution of simulated data for onecomponent power spectrum and that of MLE theory

圖4 給出了每次仿真數據的U和p2最優估計結果及其置信度范圍.可以看出，對于每組數據，最優估計結果在參數真實值附近變化，估計結果的90%置信度范圍均能涵蓋參數真實值，說明了最優估計結果的有效性.

圖4 U 和p2 的MLE 結果及其置信度范圍Fig.4 Estimated MLE results of U and p2 and its confidence intervals

圖5 進一步給出了對譜強度U和譜指數p2所有仿真數據最優參數估計的統計分析結果.可以看出，絕大部分估計參數落入橢圓內，橢圓中心(即最優參數估計的平均值)非常接近參數理論值，表明利用MLE 方法對閃爍單指數功率譜進行最優擬合，可以實現對譜參數的無偏估計.

圖5 U 和p2 的MLE 結果及其協方差橢圓Fig.5 MLE results of U and p2 and its covariance ellipse

3.2 雙指數譜估計結果分析

雙指數譜分析中，設置譜強度U=0.6，譜指數p1=2.5，p2=3.5， 菲涅爾頻率fF=2 Hz，歸一化中斷波數 μ0=5.利用相位屏方法對閃爍信號進行仿真，并對仿真的閃爍信號進行譜參數MLE，獲得的閃爍信號時間變化序列如圖6 所示，長度為5 min，計算相應的閃爍指數S4=0.55.最優估計過程同單指數譜仿真，圖7 給出了圖6中閃爍信號時間序列的譜分析和譜擬合結果，藍色曲線為仿真數據的譜分析結果，粉色曲線為增加噪聲模型后的功率譜最優估計結果，紅色曲線為不考慮噪聲時的理論譜模型結果.可以看出，通過MLE 獲得的理論譜模型可較好地對實測功率譜進行擬合.

圖6 仿真的雙指數譜GNSS 閃爍信號時間序列Fig.6 Simulated GNSS signal time series under scintillation with two-component power spectrum

圖7 仿真的GNSS 閃爍信號雙指數功率譜及MLE 理論功率譜曲線Fig.7 Two-component spectrum of simulated GNSS scintillation signals, and that of MLE theory

利用MLE 方法，對該仿真數據的功率譜進行擬合獲得的最優參數估計結果為U=0.61，p1=2.497，p2=3.480，與理論值相近.利用最優譜擬合曲線計算閃爍指數，得到S4=0.56，與實際仿真數據計算結果接近.上述結果表明MLE 方法在雙指數譜形式下估計的有效性.

圖8 給 出 了 每 次 仿 真 數 據 的U,p1,p2, μ0參數MLE 估計結果及其置信度范圍.可以看出，最優估計結果仍在參數真實值附近變化，但變化范圍要大于單參數譜情況，且絕大部分估計結果的90%置信度范圍能涵蓋參數真實值，表明在雙指數譜情況下仍可以獲得有效的最優估計結果.圖9 給出了對所有仿真數據3 個最優參數估計的統計結果.可以看出，絕大部分參數估計值落入橢圓內，橢圓中心與各參數理論值略有偏差，表明雙參數譜情況下利用MLE 方法仍可以獲得閃爍功率譜參數的近似最優無偏估計.

圖8 U, p1, p2 , μ 0 的MLE 結果及其置信度范圍Fig.8 MLE results and its confidence intervals for U, p1, p2 and μ 0

圖9 U 和p1, p2, μ 0 的MLE 結果及其協方差橢圓Fig.9 MLE of p1, p2 and μ 0 and its covariance ellipses with respect to optimized U estimation

圖10 給出了圖7 數據所得功率譜和最優估計理論譜比值的統計分布.可以看出，dRi統計分布與卡方分布變化一致.對dRi分布進一步進行Kolmogorov-Smirnov 測試，得到K-S 統計檢驗量為0.033，驗證了dRi分析結果的正確性.

圖10 雙指數譜仿真數據功率譜與MLE 理論譜比值的統計分布Fig.10 Statistic distribution of simulated data for twocomponent power law spectrum and that of MLE theory

對每個參數估計時，會存在個別參數估計結果的置信度范圍不能涵蓋參數理論值的現象，如圖8中紅線所示，說明增加參數估計個數會在一定程度上降低參數估計的可信度.但對絕大部分數據參數估計結果是有效的.說明對實際觀測數據進行分析時，可通過選擇不同初始值及增加估計次數，提高對譜參數估計的有效性.

3.3 實測數據結果分析

圖11 為三亞2011-10-24 閃爍影響下的GPS L1 信號時間變化序列，原始數據采樣頻率50 Hz[4].圖中截取了從13:21UT 開始長度為5 min 的閃爍數據，計算獲得的閃爍指數S4=0.67，表明該事件為一次強閃爍事件[29].

圖11 實測的GPS 信號電離層閃爍數據Fig.11 Real GPS scintillation data from measurement

對實測數據利用Welch 方法進行譜分析，數據分成5 個數據段，每段長度1 min.對實測功率譜利用單參數功率模型進行MLE，得到U=1.18,p2=3.59,fF=0.55，各參數90%置信度區間分別為[0.98,1.46]、[3.38, 3.81]和[0.49, 0.63].利用最優譜擬合曲線計算得到S4=0.64，與實際計算值相近.圖12 給出利用最優估計參數獲得的信號強度考慮和沒有考慮噪聲的理論功率譜與實測功率譜，可以看出理論譜可以很好地擬合實際譜.

圖12 實測GPS 信號的閃爍功率譜及其最優擬合結果Fig.12 Spectrum of real GPS scintillation data, and that of theoretical models from optimal estimation

進一步對實測功率譜與最優估計獲得的理論譜進行分析，圖13 給出了實測功率譜和最優估計理論譜比值dRi的 統計分布.可以看出，dRi統計分布符合自由度為10 的卡方分布.對dRi分布進行Kolmogorov-Smirnov 測試，得到K-S 統計檢驗量為0.042，表明對dRi分布分析的正確性.

圖13 實測功率譜與MLE 理論譜比值的統計分布Fig.13 Statistic distribution of real data and that of MLE spectrum

4 結 論

不均勻體參數是閃爍建模預測的重要輸入參量.在雙指數電離層閃爍譜模型基礎上，本文提出針對GNSS 單站常規閃爍觀測，利用對閃爍信號功率譜進行最優擬合，實現對多個譜參數進行MLE 的方法.

利用Monte Carlo 方法分別對單指數譜和雙指數譜情況進行MLE 分析.結果表明：利用MLE 方法，在單指數譜模型情況下可以實現對譜強度和譜指數的無偏估計；在雙指數譜模型情況下，仍可以實現對譜強度和譜指數的近似無偏估計.將最優估計結果帶入閃爍譜模型并計算閃爍指數S4，得到的理論結果與利用實際閃爍信號計算的結果非常接近；利用K-S 檢驗對實際數據功率譜與理論譜比值進行分析，dRi統計分布符合預期的卡方分布.單指數譜模型每次估計結果的90%置信度范圍均能涵蓋參數真實值，雙指數譜模型絕大部分估計結果的90%置信度范圍能涵蓋參數真實值，個別參數估計結果的置信度范圍不能涵蓋參數理論值的現象， 說明隨著估計參數的增多，會出現個別情況下估計結果偏差增大的現象.針對這種情況，實際分析中可以通過對數據進行多次估計并進行平均加以消除.利用實測數據對譜參數最優估計方法進一步分析，強閃爍情況下可以獲得閃爍譜參數的有效估計.同時，仿真數據和實測數據的分析結果表明，通過譜擬合進行閃爍譜參數最優估計的方法適用于中等強度閃爍和強閃爍情況.上述結果表明，通過譜擬合對電離層不均勻體參數進行MLE，是一種有效的閃爍譜參數估計方法.

現有單站電離層閃爍觀測中僅利用高頻采樣測量獲得閃爍指數，本文實現的方法為進一步利用單站觀測數據進行閃爍參數(包括譜強度、譜指數，甚至漂移速度)提取和建模提供了可能.

文中分析的強閃爍事件中沒有發生GNSS 信號失鎖情況，因此可以獲得較好的譜分析結果.若存在信號失鎖或干擾等情況，信號閃爍譜可能會發生畸變，進而影響譜擬合過程及估計的結果.同樣，對于弱閃爍事件也需進一步明確其頻譜是否存在不利于譜擬合最優估計的情況.對基于譜擬合進行閃爍參數最優估計的適用范圍進行分析，將是后繼工作的一個重點.