PCA和KLE在高采樣率GPS定位中的應用*

敖敏思 胡友健 趙 斌 葉險峰 丁開華

(1)中國地質大學信息工程學院，武漢 430074 2)中國地震局地震研究所，武漢 430071 3)武漢大學衛星導航技術研究中心，武漢430079)

PCA和KLE在高采樣率GPS定位中的應用*

敖敏思1)胡友健1)趙 斌2，3)葉險峰1)丁開華1)

(1)中國地質大學信息工程學院，武漢 430074 2)中國地震局地震研究所，武漢 430071 3)武漢大學衛星導航技術研究中心，武漢430079)

為抑制多路徑誤差和隨機噪聲以提高定位精度，引入主成分分析(PCA)和KLE變換方法對定位坐標序列隨機噪聲水平進行評價、提取和消除多天坐標序列中的多路徑誤差。對比分析結果表明，主成分系數比值能有效地反映出強隨機噪聲對某天定位坐標序列的影響，PCA方法能有效地提取和消除多天定位坐標序列中的多路徑誤差，顯著提高定位精度。KLE變換對隨機噪聲污染不敏感，對多路徑誤差的濾波能力較PCA略弱，但其對隨機噪聲以及局部異常具有較好的抑制作用。

高采樣率GPS;多路徑誤差;隨機噪聲;主成分分析;Karhunen-Loeve變換

1 引言

多路徑誤差作為GPS定位中一種非常重要的誤差源，難以采用模型化方法消除其影響。文獻［1］提出利用恒星日濾波方法提取和消除多路徑誤差，但衛星的運行和高采樣率GPS定位的坐標序列中的重復周期并非為一個恒星日，而是隨時間和不同衛星而變化的函數。為更準確地了解衛星以及高采樣率GPS定位的坐標序列中的重復周期而提高定位精度，Choi等［2］提出通過計算多顆觀測衛星的周期取平均值，作為實際周期取代理想的恒星日周期而進行濾波。

文獻［3-5］對此進行了詳細研究，得到一些有用結果。目前現代數字信號處理方法也相繼被引入到觀測值或坐標域內，達到對多路徑誤差進行提取和消除的目的［6-9］。數字信號處理方法的引入有效地提高了提取和消除多路徑誤差的效率，但從本質上來說，數字信號處理方法需要以信號特征頻率的先驗知識為基礎，同時如何盡可能地消除處理過程中參數受到的人為主觀因素的影響，仍然是值得關注的問題。此外，利用信噪比觀測值(SNR)來分離和修正載波相位觀測值中的多路徑誤差［10，11］，也取得了較好的效果，但其要求接收機能同步觀測載波相位和信噪比觀測值。

主成分分析(PCA)和Karhunen-Loeve(KLE)變換作為一種特征信息提取的方法，被廣泛應用于地球物理信號處理、GIS以及遙感領域。本文將討論PCA和KLE在高采樣率GPS定位結果評價和消除多路徑誤差方面的應用。

2 PCA和KLE的原理與方法

PCA和KLE變換的核心思想是通過構造出原始變量的一組線性組合，產生一系列的不相關的新變量，從而提取出能足夠反映原始變量信息的少部分新變量來替代原始變量。對于多天的高采樣率GPS定位坐標序列，在理想狀態下坐標值應等間隔地分布于一條直線上，而實際上，由于受多路徑和隨機噪聲等誤差的影響會產生波動。在利用PCA對n天的高采樣率GPS定位坐標序列進行處理時，設每天的歷元個數為m，將各歷元的解算坐標減去相應的坐標日均值后所得的誤差矩陣為：X(i，j)，(i=1，…，m，j=1，…，n)誤差矩陣中的元素xij表示第j天的第i個歷元的坐標誤差。為不失一般性，假設m＞n。X的協方差矩陣定義為［12-14］

式中，B為n×n的對稱矩陣，可進一步分解為

式中，由列特征向量構成的矩陣V是n×n的正交矩陣，矩陣Λ為由k個特征值構成的非零對角矩陣，實際情況中，一般k=n。由線性代數論可知，誤差矩陣X可由正交函數基V表示為：

