本發(fā)明涉及一種確定高精度垂直重力梯度的方法，特別涉及一種利用自由落體式絕對(duì)重力儀觀測數(shù)據(jù)確定垂直重力梯度的方法及系統(tǒng)。

背景技術(shù)：
一直以來，在使用自由落體式絕對(duì)重力儀進(jìn)行絕對(duì)重力測量時(shí)，需要利用高精度相對(duì)重力儀測定垂直重力梯度[1-3]。實(shí)際上，自由落體式絕對(duì)重力儀觀測數(shù)據(jù)，即自由落體運(yùn)動(dòng)的下落時(shí)間-距離序列，包含垂直重力梯度的信息。如果能夠利用絕對(duì)重力測量數(shù)據(jù)確定垂直重力梯度，一方面可以減少實(shí)施絕對(duì)重力測量的工作量和必需的儀器設(shè)備；另一方面由此獨(dú)立測定的垂直重力梯度可用于檢驗(yàn)相對(duì)重力儀的測量結(jié)果。T.M.Niebauer(1995)和H.Hanada(1996)等學(xué)者在介紹FG5自由落體式絕對(duì)重力儀時(shí)就提出了眾多誤差源，分析了觀測數(shù)據(jù)中可能存在的系統(tǒng)誤差[4-6]。K.Charles(1995)對(duì)實(shí)測數(shù)據(jù)進(jìn)行分析，探測到了FG5觀測數(shù)據(jù)中的系統(tǒng)誤差[7]。系統(tǒng)誤差一般由地面的高頻振蕩、測距激光頻率波動(dòng)以及下落棱鏡晃動(dòng)引起。系統(tǒng)誤差中的有色噪聲會(huì)使解算的垂直重力梯度結(jié)果產(chǎn)生偏差(見圖3)。隨后眾多學(xué)者研究了FG5觀測數(shù)據(jù)中系統(tǒng)誤差的探測方法，包括Lomb-Scargle譜分析法[8]和高斯鐘求和法[9]。另外，K.Charles(1995)的研究表明選取不同的下落時(shí)段，解算得到的重力加速度結(jié)果差別很大。同樣，本發(fā)明提出選取不同的下落時(shí)段，解算得到的垂直重力梯度結(jié)果差別也很大(見圖6，如果下落時(shí)段的起始數(shù)據(jù)為第i對(duì)距離-時(shí)間觀測量，終止數(shù)據(jù)為第j對(duì)距離-時(shí)間觀測量，則其對(duì)應(yīng)的編號(hào)為(i-1)*100+j，即圖中顯示了起始數(shù)據(jù)為1-30號(hào)，終止數(shù)據(jù)為601-700號(hào)的3000種時(shí)段選取情況下得到的垂直重力梯度結(jié)果。由圖4在125Hz和155Hz處出現(xiàn)了兩個(gè)峰值)。最早由R.G.Hipkin(1999)成功利用絕對(duì)重力測量數(shù)據(jù)解算出了垂直重力梯度，與相對(duì)重力儀的重力梯度觀測結(jié)果一致，并且精度達(dá)到了3μGal/m(1μGal＝10-8m/s2)[10]。然而，R.G.Hipkin在解算垂直重力梯度時(shí)沒有對(duì)觀測數(shù)據(jù)中可能存在的系統(tǒng)誤差進(jìn)行探測和消除，也沒有給出自由落體運(yùn)動(dòng)下落時(shí)段選取的依據(jù)，所以至今利用絕對(duì)重力測量數(shù)據(jù)解算垂直重力梯度仍然處于實(shí)驗(yàn)階段，而未能形成穩(wěn)定有效的數(shù)據(jù)處理方法。雖然已有眾多學(xué)者研究了對(duì)FG5觀測數(shù)據(jù)中所含誤差進(jìn)行建模的方法，但是其焦點(diǎn)一直集中于確定精度更高、更為可靠的絕對(duì)重力加速度結(jié)果，而沒有進(jìn)一步研究確定垂直重力梯度的方法。R.G.Hipkin(1999)在確定垂直重力梯度的研究中雖然得到了很好的結(jié)果，卻沒有顧及系統(tǒng)誤差和時(shí)段選取等因素的影響，更換另一套FG5觀測數(shù)據(jù)進(jìn)行處理時(shí)，結(jié)果會(huì)有明顯偏差，所以不具有普適性，沒能形成穩(wěn)定有效的數(shù)據(jù)處理方法。所以，至今還沒有形成穩(wěn)定有效的利用自由落體式絕對(duì)重力儀觀測數(shù)據(jù)確定垂直重力梯度的方法。本發(fā)明涉及的參考文獻(xiàn)如下：[1]MicroG-Lacoste:g9User’sManual,2012[OL].http://www.microglacoste.com/pdf/g9Help.pdf[2]XingLelin,LiuHui,LiJianchengetal.ComparisonofAbsoluteGravityMeasurementsObtainedwithFG5/232andFG5/214Instruments[J].Geo-specialInformationScience,2009,12(4):307-310[3]肖凡，王慶賓，張松堂等.FG5-240絕對(duì)重力觀測結(jié)果分析[J].大地測量與地球動(dòng)力學(xué)，2013，33增刊(Ⅰ)：165-167[4]T.M.Niebauer,G.S.Sasagawa,J.E.Falleretal.Anewgenerationofabsolutegravimeters[J].Metrologia,1995,32:159-180[5]H.Hanada,T.TsubokawaandS.Tsuruta.Possiblelargesystematicerrorsourceinabsolutegravimetry[J].Metrologia,1996,33:155-160[6]王勇，張為民.FG5絕對(duì)重力儀[J].地殼形變與地震，1996，16(2)：94-97[7]K.CharlesandR.Hipkin.Verticalgradientanddatumheightcorrectionstoabsolutegravimeterdataandtheeffectofstructuredfringeresiduals[J].Metrologia,1995,32:193-200[8]MartinOrlobandAlexanderBraun.ONTHEDETECTABILITYOFSYNTHETICDISTURBANCESINFG5ABSOLUTEGRAVIMETRYDATAUSINGLOMB-SCARGLEANALYSIS[J].GEOMATICA,2011,66(2):113-124[9]MartinOrlob,AlexanderBraun.ImpactEstimationandFilteringofDisturbancesinFG5AbsoluteGravimeterObservations[J].InternationalJournalofGeosciences,2013,4:302-308[10]R.G.Hipkin.Absolutedeterminationoftheverticalgradientofgravity[J].Metrologia,1999,36:47-52。

