本發(fā)明屬于時間頻率傳遞的技術領域，尤其涉及一種利用小波奇異值分解和自適應指數(shù)平滑進行gps載波相位周跳探測與修復的方法。

背景技術：

時間頻率的研究是基礎科學研究中的一個重要分支，時間頻率在科研、定位系統(tǒng)、電力系統(tǒng)、軍事和國家安全等方面處于舉足輕重的地位。世界各國大都建有自己的守時實驗室，產生本國的原子時標，并參與國際比對。

中國計量科學研究院承擔著向地方實驗室發(fā)布標準時間與頻率，并參與國際原子時比對的任務。因此，與國際原子時比對，以獲得穩(wěn)定度及準確度更高的本地原子時是至關重要的。中國計量科學研究院目前采用的時間比對方法為gps共視法，但共視法要求兩地的實驗室在同一時刻需觀測到同一顆衛(wèi)星，且持續(xù)較長時間。若兩地實驗室距離較遠，無法觀測到同一顆衛(wèi)星，則需第三個實驗室作為中繼。gps載波相位法無需考慮觀測衛(wèi)星的問題，且時頻傳遞精度更高。但gps載波相位時頻傳遞中有一個重要的問題，即周跳。周跳是指由于gps接收機失鎖，或接收載波由于其他原因被阻斷，而產生間斷的現(xiàn)象。周跳的發(fā)生對gps載波相位時頻傳遞影響較大。故而周跳探測與修復是非常重要的。本課題研究gps載波相位周跳探測與修復方法，旨在保證時頻傳遞的精度，并提高本地原子時的穩(wěn)定度和準確度。

技術實現(xiàn)要素：

本發(fā)明要解決的技術問題是，提供一種利用小波奇異值分解和自適應指數(shù)平滑進行gps載波相位周跳探測與修復的方法，利用gps觀測文件，采用小波奇異值分解方法，探測出gps載波相位傳遞中，整周跳變的現(xiàn)象以及跳變發(fā)生的歷元，并采用自適應指數(shù)平滑法以修復周跳，保證gps載波相位時頻傳遞的精度。

為解決上述問題，本發(fā)明采用如下的技術方案：

一種gps載波相位周跳探測與修復的方法，包括以下步驟：

第一步、利用matlab讀取rinex型的gps觀測文件，得到一天內觀測到的衛(wèi)星號、gps偽距觀測值和相位觀測值；

第二步、選擇偽距觀測值以及相位觀測值，由于同一個歷元，gps接收機所接收到的觀測值來自多顆衛(wèi)星，且跟蹤同一顆衛(wèi)星的時間長短也各有不同，為提高時頻傳遞的精度，需盡可能的得到持續(xù)歷元較長的觀測值，首先選擇觀測到的某一顆衛(wèi)星，然后選擇持觀測歷元的觀測值，保留相位觀測值和偽距觀測值；

第三步、采用mw組合法(melboutne-wubbena)構造初始檢測量，gps信號是調制在頻帶為l1，l2上的，設第i歷元的偽距觀測值為p1，p2，相位觀測值為寬巷模糊度為那么：

其中，λ1＝19.06cm，λ2＝24.45cm，為l1，l2的波長，f1＝1575.42mhz，f2＝1227.60mhz，為l1，l2的頻率，為寬巷波長，約為86.2cm，c為光速。

在歷元間對mw組合量做差，即可得到：

第四步、構造出的mw組合差分量是一個一維序列，利用matlab小波分解函數(shù)(wavedec函數(shù))對此mw組合差分量進行小波分解，得到各層的細節(jié)分量，以四層小波分解為例，設得到細節(jié)分量為d1，d2，d3，d4，細節(jié)分量是與mw組合差分量長度相同的一維序列；

第五步、利用各層細節(jié)分量構造hankel矩陣，對所述矩陣進行奇異值分解，得到它們的奇異值矩陣，對奇異值進行差分譜和能量比分析，選擇合適的奇異值進行降秩重構，濾除分量中的噪聲，即可判斷出發(fā)生周跳的歷元；

第六步、利用未發(fā)生周跳的多個歷元，進行自適應指數(shù)平滑方法預測出未來部分歷元的值，將已修復的歷元加入歷史歷元中，參與預測，修復下一個發(fā)生整周跳變的歷元。

作為優(yōu)選，第六步中，自適應指數(shù)平滑法通過對預測目標歷史統(tǒng)計序列進行逐層的平滑計算，找出預測目標的基本變化趨勢并以此預測，

mw組合差分量用式：計算，設選擇持續(xù)觀測歷元數(shù)為t,組合差分量值用v1，v2，v3，...，vt表示，則一次指數(shù)平滑公式為：