式中，A為m×n的矩陣。當特征向量矩陣中的特征向量以特征值降序排列時，矩陣A中的第j列即為第j個分量。矩陣A可由X表示為

由此，原誤差矩陣X的近似矩陣X'可以表示為

式中，Am×l為l個分量構成的矩陣，Vl×n為所對應主成分的系數矩陣。當對特征向量以特征值降序進行排列時，越靠前的分量則代表的是對多天坐標誤差序列協方差貢獻越大的分量。每一個分量所對應的特征值與所有特征值和的比值，即為該分量的貢獻率。也就是說，當誤差矩陣中各天誤差以多路徑誤差為主時，能以靠前的一個或幾個分量的組合來表示多路徑誤差對各天坐標序列的影響。

若利用協方差向量σ對協方差矩陣進行標準化，可得關聯矩陣C，如cij=bij/(σiσj)，其中σ定義如下：

同理，關聯矩陣可以分解為

誤差矩陣X、分量矩陣A和近似誤差矩陣X'也可分別表示為

PCA和KLE變換之間唯一的區別在于前者使用協方差矩陣B，后者采用關聯矩陣C來計算正交矩陣。

若對誤差矩陣X進行如下處理：

將式(11)中的σj定義如式(6)，則Xk的KLE變換結果與原誤差矩陣X的PCA結果相同。

3 研究數據獲取及其預處理

數據的采集在中國地質大學(武漢)64號樓樓頂進行。兩個測站REFF和MP01的觀測時間為2010年年積日361—364天(12月27—30日)連續4天觀測。測站附近無明顯遮擋及反射物。在觀測過程中，測站MP01和REFF均采用TOPCON NETG3A型GPS接收機，均配置TOPCON CR-G3型號雙頻扼流圈天線。采用1 s采樣率進行數據記錄，衛星截至高度角設置為10°(表1)。

表1 數據采集基本設置及參數(單位：m)Tab.1 Settings and parameters of data acquisition(unit：m)

表1中，儀器高度值為控制點至扼流圈天線邊緣的斜線距離。MP1和MP2分別為L1、L2載波上反映偽距碼多路徑效應的指標，采用TEQC軟件計算得出。由測站MP01和REFF的MP1、MP2指標可以發現，MP01測站處的多路徑效應影響遠遠大于REFF處位置。本文采用GAMIT/GLOBK中的TRACK模塊對MP01測站4天的觀測數據進行解算，解算時段為UTC08：55：49—UTC10：09：59，共計4 451個歷元。解算過程中，引入IGS事后精密星歷，以REFF測站為參考站，求取MP01測站4天的動態坐標序列。利用TRACK進行動態定位的關鍵在于模糊度的固定，考慮到基線長度較短，衛星、接收機、電離層、對流層等誤差能較好的通過差分的形式消除，采用L1、L2觀測值分別解算的策略對模糊度進行固定以削弱組合觀測量噪聲對解算結果的影響。在4天的4時段數據中，模糊度修復的比例在95%以上，表明在動態定位的結果中所包含的誤差以多路徑誤差和隨機誤差為主。由于TRACK模塊采用雙差觀測量進行解算，測站REFF的觀測環境造成的誤差會體現在測站MP01的解算結果中。表1中的兩個測站環視條件以及反映多路徑效應的MP1、MP2值可以認為，MP01的誤差主要來自于其西側的墻體。

為提取和消除多路徑誤差，考慮到其在時間上的相關性，采用互相關函數估計高采樣率GPS定位坐標序列的重復周期的偏移。其中，所得年積日第362、363和364較前一天分別提前250 s、249 s和249 s到達。去除提前秒數及各天均值后，從觀測時段內提取的4天的坐標時間序列如圖1所示。

由圖1不難看出，由于受到多路徑效應的影響，N、E方向坐標時間序列出現一定的波動。4天的坐標時間序列形態較為相似(尤其為E方向上)，說明其存在一定的相關性。同時，N方向的波動平緩，而E方向上波動較劇烈，也就是說，相對于NS方向EW方向的多路徑誤差更為明顯。這種情況的出現與之前對測站MP01環視條件進行分析過程中，西側東西走向的墻體是主要的反射源的結論是一致的。