技術(shù)實(shí)現(xiàn)要素：
針對(duì)現(xiàn)有方法存在的缺陷，本發(fā)明提供了一種顧及觀測數(shù)據(jù)系統(tǒng)誤差和下落時(shí)段選取標(biāo)準(zhǔn)的，利用自由落體式絕對(duì)重力儀觀測數(shù)據(jù)確定垂直重力梯度的技術(shù)方案。本發(fā)明的技術(shù)方案提供一種基于絕對(duì)重力儀的垂直重力梯度提取方法，包括以下步驟：步驟1，對(duì)多次自由落體的全部觀測數(shù)據(jù)分別解算垂直重力梯度，其中使用單次自由落體的全部觀測數(shù)據(jù)解算垂直重力梯度實(shí)現(xiàn)如下，將單次自由落體整個(gè)下落時(shí)段的下落距離觀測量(xi，ti)全部代入以下觀測方程，其中，xi為第i次檢測所得位置，ti為第i次檢測記錄的時(shí)間，第i次檢測的真實(shí)時(shí)間c為光速；以最小二乘法對(duì)下落距離觀測量進(jìn)行擬合，求解下落質(zhì)體的初始位置x0、初始速度v0、測站重力加速度g0和垂直重力梯度γ；步驟2，根據(jù)步驟1所得對(duì)多次自由落體的全部觀測數(shù)據(jù)分別解算垂直重力梯度結(jié)果，分別計(jì)算單次自由落體的下落距離擬合殘差序列；然后將多次自由落體的下落距離擬合殘差序列進(jìn)行疊加求平均，得到殘差均值序列v；所述計(jì)算單次自由落體的下落距離擬合殘差序列，包括計(jì)算下落距離擬合殘差vi，計(jì)算公式為步驟3，基于步驟2所得殘差均值序列v進(jìn)行建模得到系統(tǒng)誤差模型，獲取系統(tǒng)誤差，實(shí)現(xiàn)如下，設(shè)系統(tǒng)誤差模型形式如下，其中，vs為系統(tǒng)誤差，e為數(shù)學(xué)常數(shù)，k為指數(shù)衰減參數(shù)，t為下落時(shí)間，n為周期項(xiàng)個(gè)數(shù)，r＝1,2,…,n，fr為周期項(xiàng)的頻率，Br和Cr為周期項(xiàng)系數(shù)，等價(jià)于頻率為fr的周期項(xiàng)的振幅和相位；建模過程包括以下子步驟，步驟3.1，將殘差均值序列v按照下落距離均分為m段，每一段的殘差序列記為集合{vj}，j＝1,2,…m，計(jì)算每一段下落距離擬合殘差最大的幅值vsj＝max(vi)，vi∈{vj}，得到各段下落距離擬合殘差最大的幅值vsj及其對(duì)應(yīng)的下落時(shí)刻tsj，代入到衰減項(xiàng)模型中，采用最小二乘法進(jìn)行擬合，估計(jì)衰減參數(shù)k和周期項(xiàng)的幅值A(chǔ)；步驟3.2，從殘差均值序列v中扣除衰減項(xiàng)的影響得到周期項(xiàng)誤差vf，vf＝v/e-kt；將殘差均值序列v記為殘差序列1，將殘差均值序列v中的系統(tǒng)誤差記為殘差序列2，將周期項(xiàng)誤差相應(yīng)的殘差序列記為殘差序列3；對(duì)殘差均值序列v通過插值進(jìn)行重采樣，得到等時(shí)間間隔的殘差均值序列，采用傅里葉變換得到殘差序列3的頻譜圖；步驟3.3，由步驟3.2所得殘差序列3的頻譜圖中峰值的個(gè)數(shù)確定周期項(xiàng)的個(gè)數(shù)，并獲取各周期項(xiàng)的頻率fr，再使用最小二乘法對(duì)殘差序列3進(jìn)行線性擬合，求解出周期項(xiàng)參數(shù)Br和Cr，得到殘差序列2；步驟4，用殘差序列1減去殘差序列2得到系統(tǒng)誤差模型殘差序列，記為殘差序列4，依據(jù)殘差序列4設(shè)定閾值，通過閾值確定參與解算的下落時(shí)段；步驟5，針對(duì)各個(gè)單次自由落體整個(gè)下落時(shí)段改正下落距離觀測量并重新解算垂直重力梯度，包括從下落距離觀測量(xi，ti)中減去步驟3中得到的殘差序列2，得到改正后的下落距離觀測量利用步驟4選取的下落時(shí)段內(nèi)的改正后的下落距離觀測量解算垂直重力梯度；對(duì)多次自由落體的結(jié)算結(jié)果進(jìn)行統(tǒng)計(jì)，得到垂直重力梯度最優(yōu)估值和統(tǒng)計(jì)精度。而且，步驟4中，確定參與解算的下落時(shí)段實(shí)現(xiàn)方式為，以下落時(shí)間總長的至段內(nèi)系統(tǒng)誤差模型殘差序列的絕對(duì)值最大值為閾值，選取系統(tǒng)誤差模型殘差序列中絕對(duì)值不超過該閾值的最長的連續(xù)下落時(shí)段，作為參與解算的下落時(shí)段。而且，步驟5中，對(duì)多次自由落體的結(jié)算結(jié)果進(jìn)行統(tǒng)計(jì)，得到垂直重力梯度最優(yōu)估值和統(tǒng)計(jì)精度，實(shí)現(xiàn)如下，基于多次自由落體解算所得垂直重力梯度，得到均值和標(biāo)準(zhǔn)差，剔除與均值相差超過三倍標(biāo)準(zhǔn)差的結(jié)果，將剩余的可靠結(jié)果的均值作為垂直重力梯度的最優(yōu)估值；將標(biāo)準(zhǔn)差作為單次下落結(jié)果的精度，依據(jù)誤差傳播定律計(jì)算均值精度作為最優(yōu)估值的統(tǒng)計(jì)精度。本發(fā)明還相應(yīng)提供一種基于絕對(duì)重力儀的垂直重力梯度提取系統(tǒng)，包括以下模塊：垂直重力梯度初始解算模塊，用于對(duì)多次自由落體的全部觀測數(shù)據(jù)分別解算垂直重力梯度，其中使用單次自由落體的全部觀測數(shù)據(jù)解算垂直重力梯度實(shí)現(xiàn)如下，將單次自由落體整個(gè)下落時(shí)段的下落距離觀測量(xi，ti)全部代入以下觀測方程，其中，xi為第i次檢測所得位置，ti為第i次檢測記錄的時(shí)間，第i次檢測的真實(shí)時(shí)間c為光速；以最小二乘法對(duì)下落距離觀測量進(jìn)行擬合，求解下落質(zhì)體的初始位置x0、初始速度v0、測站重力加速度g0和垂直重力梯度γ；殘差均值序列提取模塊，用于根據(jù)對(duì)多次自由落體的全部觀測數(shù)據(jù)分別解算垂直重力梯度的結(jié)果，分別計(jì)算單次自由落體的下落距離擬合殘差序列；然后將多次自由落體的下落距離擬合殘差序列進(jìn)行疊加求平均，得到殘差均值序列v；所述計(jì)算單次自由落體的下落距離擬合殘差序列，包括計(jì)算下落距離擬合殘差vi，計(jì)算公式為系統(tǒng)誤差提取模塊，用于基于殘差均值序列v進(jìn)行建模得到系統(tǒng)誤差模型，獲取系統(tǒng)誤差，實(shí)現(xiàn)如下，設(shè)系統(tǒng)誤差模型形式如下，其中，vs為系統(tǒng)誤差，e為數(shù)學(xué)常數(shù)，k為指數(shù)衰減參數(shù)，t為下落時(shí)間，n為周期項(xiàng)個(gè)數(shù)，r＝1,2,…,n，fr為周期項(xiàng)的頻率，Br和Cr為周期項(xiàng)系數(shù)，等價(jià)于頻率為fr的周期項(xiàng)的振幅和相位；并且包括以下子模塊，幅值擬合子模塊，用于將殘差均值序列v按照下落距離均分為m段，每一段的殘差序列記為集合{vj}，j＝1,2,…m，計(jì)算每一段下落距離擬合殘差最大的幅值vsj＝max(|vi|)，vi∈{vj}，得到各段下落距離擬合殘差最大的幅值vsj及其對(duì)應(yīng)的下落時(shí)刻tsj，代入到衰減項(xiàng)模型中，采用最小二乘法進(jìn)行擬合，估計(jì)衰減參數(shù)k和周期項(xiàng)的幅值A(chǔ)；周期項(xiàng)誤差提取子模塊，用于從殘差均值序列v中扣除衰減項(xiàng)的影響得到周期項(xiàng)誤差vf，vf＝v/e-kt；將殘差均值序列v記為殘差序列1，將殘差均值序列v中的系統(tǒng)誤差記為殘差序列2，將周期項(xiàng)誤差相應(yīng)的殘差序列記為殘差序列3；對(duì)殘差均值序列v通過插值進(jìn)行重采樣，得到等時(shí)間間隔的殘差均值序列，采用傅里葉變換得到殘差序列3的頻譜圖；系統(tǒng)誤差模型生成子模塊，用于由殘差序列3的頻譜圖中峰值的個(gè)數(shù)確定周期項(xiàng)的個(gè)數(shù)，并獲取各周期項(xiàng)的頻率fr，再使用最小二乘法對(duì)殘差序列3進(jìn)行線性擬合，求解出周期項(xiàng)參數(shù)Br和Cr，得到殘差序列2；解算時(shí)段確定模塊，用于用殘差序列1減去殘差序列2得到系統(tǒng)誤差模型殘差序列，記為殘差序列4，依據(jù)殘差序列4設(shè)定閾值，通過閾值確定參與解算的下落時(shí)段；垂直重力梯度最優(yōu)估計(jì)模塊，用于針對(duì)各個(gè)單次自由落體整個(gè)下落時(shí)段改正下落距離觀測量并重新解算垂直重力梯度，包括從下落距離觀測量(xi，ti)中減去殘差序列2，得到改正后的下落距離觀測量利用步下落時(shí)段內(nèi)的改正后的下落距離觀測量解算垂直重力梯度；對(duì)多次自由落體的結(jié)算結(jié)果進(jìn)行統(tǒng)計(jì)，得到垂直重力梯度最優(yōu)估值和統(tǒng)計(jì)精度。而且，解算時(shí)段確定模塊中，確定參與解算的下落時(shí)段實(shí)現(xiàn)方式為，以下落時(shí)間總長的至段內(nèi)系統(tǒng)誤差模型殘差序列的絕對(duì)值最大值為閾值，選取系統(tǒng)誤差模型殘差序列中絕對(duì)值不超過該閾值的最長的連續(xù)下落時(shí)段，作為參與解算的下落時(shí)段。而且，垂直重力梯度最優(yōu)估計(jì)模塊中，對(duì)多次自由落體的結(jié)算結(jié)果進(jìn)行統(tǒng)計(jì)，得到垂直重力梯度最優(yōu)估值和統(tǒng)計(jì)精度，實(shí)現(xiàn)如下，基于多次自由落體解算所得垂直重力梯度，得到均值和標(biāo)準(zhǔn)差，剔除與均值相差超過三倍標(biāo)準(zhǔn)差的結(jié)果，將剩余的可靠結(jié)果的均值作為垂直重力梯度的最優(yōu)估值；將標(biāo)準(zhǔn)差作為單次下落結(jié)果的精度，依據(jù)誤差傳播定律計(jì)算均值精度作為最優(yōu)估值的統(tǒng)計(jì)精度。與現(xiàn)有技術(shù)相比本發(fā)明具有以下優(yōu)點(diǎn)和有益效果：1、減小絕對(duì)重力測量的成本。當(dāng)前實(shí)施絕對(duì)重力測量必須使用高精度相對(duì)重力儀測定測站點(diǎn)的垂直重力梯度，才能解算出高精度重力加速度，并將重力加速度正確歸算到參考高度上。本發(fā)明直接利用絕對(duì)重力儀觀測數(shù)據(jù)確定重力垂直梯度，可以避開利用相對(duì)重力儀測量垂直重力梯度的環(huán)節(jié)，節(jié)省了測量所需人力成本和儀器設(shè)備。2、方案的理論依據(jù)充分，穩(wěn)定有效。之前采用絕對(duì)重力測量數(shù)據(jù)確定垂直重力梯度的方案并沒有討論系統(tǒng)誤差的探測和建模以及下落時(shí)段的選取，而這兩項(xiàng)因素對(duì)結(jié)果有較大影響(見圖3和圖6)。本發(fā)明對(duì)這兩項(xiàng)因素進(jìn)行了充分考慮，因而利用該方案確定垂直重力梯度是穩(wěn)定有效的。附圖說明圖1為本發(fā)明實(shí)施例方法的流程圖。圖3為本發(fā)明實(shí)施例周期性有色噪聲對(duì)垂直重力梯度解算結(jié)果的影響示意圖。圖6為本發(fā)明實(shí)施例選取不同下落時(shí)段解算的垂直重力梯度的結(jié)果示意圖。圖2為本發(fā)明實(shí)施例殘差序列1，即多個(gè)殘差序列疊加得到的殘差均值序列示意圖。圖4為本發(fā)明實(shí)施例探測周期性噪聲所得頻譜圖。圖5為本發(fā)明實(shí)施例殘差序列2示意圖。圖7為本發(fā)明實(shí)施例系統(tǒng)誤差模型殘差(殘差序列4)與下落時(shí)段的選取示意圖。圖8為本發(fā)明實(shí)施例解算垂直重力梯度的結(jié)果統(tǒng)計(jì)圖。具體實(shí)施方式以下結(jié)合附圖和實(shí)施例對(duì)本發(fā)明技術(shù)方案進(jìn)行具體描述。本發(fā)明通過將自由落體式觀測數(shù)據(jù)按照常規(guī)方法解算垂直重力梯度時(shí)所產(chǎn)生的下落距離觀測量擬合殘差序列進(jìn)行疊加，獲取了下落距離觀測量中的系統(tǒng)誤差。進(jìn)一步對(duì)系統(tǒng)誤差進(jìn)行分析、建模，利用系統(tǒng)誤差模型改正觀測量，并依據(jù)模型殘差選取可靠的下落時(shí)段，重新解算垂直重力梯度。由此本發(fā)明提出了一種利用自由落體式絕對(duì)重力儀觀測數(shù)據(jù)解算垂直重力梯度的方法。下面以處理FG5絕對(duì)重力儀在某測站上的一套觀測數(shù)據(jù)為例，并結(jié)合圖1說明本發(fā)明實(shí)施例的具體步驟。采用其他以距離和時(shí)間為觀測量的絕對(duì)重力儀時(shí)可以同樣處理。具體實(shí)施時(shí)，可采用計(jì)算機(jī)軟件方式實(shí)現(xiàn)流程的自動(dòng)運(yùn)行。實(shí)施例的流程包括以下步驟：步驟1，對(duì)多次自由落體的全部觀測數(shù)據(jù)分別解算垂直重力梯度，其中使用單次自由落體的全時(shí)段觀測數(shù)據(jù)解算垂直重力梯度實(shí)現(xiàn)如下：將單次自由落體的全部觀測數(shù)據(jù)用于解算，即整個(gè)下落時(shí)段的下落距離-時(shí)間觀測量(xi，ti)全部代入觀測方程中，其中，xi為第i次檢測所得位置，ti為第i次檢測記錄的時(shí)間，而第i次檢測的真實(shí)時(shí)間需要由ti經(jīng)過光路延遲改正得到，即c為光速。以最小二乘法對(duì)下落距離觀測量進(jìn)行擬合，通過迭代計(jì)算可以求解下落質(zhì)體的初始位置x0，初始速度v0，測站重力加速度g0和垂直重力梯度γ等參數(shù)。這一步中涉及的數(shù)據(jù)處理方法在R.G.Hipkin(1999)的研究中已經(jīng)做了詳細(xì)說明，本發(fā)明不予贅述。實(shí)施例中，單次自由落體實(shí)驗(yàn)持續(xù)時(shí)間大約為0.21s，在此期間內(nèi)FG5絕對(duì)重力儀利用激光干涉測距儀和銣原子鐘可以測得落體做自由落體運(yùn)動(dòng)到達(dá)700個(gè)等距離間隔處的時(shí)間。初次解算時(shí)，700對(duì)距離-時(shí)間觀測量(xi,ti)(i＝1,2,…,700)全部代入觀測方程用于解算。步驟2，根據(jù)步驟1所得對(duì)多次自由落體的全部觀測數(shù)據(jù)分別解算垂直重力梯度結(jié)果，分別計(jì)算單次自由落體的下落距離擬合殘差序列；然后將多次自由落體的下落距離擬合殘差序列進(jìn)行疊加求平均，得到殘差均值序列。本發(fā)明進(jìn)一步提出，步驟2中通過疊加距離擬合殘差序列獲取殘差均值序列，并認(rèn)為其中的主要成分為系統(tǒng)誤差：分別計(jì)算單次自由落體的下落距離擬合殘差序列，再將多次下落距離擬合殘差序列進(jìn)行疊加，求得殘差均值序列v(記為殘差序列1)。由于理論上偶然誤差在取均值的過程中得以大幅度削弱，可以認(rèn)為殘差序列1中的主要成分就是距離觀測量的系統(tǒng)誤差(記為殘差序列2)。實(shí)施例中，確定包括垂直重力梯度在內(nèi)的參數(shù)后，計(jì)算下落距離擬合殘差vi，計(jì)算公式為(i＝1,2,…,700)。將多次自由落體實(shí)驗(yàn)得到的距離擬合殘差序列進(jìn)行疊加求平均，得到的殘差均值序列(殘差序列1)中主要包含系統(tǒng)誤差(相應(yīng)殘差序列記為殘差序列2)。本例中FG5提供了2400次自由落體的觀測數(shù)據(jù)，疊加求平均后得到的殘差序列1如圖2所示。步驟3，根據(jù)距離擬合所得殘差均值序列進(jìn)行建模得到系統(tǒng)誤差模型。其中包含以下子步驟：3.1對(duì)殘差均值序列分段求幅值，通過最小二乘法確定衰減項(xiàng)；3.2使用三次樣條插值法對(duì)殘差均值序列重采樣，獲取等時(shí)間間隔的殘差均值序列；3.3使用傅里葉變換探測周期項(xiàng)的個(gè)數(shù)以及各周期項(xiàng)的頻率，采用最小二乘法擬合得到系統(tǒng)誤差模型中的周期項(xiàng)參數(shù)，得到系統(tǒng)誤差模型。具體實(shí)施時(shí)，模型可當(dāng)依據(jù)殘差的特征進(jìn)行選取，由于儀器開始測量階段誤差較大，而一般的外界震動(dòng)、激光不穩(wěn)定帶來的干擾多具有周期性，所以模型中應(yīng)當(dāng)含有衰減項(xiàng)和周期項(xiàng)。本發(fā)明進(jìn)一步提出，衰減項(xiàng)采用指數(shù)衰減項(xiàng)，普適性更佳：根據(jù)殘差序列1進(jìn)行建模，用以獲取殘差序列2。參見圖3，采用余弦噪聲進(jìn)行模擬實(shí)驗(yàn)，振幅分別取為0.1nm、0.5nm和1.0nm，頻率取值從1Hz到150Hz，將噪聲加入到模擬觀測值中，分析其對(duì)解算結(jié)果的影響。本發(fā)明實(shí)施例針對(duì)圖2中殘差序列1的特征，為將殘差均值序列v中呈規(guī)律性變化(衰減趨勢和周期變化)的成分作為系統(tǒng)噪聲提取出來，選用的模型為模型中包含隨時(shí)間變化的指數(shù)衰減項(xiàng)和周期性誤差項(xiàng)。其中，vs為系統(tǒng)誤差，e為數(shù)學(xué)常數(shù)，k為指數(shù)衰減參數(shù)，t為下落時(shí)間，n為周期項(xiàng)個(gè)數(shù)，r＝1,2,…,n，fr為周期項(xiàng)的頻率，Br和Cr為周期項(xiàng)系數(shù)，等價(jià)于頻率為fr的周期項(xiàng)的振幅和相位。實(shí)施例中，為對(duì)殘差系列1進(jìn)行建模，提取殘差序列2，要確定系統(tǒng)誤差的模型參數(shù)，還需要實(shí)施如下子步驟：3.1殘差均值序列v(殘差序列1)按照下落距離均分為m段(具體實(shí)施時(shí)本領(lǐng)域技術(shù)人員可自行設(shè)定m的取值，一般為7-10)，每一段的殘差序列記為集合{vj}(j＝1,2,…m)，計(jì)算每一段下落距離擬合殘差最大的幅值vsj＝max(|vi|)(vi∈{vj})。得到各段下落距離擬合殘差最大的幅值vsj及其對(duì)應(yīng)的下落時(shí)刻tsj，代入到衰減項(xiàng)模型中，采用最小二乘法進(jìn)行擬合，估計(jì)衰減參數(shù)k和周期項(xiàng)的幅值A(chǔ)。實(shí)施例中，m取值為10，由7對(duì)距離-時(shí)間觀測量(vsj,tj)解算得到k＝5.4652。3.2從殘差均值序列v(殘差序列1)中扣除衰減項(xiàng)的影響得到周期項(xiàng)誤差vf(記為殘差序列3)，即vf＝v/e-kt。殘差序列3為多項(xiàng)周期性誤差疊加的結(jié)果。由于自由落體式絕對(duì)重力儀提供的每次下落過程中下落距離觀測量并非等時(shí)間間隔采樣，所以得到的殘差均值序列也不是等時(shí)間間隔采樣，本發(fā)明通過插值進(jìn)行重采樣，得到等時(shí)間間隔的殘差均值序列，插值方法選用三次樣條插值法。然后，采用傅里葉變換得到殘差序列3的頻譜圖。頻譜圖中的一些頻率值對(duì)應(yīng)的振幅較大，即圖中的峰值，表明殘差序列3中具有這一頻率成分。依據(jù)采樣定理，要探測出頻率為f的成分，對(duì)原有數(shù)據(jù)進(jìn)行重采樣的頻率必須高于2f。在圖3中分析了不同振幅和不同頻率的余弦噪聲對(duì)解算結(jié)果的影響，可以看出周期性噪聲對(duì)解算結(jié)果的影響隨振幅的增大而增大，隨頻率的增大而減小，即低頻噪聲的影響較大。實(shí)施例中，重采樣頻率選為1000Hz，用于探測頻率在500Hz以內(nèi)的周期項(xiàng)。3.3由3.2得到的頻譜圖中峰值的個(gè)數(shù)確定周期項(xiàng)的個(gè)數(shù)，并獲取各周期項(xiàng)的頻率fr，再使用最小二乘法對(duì)殘差序列3進(jìn)行線性擬合，求解出周期項(xiàng)參數(shù)Br和Cr。至此，模型參數(shù)全部得以確定，完成了根據(jù)殘差均值序列(殘差序列1)的建模，得到了系統(tǒng)誤差模型，得到了殘差序列2。實(shí)施例中，采用vf經(jīng)過傅里葉變換得到的結(jié)果確定周期項(xiàng)個(gè)數(shù)以及每個(gè)周期項(xiàng)的頻率。參見圖4，可以發(fā)現(xiàn)在殘差序列1中存在頻率為125Hz和155Hz的兩個(gè)周期項(xiàng)。依據(jù)模型對(duì)殘差序列3進(jìn)行擬合，得到周期項(xiàng)參數(shù)Bi和Ci，進(jìn)而確定了vf，即殘差序列3的模型。殘差序列2基于系統(tǒng)誤差模型被確定為vs＝e-5.4652t[0.0524cos(2π×125t)-0.0526sin(2π×125t)+0.00856cos(2π×155t)-0.14sin(2π×155t)]，單位為納米，如圖5所示系統(tǒng)誤差的建模結(jié)果，包含指數(shù)衰減項(xiàng)和兩個(gè)周期項(xiàng)的疊加結(jié)果。殘差序列2即由系統(tǒng)誤差模型確定的離散序列。步驟4，依據(jù)系統(tǒng)誤差模型殘差序列設(shè)定閾值，確定參與解算的下落時(shí)段。本發(fā)明進(jìn)一步提出，步驟4中依據(jù)系統(tǒng)誤差的模型殘差選取參與解算的下落時(shí)段，并提出時(shí)段選取的標(biāo)準(zhǔn)：以下落時(shí)間總長的至段內(nèi)模型殘差的絕對(duì)值最大值為閾值，選取模型殘差的絕對(duì)值不超過該閾值的最長的連續(xù)下落時(shí)段：具體選取下落時(shí)段，一般認(rèn)為在落體開始下落和下落結(jié)束階段的觀測數(shù)據(jù)不可靠，所以應(yīng)當(dāng)選取下落中間時(shí)段的觀測數(shù)據(jù)。用殘差序列1減去殘差序列2得到系統(tǒng)誤差模型殘差序列(記為殘差序列4)。殘差均值序列包含系統(tǒng)誤差vs、偶然誤差和粗差，偶然誤差和粗差構(gòu)成剩余殘差即殘差序列4。首先選取殘差序列4中部一小段殘差(具體選取部分可由本領(lǐng)域人員設(shè)定，建議選取下落時(shí)間總長的到段，認(rèn)為這一段必須參與解算，否則參與解算的數(shù)據(jù)量嚴(yán)重不足)，計(jì)算這一段殘差的絕對(duì)值最大值vmax。以vmax為閾值，在殘差序列4中選取殘差絕對(duì)值不超過vmax的最長的連續(xù)序列，該序列對(duì)應(yīng)的下落時(shí)段即為本發(fā)明提出的應(yīng)該參與解算的下落時(shí)段。實(shí)施例中，計(jì)算殘差序列1與殘差序列2的差值，得到系統(tǒng)誤差模型殘差(殘差序列4)，如圖7所示，圖中兩條橫線為選取時(shí)段采用的閾值，兩條豎線之間即為依據(jù)設(shè)定的閾值應(yīng)當(dāng)選取的下落時(shí)段。選取殘差序列4中下落時(shí)間總長的到段，計(jì)算這一段殘差的絕對(duì)值最大值vmax。以vmax為閾值，在殘差序列4中選取最長的連續(xù)下落時(shí)段，使得該時(shí)段內(nèi)所有殘差的絕對(duì)值都不超過vmax。本例中，先選取殘差序列4中下落時(shí)段0.071s-0.143s對(duì)應(yīng)的殘差段，計(jì)算得到這段殘差的絕對(duì)值最大值0.113nm，依據(jù)該值畫出圖7中兩條橫線，依據(jù)橫線與殘差序列4的交點(diǎn)確定圖中兩條豎線。將圖7中兩條豎線之間的時(shí)段作為參與解算的合適時(shí)段。本例中選取的時(shí)段為0.0534s-0.2111s，相當(dāng)于使用了觀測量(xi,ti)(i＝42,43,…,680)。步驟5，利用系統(tǒng)誤差模型改正下落距離觀測量，利用選取的下落時(shí)段內(nèi)的觀測數(shù)據(jù)解算垂直重力梯度。對(duì)同一測站上的多次下落實(shí)驗(yàn)的結(jié)果進(jìn)行統(tǒng)計(jì)，得到垂直重力梯度最優(yōu)估值和統(tǒng)計(jì)精度。本發(fā)明進(jìn)一步提出，采用系統(tǒng)誤差模型改正下落距離觀測量，然后采用步驟4中選取的下落時(shí)段內(nèi)的觀測數(shù)據(jù)解算垂直重力梯度。對(duì)垂直重力梯度進(jìn)行統(tǒng)計(jì)時(shí)，以三倍標(biāo)準(zhǔn)差為閾值剔除有明顯偏差的結(jié)果。以剩余結(jié)果的平均值作為垂直重力梯度的最優(yōu)估值，以均值的標(biāo)準(zhǔn)差作為結(jié)果的統(tǒng)計(jì)精度：使用步驟3中確定的系統(tǒng)誤差模型，改正下落距離觀測量。選用步驟4中確定的下落時(shí)段，按照步驟1中的方法全部代入觀測方程，以最小二乘法對(duì)下落距離觀測量進(jìn)行擬合，同時(shí)解算重力加速度和垂直重力梯度。每次設(shè)站觀測一般可以獲取多次自由落體的觀測數(shù)據(jù)。由每次自由落體的觀測數(shù)據(jù)可以獨(dú)立地解算出一個(gè)垂直重力梯度。對(duì)多個(gè)垂直重力梯度進(jìn)行統(tǒng)計(jì)，得到均值和標(biāo)準(zhǔn)差。剔除與均值相差超過三倍標(biāo)準(zhǔn)差的結(jié)果，將剩余的可靠結(jié)果的均值作為測站垂直重力梯度的最優(yōu)估值，將其標(biāo)準(zhǔn)差作為單次下落結(jié)果的精度，依據(jù)誤差傳播定律計(jì)算均值精度作為最優(yōu)估值的統(tǒng)計(jì)精度。實(shí)施例中，首先，從下落距離觀測量(xi，ti)中減去步驟3中得到的殘差序列2，即將下落距離觀測量中的系統(tǒng)殘差扣除，得到改正后觀測數(shù)據(jù)，改正后的下落距離觀測量為：然后，使用步驟4確定的下落時(shí)段內(nèi)的改正后觀測數(shù)據(jù)同時(shí)解算重力加速度和垂直重力梯度。由本例中多次自由落體的全部觀測數(shù)據(jù)可以分別解算出2400個(gè)獨(dú)立的垂直重力梯度，對(duì)2400個(gè)結(jié)果進(jìn)行統(tǒng)計(jì)得到垂直重力梯度估值。具體統(tǒng)計(jì)方法為：先求取2400個(gè)垂直重力梯度的均值γ和標(biāo)準(zhǔn)差σ，剔除與均值相差超過3倍標(biāo)準(zhǔn)差的結(jié)果，那么剩余的P個(gè)結(jié)果γp均滿足|γp-γ|＜3σ。對(duì)剩余結(jié)果再次求均值和標(biāo)準(zhǔn)差均值則視為垂直重力梯度的最優(yōu)估值，標(biāo)準(zhǔn)差則視為單次下落數(shù)據(jù)解算結(jié)果的精度。進(jìn)一步由誤差傳播定律可以求得均值精度作為垂直重力梯度均值的統(tǒng)計(jì)精度，統(tǒng)計(jì)結(jié)果如圖8所示。實(shí)施例中，依照R.G.Hipkin(1999)中的方法求得垂直重力梯度的均值為238.86μGal/m，單次下落數(shù)據(jù)解算的結(jié)果精度為376.35μGal/m，均值的統(tǒng)計(jì)精度σaverage為7.88μGal/m。使用相對(duì)重力儀測量垂直重力梯度的方法比較成熟，其結(jié)果也相對(duì)可靠。本例中由相對(duì)重力儀測量的重力垂直梯度為289.0μGal/m，精度優(yōu)于3μGal/m。以相對(duì)重力儀的測量結(jié)果為參考標(biāo)準(zhǔn)，可以認(rèn)為利用本發(fā)明提出的數(shù)據(jù)處理方法解算得到的垂直重力梯度更準(zhǔn)確。綜上所述，本發(fā)明提供了一種穩(wěn)定有效的利用FG5絕對(duì)重力儀觀測數(shù)據(jù)確定垂直重力梯度的方法。主要特征表現(xiàn)在對(duì)觀測量中的系統(tǒng)噪聲進(jìn)行了探測、分析和建模，同時(shí)提出了下落時(shí)段的選取標(biāo)準(zhǔn)。最終通過對(duì)實(shí)測數(shù)據(jù)進(jìn)行處理，驗(yàn)證了該方法的可行性和正確性。具體實(shí)施時(shí)，還可以采用模塊化方式提供相應(yīng)系統(tǒng)。本發(fā)明實(shí)施例提供一種基于絕對(duì)重力儀的垂直重力梯度提取系統(tǒng)，包括以下模塊：垂直重力梯度初始解算模塊，用于對(duì)多次自由落體的全部觀測數(shù)據(jù)分別解算垂直重力梯度，其中使用單次自由落體的全部觀測數(shù)據(jù)解算垂直重力梯度實(shí)現(xiàn)如下，將單次自由落體整個(gè)下落時(shí)段的下落距離觀測量(xi，ti)全部代入以下觀測方程，其中，xi為第i次檢測所得位置，ti為第i次檢測記錄的時(shí)間，第i次檢測的真實(shí)時(shí)間c為光速；以最小二乘法對(duì)下落距離觀測量進(jìn)行擬合，求解下落質(zhì)體的初始位置x0、初始速度v0、測站重力加速度g0和垂直重力梯度γ；殘差均值序列提取模塊，用于根據(jù)對(duì)多次自由落體的全部觀測數(shù)據(jù)分別解算垂直重力梯度的結(jié)果，分別計(jì)算單次自由落體的下落距離擬合殘差序列；然后將多次自由落體的下落距離擬合殘差序列進(jìn)行疊加求平均，得到殘差均值序列v；所述計(jì)算單次自由落體的下落距離擬合殘差序列，包括計(jì)算下落距離擬合殘差vi，計(jì)算公式為系統(tǒng)誤差提取模塊，用于基于殘差均值序列v進(jìn)行建模得到系統(tǒng)誤差模型，獲取系統(tǒng)誤差，實(shí)現(xiàn)如下，設(shè)系統(tǒng)誤差模型形式如下，其中，vs為系統(tǒng)誤差，e為數(shù)學(xué)常數(shù)，k為指數(shù)衰減參數(shù)，t為下落時(shí)間，n為周期項(xiàng)個(gè)數(shù)，r＝1,2,…,n，fr為周期項(xiàng)的頻率，Br和Cr為周期項(xiàng)系數(shù)，等價(jià)于頻率為fr的周期項(xiàng)的振幅和相位；并且包括以下子模塊，幅值擬合子模塊，用于將殘差均值序列v按照下落距離均分為m段，每一段的殘差序列記為集合{vj}，j＝1,2,…m，計(jì)算每一段下落距離擬合殘差最大的幅值vsj＝max(|vi|)，vi∈{vj}，得到各段下落距離擬合殘差最大的幅值vsj及其對(duì)應(yīng)的下落時(shí)刻tsj，代入到衰減項(xiàng)模型中，采用最小二乘法進(jìn)行擬合，估計(jì)衰減參數(shù)k和周期項(xiàng)的幅值A(chǔ)；周期項(xiàng)誤差提取子模塊，用于從殘差均值序列v中扣除衰減項(xiàng)的影響得到周期項(xiàng)誤差vf，vf＝v/e-kt；將殘差均值序列v記為殘差序列1，將殘差均值序列v中的系統(tǒng)誤差記為殘差序列2，將周期項(xiàng)誤差相應(yīng)的殘差序列記為殘差序列3；對(duì)殘差均值序列v通過插值進(jìn)行重采樣，得到等時(shí)間間隔的殘差均值序列，采用傅里葉變換得到殘差序列3的頻譜圖；系統(tǒng)誤差模型生成子模塊，用于由殘差序列3的頻譜圖中峰值的個(gè)數(shù)確定周期項(xiàng)的個(gè)數(shù)，并獲取各周期項(xiàng)的頻率fr，再使用最小二乘法對(duì)殘差序列3進(jìn)行線性擬合，求解出周期項(xiàng)參數(shù)Br和Cr，得到殘差序列2；解算時(shí)段確定模塊，用于用殘差序列1減去殘差序列2得到系統(tǒng)誤差模型殘差序列，記為殘差序列4，依據(jù)殘差序列4設(shè)定閾值，通過閾值確定參與解算的下落時(shí)段；垂直重力梯度最優(yōu)估計(jì)模塊，用于針對(duì)各個(gè)單次自由落體整個(gè)下落時(shí)段改正下落距離觀測量并重新解算垂直重力梯度，包括從下落距離觀測量(xi，ti)中減去殘差序列2，得到改正后的下落距離觀測量利用步下落時(shí)段內(nèi)的改正后的下落距離觀測量解算垂直重力梯度；對(duì)多次自由落體的結(jié)算結(jié)果進(jìn)行統(tǒng)計(jì)，得到垂直重力梯度最優(yōu)估值和統(tǒng)計(jì)精度。本發(fā)明提供了本領(lǐng)域技術(shù)人員能夠?qū)崿F(xiàn)的技術(shù)方案。以上實(shí)施例僅供說明本發(fā)明之用，而非對(duì)本發(fā)明的限制，有關(guān)技術(shù)領(lǐng)域的技術(shù)人員，在不脫離本發(fā)明的精神和范圍的情況下，還可以做出各種變換或變型，因此所有等同的技術(shù)方案，都落入本發(fā)明的保護(hù)范圍。