式中：為第t周期的一次指數(shù)平滑值，α為加權系數(shù)，0＜α＜1，

預測模型為

為預測值，

二次指數(shù)平滑即為對一次指數(shù)平滑的結果做指數(shù)平滑。

采用三次指數(shù)平滑，即對二次指數(shù)平滑的結果做指數(shù)平滑，計算公式為：

三次指數(shù)平滑法的預測模型為：

其中：

使用誤差平方和：表征預測精度，

采用0.618優(yōu)選法確定最終的預測平滑系數(shù)α，0＜α＜1，首先選擇α的值為0.618，使用三次指數(shù)平滑模型進行預測，利用預測出的結果計算誤差平方和，第二次選擇(1-0.618+0＝0.382)為α的值，同樣的方法進行預測，計算誤差平方和，若第二次預測的誤差平方和小于第一次預測，則去除0.618以上的部分，反之去除0.618一下的部分，對得到的新的區(qū)間繼續(xù)進行0.618優(yōu)選法選擇α的值，直到選擇出最優(yōu)的α。

本發(fā)明的特征如下：

(1)能夠準確的判斷出多個歷元內是否有周跳發(fā)生，和周跳發(fā)生的歷元。

(2)能夠實時的進行周跳探測。

(3)可探測出一周的小周跳。

與現(xiàn)有技術相比，本發(fā)明具有以下有益效果：

本發(fā)明提出了一種基于小波奇異值分解的gps載波相位周跳探測的方法，該方法與現(xiàn)行的電離層殘差法、多項式擬合法和mw組合法等相比，可以準確探測出較小周跳的發(fā)生，漏檢率大大降低，為gps載波相位時頻傳遞精度提供了保障。采用自適應指數(shù)平滑法預測發(fā)生周跳歷元的值以修復周跳，結果表明，預測值與實際值誤差較小，精度較高。

附圖說明

圖1gps載波相位周跳探測流程圖；

圖2mw組合法探測周跳結果圖；

圖3奇異值能量比示意圖；

圖4奇異值差分譜圖；

圖5基于小波分解和奇異值分解探測周跳結果圖；

圖6自適應指數(shù)平滑法預測結果圖。

具體實施方式

以下結合具體實施例，并參照附圖，對本發(fā)明進一步詳細說明。

如圖1所示，本發(fā)明實施例提供一種利用小波奇異值分解和自適應指數(shù)平滑進行gps載波相位周跳探測與修復的方法，包括以下步驟：

步驟1、對rinex文件進行處理。

使用matlab讀取rinex文件，得到一天內觀測到的所有衛(wèi)星和其觀測值等信息。觀測值一般包括相位觀測值、c/a碼偽距觀測值、p碼偽距觀測值、多普勒頻率等。

步驟2、從步驟1所得到的衛(wèi)星中選擇某一顆衛(wèi)星，然后選擇持續(xù)觀測歷元的觀測值，保留相位觀測值和p碼偽距觀測值。

步驟3、構造mw組合差分量。

gps信號是調制在頻帶為l1，l2上的，設第i歷元的偽距觀測值為p1，p2，相位觀測值為第i個歷元的寬巷模糊度可用下式表示：

其中，λ1＝19.06cm，λ2＝24.45cm，為l1，l2的波長，f1＝1575.42mhz，f2＝1227.60mhz，為l1，l2的頻率，為寬巷波長，約為86.2cm，c為光速。

在歷元間對其做差，即可得到：

以此作為后續(xù)周跳探測的檢測量。此檢測量為一維序列，本發(fā)明以長度為400的序列為例。

步驟4、利用matlab函數(shù)(wavedec)對mw組合差分量進行四層小波分解。得到四層細節(jié)分量，分別為d1，d2，d3，d4，四層細節(jié)分量長度都為399。設d1＝(x1，x2，...，x399)

步驟5、利用細節(jié)分量構造hankel矩陣。

對于一個一維的信號序列，為對其進行奇異值分解處理，必須首先構造一個矩陣。本文采用hankel矩陣的生成方式。hankel矩陣的構造形式如下：

hankel矩陣的特點為，其反對角線上的元素相同。若利用時間序列構造hankel矩陣，則下一行矢量元素僅比上一行元素滯后一個時間點。

對于mw組合差分量進行小波分解重構得到細節(jié)分量，其長度為n，利用此序列構造hankel矩陣，若n為偶數(shù)，則令此矩陣的行數(shù)m＝n/2+1，列數(shù)n＝n/2+1；若n為奇數(shù)，則令矩陣行數(shù)m＝(n+1)/2，列數(shù)n＝(n+1)/2。

本發(fā)明以長度400的序列為例，細節(jié)分量長度為399。構造的hankel矩陣行數(shù)為200，列數(shù)為200。d1，d2，d3，d4構造的矩陣為m1，m2，m3，m4。以第一層分量d1＝(x1，x2，...，x399)為例，構造出的m1為：