圖1 去除提前秒數及各天均值4天的坐標序列Fig.1 Time series of coordinate on four days with elimination of shift seconds and means

4 隨機噪聲污染的評價

傳統的定位結果評價一般采用標準差或四分位差值(IQR)，其能較好地評價整體定位結果，但不能很好地反映出某一誤差源對定位結果的影響。為此，本文提出利用PCA和KLE變換的方法評價隨機誤差對定位結果的影響。

以高采樣率GPS定位結果為例，當坐標時間序列主要受到多路徑誤差和隨機誤差的影響時，由于多路徑誤差在時間上的相關性，在對4天的坐標時間序列所構成的矩陣進行PCA和KLE變換時，所提取的主成分即為在各天時間序列中的具有較強相關性特征的多路徑誤差。因此，單獨各天對于主成分的響應，即可認為是其受多路徑誤差影響所致。當各天多路徑誤差為坐標時間序列主要誤差源時，各天對于主成分的響應的大小也應差異較小;反之，當某天隨機誤差的影響較大時，將影響當天對于主成分的響應。由式(5)可知，主成分在各天所對應的系數，即可衡量各天對主成分的響應。本例中，4天N、E方向上PCA和KLE變換所得第一主成分及其響應系數如圖2和表2所示。

從圖2不難發現，在沒有強烈異常值影響的情況下，PCA和KLE變換所得到的第一分量具有較強的相似性。

圖2 根據PCA和KLE變換得到的N、E方向第一主分量Fig.2 First principal components in north，east directions from PCA and KLE

表2 利用PCA和KLE得到的4天N、E方向上的主成分響應系數Tab.2 Principal component coefficients on four days in north，east directions from PCA and KLE

由表2可知，當坐標時間序列沒有受到較大的隨機噪聲污染時，PCA和KLE變換所得到的第一主成分的響應系數差異均較小。而KLE變換所得的主成分系數差異較小，其原因在于KLE變換中的標準化過程，對異常的坐標時間序列進行了抑制。為進一步體現隨機誤差對第一主成分相應系數的影響，對第四天的左邊時間序列人為加入不同程度的高斯隨機噪聲。同時，利用第4天對于主成分的響應系數與其他3天的響應系數的均值之比衡量第4天與其他3天之間的響應系數差異。當第四天分別加入不同標準差的高斯白噪聲時，其響應系數比值變化如圖3所示。

由圖3可知，當對第四天N、E方向上不加入或加入較弱高斯噪聲時，響應系數比值約為1，說明其與其他3天響應系數相當;而隨著被高斯誤差污染的程度加劇，響應系數比值逐漸增大。從圖3中來看，當第四天受到原始噪聲的1.5倍標準差的高斯噪聲污染時，該天N、E方向上所得到的主成分的響應系數分別上升到1.5和2.2;隨著高斯白噪聲進一步加入，相應系數也急劇地上升。因此，利用PCA變換后主成分系數的比值，能較好地判斷某天高采樣率GPS定位結果是否受到強隨機誤差的污染。而采用KLE變換時，在所得第一主分量與PCA結果相似的情況下，其響應系數比值則基本不隨所加入高斯誤差的標準差變化而變化，這種現象的出現進一步驗證了KLE變換過程中‘協方差矩陣標準化’這一步驟對異常值起到的抑制作用。

圖3 高斯噪聲的標準差與響應系數比值之間的關系Fig.3 Relationship between response coefficients and standard deviation of added Gaussian noise

5 多路徑誤差的提取和消除

