一種GNSS位置時間序列周期特性挖掘方法與流程

本發明屬于連續運行全球定位導航系統
技術領域：
，涉及一種GNSS位置時間序列周期特性挖掘方法。
背景技術：
：近來的研究表明，GNSS坐標時間序列呈現出明顯的周期性變化(Dongetal.,2002；VanDametal.,2010,2012；Jiangetal.,2013)。傳統的GNSS坐標序列模型認為GNSS坐標序列僅包含周年、半周年項，忽視了其他周期性變化，對一些地球物理現象甚至可能做出錯誤的解釋。如何準確、高效的識別GNSS坐標序列的周期特性，是GNSS坐標時間序列中的關鍵問題之一。目前對GNSS位置時間序列周期特性挖掘一些缺陷：1)什么是周期性的物理起源，GNSS坐標序列中包含哪些周期性信號，缺乏深入的研究；2)目前大部分的周期特性挖掘研究以經典的諧函數為基礎，通過最小二乘方法進行估計，存在一定的局限性；3)大部分研究沒有對坐標序列的周期特性進行深入分析，僅僅確定獲得的周期新信號的周期(主頻)，沒有對其物理屬性、影響因素、起源進行探討。技術實現要素：為了解決上述技術問題，本發明提供了一種GNSS位置時間序列周期特性挖掘方法。本發明所采用的技術方案是：一種GNSS位置時間序列周期特性挖掘方法，其特征在于，包括以下步驟：步驟1：針對原始GNSS觀測值，解算獲取GNSS測站單日松弛解；通過公共基站進行不同解加權進行聯合解算，獲得GNSS測站坐標時間序列及速度參數；步驟2：對獲取的GNSS測站坐標時間序列建立殘差時間序列模型；步驟3：對獲取的GNSS測站坐標時間序列進行預處理；步驟4：對步驟3獲得的GNSS測站坐標時間序列進行地表負載引起的測站位移進行糾正；步驟5：消除GNSS測站非構造運動引起的測站非線性變化及對應的周期性波動，對步驟4獲得的GNSS測站坐標時間序列進行糾正，獲得糾正后的測站坐標時間序列；步驟6：確定小波子序列的主頻率范圍，采用小波分析的方法對糾正后的測站坐標時間序列進行周期特性分析；步驟7：對步驟6中獲得的各子序列進行周期特性分析，分析各子序列的頻譜圖，并根據頻譜圖確定各子序列ENU三坐標序列頻譜圖主峰值對應的頻率，即確定子序列的主頻，再與確定的主頻與步驟6中確定的小波子序列主頻率范圍進行比對，最終確定各子序列的主峰的頻率；步驟8：對確定周期主頻頻段后的子序列進行處理，采用極大似然估計方法進一步確定子序列的噪聲模型、周期性變化振幅，以獲得GNSS坐標序列周期特性具體參數值。與現有的技術相比，本發明具有特點：本發明的創新之處在于，一方面，克服了傳統GNSS利用時間序列諧函數模型進行周期性信號估計的局限性，并對GNSS坐標序列進行粗差、階躍、地表負載及共模誤差糾正，以降低測站的非線性運動及隨機誤差的影響；另一方面本發明對引入小波分析、時頻分析、極大似然估計等方法，提供了一種時頻分析的周期性信號分析方法。該方法依據原始序列中不同信號的固有頻率不相同，從而將原始時間序列分解成不同頻段的新的時間序列，將信號進行分解和重構，對每一層進行處理，從而獲得所需要的部分；此外，本發明專利獎噪聲模型估計方法引入到周期特性分析中，對經過小波分析分解后的信號采用頻率分析結合極大似然估計，更加準確的識別各子序列的周期特性，包括主頻、周年信號振幅、噪聲模型等，有助于真實反映坐標序列的周期性變化，進一步提高GNSS站坐標的精度與可靠性，獲得高精度的位置與速度參數，為深入了解相關地球物理現象的影響機制和變化規律提供依據。附圖說明圖1是本發明實施例的流程圖；圖2是本發明實施例的小波分解原理示意圖。具體實施方式為了使本發明的目的、技術方案和有益效果更加清楚明白，下面結合附圖及具體實施方式，進一步說明本發明。應當理解，一下描述的集體實施方式僅用于解釋本發明，并不用于限定本發明。請見圖1，本發明提供的一種GNSS位置時間序列周期特性挖掘方法，包括以下步驟：步驟1：針對GNSS觀測值及相關文件(星歷文件、表文件等)，本實施例采用高精度GNSS數據后處理軟件以及相應的解算模型解算，分別獲取GNSS測站單日松弛解，通過公共基站進行不同解加權進行聯合解算，獲得GNSS測站坐標時間序列；其中步驟1中GNSS測站坐標時間序列解算方法主要包括2種：1)雙差觀測值進行處理，根據觀測值之間的相關性，可以消除或減弱GNSS觀測手段本身的誤差(如衛星接收機鐘差、電離層延遲等)，從而達到提高精度的目的，以消除隨機誤差引入的周期性新號。基于雙差觀測值的解算軟件主要有GAMIT。2)直接對載波非差觀測量以單點定位(PrecisePointPositioning,PPP)的方式進行處理分析。基于非差觀測值的解算方法具有可用觀測值多，不同測站的觀測值不相關的優點，能更加準確反映測站的真實運動。基于非差觀測值的解算軟件主要有GIPSY、BERNESE、PANDA。3)對雙差、非差的中間解(松弛解)進行聯合解算，可以有效的消除(或減弱)不同軟件及模型引入的隨機及系統誤差，如軟件及其模型的不完善、模型存在的系統偏差等，提高坐標序列的精度。其中步驟1中本發明專利采用GAMIT、GIPSY進行加權進行聯合解算(采用SOPAC給出的GAMIT/GIPSY＝1:2.4經驗權計算)，以獲得最終的GNNS測站單日解位置時間序列(X、Y、Z以及協方差矩陣)。其中步驟1解算過程中所提及解算模型見下表1。表1聯合解解算策略采用的模型固體潮模型IERS2003convention海潮模型FES04model(CMC)極潮模型IERS2003convention投影函數GlobalMappingFunction天線模型絕對IGS相位中心本步驟屬于現有技術，具體可通過現有技術中成熟的GAMIT/GLOBK、Bernese、GIPSY、PANDA、QOCA等高精度數據處理軟件或IGS分析中心獲取數據。不同數據處理軟件由于算法不的完善、模型系統偏差等往往會引入不可避免的解算誤差，本發明的新穎之處在于才用了多種解算軟件(GAMIT、GIPSY等)進行解算，通過公共基站進行不同解加權進行聯合解算，能有效的消除單一軟件解算的模型系統誤差，進一步提高解的可靠性。步驟2：對獲取的測站坐標時間序列建立以下殘差時間序列模型：其中：x(t)E/N/U為時刻t對應的GNSS測站坐標觀測值，包括ENU三坐標分量；x0為測站初始位置，vx為線性速度即趨勢項，t代表時間；代表周期性項階數，為跳變改正項，gj表示跳變振幅，Tgj為跳變發生的時刻即歷元，ng表示跳變個數，j為跳變編號，這里假定發生偏移的時刻Tg已知，H為海維西特階梯函數，在跳變前H值為0，跳變后H值為1；εx(t)為隨機誤差，Ox為粗差)。步驟3：對獲取的GNSS測站坐標時間原始序列繼續進行預處理，包括去均值、粗差剔除(Ox)、階躍改正去均值的目的是便于對數據進行分析，去粗差、去階躍主要為了消除粗差對周期信號的干擾；與傳統方法不同的是，考慮到傳統基于諧函數結合最小二乘擬合去除趨勢項、周年、半周年項的局限性，在對GNSS時間序列進行周期性信號探測時，若事先扣除上述項，會使得探測出的周期性信號存在一定的偏差，因此本發明實施的新穎之處在顧及了傳統的諧函數存在的局限性，依據原始序列中不同信號的固有頻率不相同，從而將原始時間序列分解成不同頻段的新的時間序列，將信號進行分解和重構，對每一層進行處理，從而獲得所需要的部分，有助于真實反映坐標序列的周期性變化。其中步驟3中粗差剔除采用四分位數間距法(interquartilerange,IQR)，四分位數間距由P25、P50、P75將一組變量值等分為四部分，P25稱下四分位數(Q1)，P75稱上四分位數(Q3)，將P75與P25之差定義為四分位數間距(IQR)。分別計算a(Q1-3*IQR)、b(Q3+3*IQR)的值，原始序列中位于(a,b)區間之外的值，則為粗差，IQR方法具有較好的抗差性，能有效剔除GNSS坐標序列中存在的粗差。對于步驟3中GNSS測站坐標時間序列中存在的階躍，采用如下方法進行改正：1)對已知的Offset根據SOPAC發布的發生的時刻及影響進行改正，2)SOPAC未公布的offset采用MedianInterannualDifferenceAdjustedforSkewness(MIDAS)method進行改正。步驟4：對步驟3獲得的GNSS測站坐標時間序列進行地表負載引起的測站位移進行糾正。采用mload程序分別計算大氣壓、非潮汐海洋、積雪和土壤水負荷引起的臺站位移，通過負載改正，提高坐標序列的精度，消除部分非構造信號。在計算地表負載引起的測站位移中采用如下數據模型：大氣負荷采用NCEP的全球表面大氣壓力數據，時間分辨率為6小時，空間分辨率為2.5°×2.5°；非海洋潮汐負荷采用ECCO海洋底壓力模型，時間分辨率為12小時，空間分辨率為空間分辨率為1°×1°；積雪和土壤水負荷數據來自NCEP。步驟5：進行共模誤差剔除，以消除測站非構造運動引起的測站非線性變化及對應的周期性波動。共模誤差分離采用PCA方法進行剔除，其具體實現過程為：假定GNSS臺站獲得的三維坐標觀測值形成一個n×m(n>m，n為觀測數或歷元數，m為觀測類型)的數據矩陣X，其協方差陣為CX，則CX＝XTX。數據矩陣如下：其中：(m×1維列向量)為其協方差陣的特征向量，λi為對應的特征值，令其中σi為正的奇異值，i＝1,2…r。則有：假定其中是n×1列向量，U為n×n向量矩陣，V為m×m向量矩陣，則有：X＝UΣVTCX＝VΛVT即V構成X的正交基底，矩陣X展開可得：ak(ti)可由下式求出：式中ak(t)是第k個主成分(k＝1,2,3)，vk(x)是對應主成分的響應特征矩陣，分別代表時間特征和空間響應，取前k個主分量(k≤4，累積貢獻率大于80％)計算得到的共模誤差為：對原始序列進行糾正，獲得糾正后坐標序列(S)。步驟6：經步驟5后，獲得糾正后的坐標序列(S)，對采用小波分析的方法其進行周期特性分析。首先確定小波子序列的主頻率范圍(附表2)，然后依據原始序列中不同信號的固有頻率不相同，從而將原始時間序列進行小波分解(請見附圖2)成不同頻段的新的時間序列，將信號進行分解和重構，對每一層進行處理，從而獲得各子序列(D1，D2，…Dn，A1，A2，…An)。；表2步驟7：對步驟6獲得的各子序列進行租期特性分析，采用FrequencyAnalysisMappingOnUnusualSampling(FAMOUS)方法分析各子序列的頻譜圖，確定其頻譜主峰值對應的頻率，確定各子序列的主頻，并與步驟6中確定的主頻范圍進行比對，最終確定各子序列的主峰的頻率。步驟8：對確定周期屬性后的子序列處理，采用極大似然估計方法對確定周期主頻頻段后的子序列進行處理，進一步確定其周期特性(噪聲模型、周期性變化振幅等),，以獲得GNSS坐標序列周期特性具體參數值。針對現有技術存在的問題，本發明考慮到已有方法的局限性，進一步對GNSS位置時間序列周期特性進行探討，考慮了傳統的基于諧函數GNSS位置時間序列模型存在的局限性，引入小波分析、時頻分析、極大似然估計等方法，提供了一種GNSS位置時間序列周期特性挖掘方法，高效、準確的分離出GNSS位置時間序列周期特性，有助于真實反映坐標序列的周期性變化，進一步提高GNSS站坐標的精度與可靠性，獲得高精度的位置與速度參數，為深入了解相關地球物理現象的影響機制和變化規律提供依據。應當理解的是，本說明書未詳細闡述的部分均屬于現有技術。應當理解的是，上述針對較佳實施例的描述較為詳細，并不能因此而認為是對本發明專利保護范圍的限制，本領域的普通技術人員在本發明的啟示下，在不脫離本發明權利要求所保護的范圍情況下，還可以做出替換或變形，均落入本發明的保護范圍之內，本發明的請求保護范圍應以所附權利要求為準。當前第1頁1 2 3