步驟6、對構造成的hankel矩陣進行奇異值分解。

奇異值分解(singularvaluedecomposition,svd)是一種矩陣的正交化分解的方法。設m是一個m×n階的實矩陣，則必定存在一個正交矩陣u∈r^m×m和另一個正交矩陣v∈r^n×n，使得

m＝udv^t(4)

式中，d是一個半正定對角矩陣，且d∈r^m×n，稱為奇異值矩陣，矩陣d可表示為：

式中，矩陣s＝diag(σ1，σ2，σ3，…，σp)，其中p＝min(m，n)，且σ1≥σ2≥σ3≥…≥σp，對角線上的元素為矩陣m的奇異值，奇異值的數(shù)目及為此矩陣的秩。

本發(fā)明中進行奇異值分解的矩陣為200×200的矩陣，得到其奇異值矩陣d大小為200×200，即得到200個奇異值。

步驟7、奇異值能量比和差分譜分析。

將奇異值由大到小排列成一個序列，每個奇異值的能量比和差分定義為：

δσi＝σi+1-σi(6)

步驟8、降秩重構

結合差分譜和能量比，保留快速下降及之前的奇異值進行降秩重構，即可去除噪聲。

以第一層分量為例，如圖3、4所示，本例保留前60個奇異值，原奇異值矩陣d對角線上有200個不為零的數(shù)據(jù)，即矩陣m有200個奇異值。現(xiàn)保留前60個較大的奇異值重新構造奇異值矩陣d，奇異值依舊排列在對角線上，對角線的后140個值為零。再利用式:m＝udv^t重構出新的矩陣m。再根據(jù)hankel矩陣m的特點，還原出長度為399的細節(jié)分量。

重構后的細節(jié)分量如圖5所示，100、200、300歷元有峰值，即可判斷出周跳發(fā)生在此三個歷元上。

步驟9、指數(shù)平滑法修復周跳。

本發(fā)明中，以第一個周跳為例，100歷元發(fā)生周跳，則利用前99個數(shù)據(jù)進行預測，預測出的100歷元值取代實際值，即可修復周跳。本發(fā)明采用自適應指數(shù)平滑法進行周跳修復。

指數(shù)平滑法通過對預測目標歷史統(tǒng)計序列進行逐層的平滑計算，消除由于隨機因素造成的影響，找出預測目標的基本變化趨勢并以此預測。

mw組合差分量用式：計算，設選擇持續(xù)觀測歷元數(shù)為t,組合差分量值用v1，v2，v3，...，vt表示，則一次指數(shù)平滑公式為：

式中：為第t周期的一次指數(shù)平滑值，α為加權系數(shù)，0＜α＜1。

預測模型為

為預測值。

二次指數(shù)平滑即為對一次指數(shù)平滑的結果做指數(shù)平滑。

本發(fā)明采用三次指數(shù)平滑，即對二次指數(shù)平滑的結果做指數(shù)平滑。計算公式為：

三次指數(shù)平滑法的預測模型為：

其中：

本發(fā)明使用誤差平方和：表征預測精度。

步驟10、確定最終的平滑系數(shù)。

本發(fā)明采用0.618優(yōu)選法確定最終的預測平滑系數(shù)α。0＜α＜1,首先選擇α的值為0.618，使用三次指數(shù)平滑模型進行預測，利用預測出的結果計算誤差平方和。第二次選擇(1-0.618+0＝0.382)為α的值，同樣的方法進行預測，計算誤差平方和，若第二次預測的誤差平方和小于第一次預測，則去除0.618以上的部分，反之去除0.618一下的部分，對得到的新的區(qū)間繼續(xù)進行0.618優(yōu)選法選擇α的值，直到選擇出最優(yōu)的α。

利用最優(yōu)的α的值進行三次指數(shù)平滑預測，預測到的第100個值取代原始值，即可修復周跳，如圖6所示。

本發(fā)明的特點在于：

1、利用小波分解，得到mw組合差分量的各層細節(jié)分量更好的反映出周跳發(fā)生的現(xiàn)象。

2、利用svd分解對細節(jié)分量進行降秩重構，一定程度上減少了噪聲對于周跳探測的影響，大大降低了漏檢率。

3、利用自適應指數(shù)平滑法預測，以修復周跳。

以上實施例僅為本發(fā)明的示例性實施例，不用于限制本發(fā)明，本發(fā)明的保護范圍由權利要求書限定。本領域技術人員可以在本發(fā)明的實質和保護范圍內，對本發(fā)明做出各種修改或等同替換，這種修改或等同替換也應視為落在本發(fā)明的保護范圍內。