傳統的高采樣率GPS地球物理應用中，一般采用恒星日濾波來消除多路徑誤差。其基本思想是，根據多路徑誤差的恒星日重復性，通過濾波疊加的方法來對其進行提取和消除。其步驟主要包括：1)對高采樣率GPS定位結果進行低通濾波，消除高頻隨機噪聲;2)對濾波后的結果，根據恒星日周期進行等權疊加;3)利用提取的多路徑誤差模型修正高采樣率GPS定位結果。為進一步提高定位精度，提出利用PCA和KLE方法取代等權疊加的算法，從低通濾波后的結果中提取多路徑誤差。在進行PCA和KLE變換時，選擇累積貢獻率達到80%的分量作為代表多路徑誤差的分量，進行提取和消除。一般情況下，第一個主分量或最多兩個主分量的累計貢獻率即達到此標準。仍以前面采用的研究數據為例，采樣窗口長度為7個歷元的均值濾波消除高頻率隨機噪聲。N、E方向均選取前兩個主分量，PCA和KLE變換的累積貢獻率分別為93.9%、98.3%和93.2%、98.3%。采用恒星日濾波和PCA濾波消除多路徑誤差后的高采樣率GPS定位結果N、E方向的誤差分布如圖4所示。其中，任一歷元誤差的L1范數分布定義為殘差絕對值的和除以天數。

由于PCA和KLE變換結果幾乎相同，圖4中僅繪制PCA和恒星日濾波的殘余誤差L1范數分布。由圖4中N、E方向上的殘余誤差L1范數分布可知，PCA濾波和恒星日濾波的殘余誤差差異不大，兩種方法所提取和濾除的多路徑誤差之間存在一致性。但是從整體上來看，PCA濾波的殘余誤差的分布略微低于恒星日濾波，這說明PCA濾波方法利用多天坐標誤差序列之間的相關性來確定其加權值的算法要優于傳統恒星日濾波的等權值疊加算法。在N方向上，誤差總體上較E方向小，絕大部分歷元時刻的誤差L1范數在2.5mm以內。在4天的約第2 500歷元之前時段，誤差的分布均較為平緩，而在2 500歷元之后時段，誤差分布出現了較大的波動。在E方向上，殘余誤差分布的范圍更大一些。在4天的約2 000歷元之前，誤差的分布均較為平緩，而在2 500歷元之后，誤差分布開始出現了波動。從殘余誤差的分布結合圖1中的坐標誤差序列來看，參與誤差分布基本上與4天坐標誤差序列也存在相關性。在N方向上，4天的坐標誤差序列在2 500歷元之前，變換也較為平緩;而在2500歷元之后以及3 300～3 700歷元時段，坐標誤差序列出現跳變或無規律的異常。在E方向上這種情況更加顯著，在第2500歷元左右以及3 200至3 400歷元時段，坐標誤差序列出現無規則的抖動或者跳變。其中，在N、E方向上的某些歷元時刻，PCA濾波和恒星日濾波的殘余誤差均較大，這種情況產生的原因主要是在計算過程中，原始誤差序列在此歷元時刻附近的一個或多個局部異常值起到了主導作用。例如，N方向上的第2 467歷元和2 667歷元(PCA濾波殘余誤差為2.77 mm和2.64 mm，恒星日濾波殘余誤差為2.86 mm和2.85 mm)以及E方向上第2 595歷元時刻(PCA濾波殘余誤差為4.77 mm，恒星日濾波殘余誤差為4.72 mm)，結合圖1看，第一天和第四天均在該段時刻內出現局部性的無規則抖動而形成異常值，造成坐標誤差序列之間重復性和相關性降低從而影響濾波的結果。分別采用恒星日濾波、PCA和KLE變換所得到的N、E方向上共計3 703個歷元的時間序列的平均L1、L2誤差(表3)。

圖4 N、E方向PCA和恒星日濾波殘余誤差L1范數分布Fig.4 North and east L1 norm scatter of coordinate series from PCA and Sidereal filtering

