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