表3中，L2范數分布定義為殘差的平方和除以天數的平方根。由表3不難發現，3種濾波模型均能有效地降低高采樣率GPS定位中的多路徑誤差。從誤差的L2范數分布來看，PCA濾波結果N、E方向上分別下降至40.3%和22.0%，KLE濾波結果N、E方向上平均每天下降至42.6%和22.7%，恒星日濾波則平均每天下降至42.9%和22.7%。表3中的運行時間在2GHz主頻、1.99G內存的Matlab (R2009a)Windows環境下統計得到。總體上來說，3種方法的運行時間差異不大。PCA和KLE基本相同;而恒星日濾波略微低于PCA和KLE方法(約為0.01s)。其原因在于PCA和KLE方法需要計算各天之間的協方差參數確定其權值，而恒星日濾波直接對各天坐標序列進行等權疊加，減少了程序的計算時間。PCA和KLE方法中權值的計算時間取決于坐標誤差序列矩陣的數據量，而矩陣數據量又取決于數據的天數以及時段長度。在一般的地球物理應用中，天數和時段長度一般都不會很長，因此PCA、KLE和恒星日濾波方法在運算時間上不會很大。為進一步討論幾種濾波方法的特征，對4天濾波后的結果進行功率譜分析，由于4天功率譜分析結果相似，僅列出第4天N、E方向的功率譜(圖5)。

表3 恒星日濾波、PCA和KLE處理后N、E方向平均誤差分布及算法運行時間Tab.3 Mean norm scatter and calculation time of PCA and KLE

圖5 年積日364天N、E方向濾波結果的功率譜Fig.5 Power spectrum of north，east components for doy 364

圖5中，根據采樣定理，由于坐標序列的采樣間隔為1s，功率譜曲線的最高頻率截至0.5 Hz。由圖5中的原始坐標誤差序列的功率譜曲線可知，多路徑誤差主要集中在0.000 3～0.01 Hz的頻帶內。當頻率高于0.01 Hz后，N、E方向上的功率譜幅度分別下降約2個和4個數量級。在N和E方向上恒星日濾波、PCA和KLE變換均能有效地減少多路徑誤差的影響，提高高采樣率GPS的定位精度。3種濾波方法主要作用于0.1 Hz以下的頻段。而3種濾波方法的差異則也主要體現在0.000 3～0.01 Hz的頻段。其中，KLE濾波和恒星日濾波的結果較為相似，尤其在E方向上，兩種方法的功率譜曲線基本重合，而PCA濾波方法則略微優于KLE和恒星日濾波。PCA濾波方法的這種優勢與表3中的L1、L2誤差范數是一致的。

5 結語

在顧及多路徑誤差在時間上具有重復性特征的基礎上，提出了將PCA和KLE變換的方法引入高采樣率GPS定位中，對隨機噪聲水平進行評價以及提取和消除多天坐標序列中的多路徑誤差。通過對實際觀測數據的處理以及與恒星日濾波方法進行對比分析的結果表明，利用PCA方法中的主成分系數比值能有效地反映出強隨機噪聲對某天定位坐標序列的影響，同時，該方法能有效地提取和消除多天定位坐標序列中的多路徑誤差，從而提高定位精度。在KLE變換與PCA方法相比較，其主成分系數比值對隨機噪聲污染不敏感，對多路徑誤差的消除能力與PCA相比略弱，但其對隨機噪聲以及局部異常能起到較好的抑制作用。在實際應用中，若未知多路徑誤差是否為主要誤差源或已知多路徑誤差為主要誤差源的情況下，可以利用PCA方法對隨機噪聲水平進行評價或直接進行濾波;若已知某天數據隨機噪聲水平較高的情況下，則可以利用KLE方法進行濾波以提高定位精度。

1 Bock Y，et al.Instantaneous geodetic positioning at medium distances with the global positioning system［J］.J Geophy Res.，2000，105(B12)：28 223-28 253.

2 Choi K，et al.Modified sidereal filtering：implications for high-rate GPS positioning［J］.Geophys Res Lett.，2004，31：L22608.

3 Larson K M，Bilich A and Axelrad P.Improving the precision of high-rate GPS［J］.J Geophys Res.，2007，112：B05422.

4 Ragheb A E，Clarke P J，Edwards S J.GPS sidereal filtering：coordinate-and carrier-phase-level strtegies［J］.J Geod.，2007，81(5)：325-335.

5 Ping Z，et al.Sidereal filtering based on single differences for mitigating GPS multipath effects on short baselines［J］.J Geod.，2010，84(2)：145-158.

6 Satirapod C and Rizos C.Multipath mitigation by wavelet analysis for GPS base station applications［J］.Serv Rev.，2005，38(295)：2-10.

7 Chang C and Juang J.An adaptive multipath mitigation filter for GNSS applications［R］.Journal on Advances in Signal Processing，2008，214815.

8 Zheng D，et al.Filtering GPS timeseries using a Vondrak filter and cross-validation［J］.J Geod.，2005，79(6/7)：363 -369.

9 戴吾蛟，等.基于經驗模式分解的濾波去噪法及其在GPS多路徑效應中的應用［J］.測繪學報，2006，35(4)：321-327.(Dai Wujiao，et al.EMD filter method and its application in GPS multipath［J］.Acta Geodaetica et Catographica Sinica，2006，35(4)：321-327)

10 Bilich A，Larson K M and Axelrad P.Modeling GPS phase multipath with SNR：Case study from the Salar de Uyuni，Boliva［J］.J Geophys Res.，2008(113)：B04401.

11 吳雨航，等.利用信噪比削弱多路徑誤差的方法研究［J］.武漢大學學報(信息科學版)，2008，33(8)：842-845.(Wu Yuhang，et al.Mitigation of multipath effect using SNR values［J］.Geomatics and Information Science of Wuhan University，2008，33(8)：842-845)

12 Dong D，et al.Spatiotemporal filtering using principal component analysis and Karhunen-Loeve expansion approaches for regional GPS network analysis［J］.J Geophys Res.，2006，(111)：B03405.

13 伍吉倉，孫亞鋒，劉朝功.連續GPS站坐標序列共性誤差的提取與形變分析［J］.大地測量與地球動力學，2008，(4)：97-101.(Wu Jicang，Sun Yafeng and Liu Chaogong.Extraction of common mode error for continous GPS network and deformation analysis［J］.Journal of Geodesy and Geodynamics，2008，(4)：97-101)

14 胡守超，伍吉倉，孫亞鋒.區域GPS網三種時空濾波方法的比較［J］.大地測量與地球動力學，2009，(3)：95-99.(Hu Shouchao，Wu Jicang and Sun Yafeng.Comparison among three spatiotemporal filtering methods for regional GPS network analysis［J］.Journal of Geodesy and Geodynamics，2009，(30：95-99)

APPLICATION OF PCA AND KLE TO HIGH-RATE POSITIONING

Ao Minsi1)，Hu Youjian1)，Zhao Bin2，3)，Ye Xianfeng1)and Ding Kaihua1)

(1)Faculty of Information and Engineering，China University of Geosciences，Wuhan 430074 2)Institute of Seismology，China Earthquake Administration，Wuhan 430071 3)Research Center of GNSS，Wuhan University，Wuhan430079)

Multipath error and random noise are two important error sources in high-rate GPS positioning.In order to characterize and mitigate these errors，and improve the accuracy of high-rate GPS positioning，the Principal Component Analysis(PCA)and Karhunen-loeve Expansion(KLE)are introduced to evaluate the random noise level of the time series of coordinate on some days，extract and mitigate the multipath error from coordinates time series of multiple days.The data processing，comparison and analysis based on real data show that the principal component coefficients ratio of PCA can effectively reflect the random noise level of the time series of coordinates on certain day，meanwhile，the multipath error of the time series of coordinates on many days can be extracted and mitigated with introduction of PCA.Compared to PCA，the KLE approach performs slightly weaker on accurate improvement，but the random noise and local anomaly can be suppressed by KLE better since it is not sensitive to the random noise.

high-rate GPS;multipath error;random noise;principal component analysis;Karhunen-Loeve expansion

1671-5942(2011)06-0145-06

2011-06-15

國家自然科學基金(40974002);國家重大科技基礎設施建設項目(發改高技［2007］1911)

敖敏思，男，1985年生，博士研究生，研究方向為GNSS技術及其地學應用、變形監測.E-mail：aominsi@163.com

P228.42

A