王文超 王夏莉 王 斌 耿曉燕 劉煥玲

1 國家自然資源部大地測(cè)量數(shù)據(jù)處理中心，西安市友誼東路334號(hào)，7100542 中國測(cè)繪科學(xué)研究院，北京市蓮花池西路28號(hào)，100036

帶有齒輪系統(tǒng)及量測(cè)螺旋系統(tǒng)的LCR-G型和Burris型相對(duì)重力儀由于其特殊的結(jié)構(gòu)，會(huì)存在周期誤差，該誤差主要由齒輪偏心、齒輪齒距分度誤差、量測(cè)螺桿螺紋距誤差等機(jī)制產(chǎn)生[1]。短周期項(xiàng)對(duì)重力基準(zhǔn)網(wǎng)解算精度的影響平均約在3 μGal以下，而中長(zhǎng)周期項(xiàng)的影響最大超過20 μGal[2-3]。因此，重力基準(zhǔn)網(wǎng)的解算必須考慮周期項(xiàng)的影響。目前，重力基準(zhǔn)網(wǎng)對(duì)周期項(xiàng)的解算主要是利用其周期項(xiàng)的展開式。一般情況下需考慮7個(gè)周期項(xiàng)，但是當(dāng)7個(gè)周期項(xiàng)全部考慮時(shí)，由于參數(shù)過多會(huì)發(fā)生過擬合現(xiàn)象，且7個(gè)周期項(xiàng)全部解算精度要達(dá)到2 μGal以上大約需要200多個(gè)觀測(cè)值[3-4]，但在實(shí)際觀測(cè)中很難保證每臺(tái)儀器均獲得200個(gè)以上的觀測(cè)量。因此，相對(duì)重力儀周期項(xiàng)的選取是重力基準(zhǔn)網(wǎng)解算的一個(gè)難題。

目前相對(duì)重力儀周期項(xiàng)的選取主要采用以下2種方法：1)根據(jù)誤差傳播定律，得出周期項(xiàng)的取舍標(biāo)準(zhǔn)[5]，然后對(duì)周期項(xiàng)進(jìn)行篩選；2)采用逐步回歸分析法[6]，利用多元回歸方程中各變量方差的貢獻(xiàn)進(jìn)行顯著性檢驗(yàn)，根據(jù)檢驗(yàn)結(jié)果對(duì)周期項(xiàng)進(jìn)行保留或剔除。但這2種方法在篩選過程中均需要大量的觀測(cè)量來支撐。在機(jī)器學(xué)習(xí)領(lǐng)域，有學(xué)者使用正則化方法[7-8]對(duì)多項(xiàng)式回歸方程的冗余項(xiàng)進(jìn)行剔除，結(jié)果表現(xiàn)良好，能有效防止過擬合現(xiàn)象。因此，本文基于2000國家重力基本網(wǎng)及2020最新國家重力基準(zhǔn)網(wǎng)的觀測(cè)數(shù)據(jù)，使用正則化方法，以LCR-G型相對(duì)重力儀為例，研究國家重力基準(zhǔn)網(wǎng)中相對(duì)重力儀周期項(xiàng)參數(shù)的確定方法。

1 周期誤差標(biāo)定的原理

1.1 周期誤差數(shù)學(xué)模型

相對(duì)重力儀周期項(xiàng)的函數(shù)模型[9]為：

(1)

式中，An為各周期項(xiàng)的振幅，φn為各周期項(xiàng)的初相，R為儀器讀數(shù)，Tn為相對(duì)重力儀的傳動(dòng)周期，一般情況下，n≤7。

由于該函數(shù)模型的形式不適合平差過程中的誤差方程組建，因此，對(duì)其進(jìn)行推導(dǎo)變換。令Xn=Ancosφn、Yn=-Ansinφn，將Xn和Yn作為誤差參數(shù)代入式(1)得：

(2)

在重力基準(zhǔn)網(wǎng)平差解算過程中，相對(duì)重力儀的參數(shù)結(jié)果分為以下幾種情況：

1)如果Xn=0、Yn=0，則振幅An為0，該周期項(xiàng)可以忽略；

1.2 平差模型

重力基準(zhǔn)網(wǎng)平差采用間接平差模型，將經(jīng)過固體潮改正、海潮改正、氣壓改正、儀器高改正、零漂改正后測(cè)段起、止點(diǎn)的最后觀測(cè)值和儀器讀數(shù)作為觀測(cè)量，其誤差方程為[4]：

(3)

式中，gi、gj分別為測(cè)站i、j點(diǎn)平差后的重力值，gRZi、gRZj分別為測(cè)站i、j點(diǎn)經(jīng)過5項(xiàng)改正的相對(duì)聯(lián)測(cè)最后觀測(cè)值，Ri、Rj分別為儀器在測(cè)站i、j點(diǎn)的觀測(cè)讀數(shù)，CK為重力儀的M次多項(xiàng)式格值函數(shù)的K次格值改正因子，Xn、Yn為重力儀周期誤差參數(shù)，Tn為周期誤差的周期，其周期值大小見表1，1個(gè)計(jì)數(shù)單位(1格)相當(dāng)于mGal。

表1 LCR-G型重力儀傳動(dòng)周期

2 周期項(xiàng)確定方法

2.1 正則化函數(shù)模型

正則化函數(shù)模型[7-8]公式為：

Rλ(β)=

(4)

式中，N為觀測(cè)量的個(gè)數(shù)，β為待求量，xi和yi為觀測(cè)量，λ為懲罰因子，α為正則化參數(shù)，Pα(β)為懲罰項(xiàng)，其表示形式如下：

(5)

式中，p為待求參數(shù)的個(gè)數(shù)。α=1時(shí)，式(5)為L(zhǎng)asso算法模型；α=0時(shí)，式(5)為Ridge算法模型；α∈(0,1)時(shí)，式(5)為ElasticNet算法模型。

通過調(diào)整懲罰因子λ和正則化參數(shù)α，確定出剔除周期項(xiàng)的最優(yōu)參數(shù)值并對(duì)無效的周期項(xiàng)進(jìn)行剔除；然后再利用周期項(xiàng)取舍標(biāo)準(zhǔn)，剔除誤差超限的周期項(xiàng)，從而完成對(duì)周期項(xiàng)參數(shù)的篩選；最后通過平差對(duì)周期項(xiàng)進(jìn)行精確解算，計(jì)算出儀器周期項(xiàng)的振幅及初相，并對(duì)其結(jié)果進(jìn)行分析。

2.2 周期項(xiàng)取舍標(biāo)準(zhǔn)

根據(jù)誤差傳播定律，周期項(xiàng)的取舍標(biāo)準(zhǔn)為[2-3]：

(6)

式中，An為周期項(xiàng)的振幅，MA和Mφ分別為周期項(xiàng)的振幅中誤差和相位中誤差[4]：

(7)

式中，MX、MY分別為周期項(xiàng)的點(diǎn)位中誤差，A為對(duì)應(yīng)周期的振幅。

當(dāng)初相φ的求定誤差大于29°或者振幅誤差比小于2時(shí)，考慮剔除相應(yīng)的周期項(xiàng)。

3 實(shí)驗(yàn)分析

3.1 數(shù)據(jù)介紹

本文以2000國家重力基本網(wǎng)和2020最新國家重力基準(zhǔn)網(wǎng)中G796、G922兩臺(tái)LCR-G型相對(duì)重力儀觀測(cè)數(shù)據(jù)為例，對(duì)其周期項(xiàng)進(jìn)行參數(shù)解算，相應(yīng)的觀測(cè)量統(tǒng)計(jì)信息見表2。由表2可以看出，兩臺(tái)儀器的觀測(cè)量個(gè)數(shù)均比待求量個(gè)數(shù)多，保證了足夠的多余觀測(cè)量來解算相對(duì)重力儀的周期參數(shù)；最大段差、最小段差、平均段差分別在同一量級(jí)，可以保證所選觀測(cè)數(shù)據(jù)具有同一性。

表2 重力儀觀測(cè)量統(tǒng)計(jì)

3.2 周期項(xiàng)選取實(shí)驗(yàn)分析

首先，不考慮周期項(xiàng)對(duì)重力基準(zhǔn)網(wǎng)進(jìn)行平差，將得到的各個(gè)重力點(diǎn)上的重力值及所有儀器的格值一次項(xiàng)作為輸入數(shù)據(jù)，組建由G796、G922兩臺(tái)儀器周期項(xiàng)構(gòu)成的多元線性方程組，并采用正則化算法進(jìn)行解算，正則化參數(shù)α分別取0.1、0.5、1.0，懲罰因子λ取0～100，由解算結(jié)果計(jì)算的均方誤差如圖1所示。從圖1可以看出，當(dāng)懲罰因子λ大于80時(shí)，均方誤差均開始增大；正則化參數(shù)α的取值對(duì)周期項(xiàng)的均方誤差幾乎沒有影響。由此得出，在解算儀器周期項(xiàng)時(shí)，懲罰因子λ取80，正則化參數(shù)α的值可以任意選取。

圖1 重力儀周期項(xiàng)均方誤差Fig.1 The mean square error of periodicterms of gravity meters

采用正則化方法，懲罰因子λ取80，正則化參數(shù)α分別取0.1、0.5、1.0，解算的G796、G922儀器周期參數(shù)結(jié)果見圖2。可以看出，正則化方法對(duì)周期項(xiàng)的剔除效果明顯，正則化參數(shù)α的選取對(duì)周期項(xiàng)的篩選幾乎不產(chǎn)生影響。

圖2 重力儀周期項(xiàng)振幅Fig.2 The amplitude values of periodic terms of gravity meters

采用式(6)對(duì)由正則化方法解算的周期項(xiàng)再次進(jìn)行判定，去掉誤差超限的周期項(xiàng)，最終的篩選結(jié)果見表3。

表3 重力儀周期項(xiàng)參數(shù)篩選結(jié)果

3.3 結(jié)果分析

由于固定線性項(xiàng)和不固定線性項(xiàng)的標(biāo)定結(jié)果基本一致[10]，因此，采用不固定線性項(xiàng)的方式來求解標(biāo)定參數(shù)。將表3中篩選出的有效周期項(xiàng)確定的參數(shù)代入國家重力基準(zhǔn)網(wǎng)的誤差方程，重新進(jìn)行平差得到周期項(xiàng)的振幅和相位，結(jié)果見表4。由表4可知，通過正則化方法確定的周期項(xiàng)，振幅中誤差優(yōu)于4.2 μGal，相位中誤差優(yōu)于27.7，解算結(jié)果精度均在周期項(xiàng)取舍標(biāo)準(zhǔn)要求以內(nèi)。2000國家重力基本網(wǎng)G796重力儀周期項(xiàng)振幅最大可達(dá)12 μGal，G922重力儀周期項(xiàng)振幅最大可達(dá)8.0 μGal；2020最新國家重力基準(zhǔn)網(wǎng)G796重力儀周期項(xiàng)振幅最大可達(dá)10.7 μGal，G922重力儀周期項(xiàng)振幅最大可達(dá)7.5 μGal。兩臺(tái)儀器的振幅在近20 a內(nèi)均有所下降，可能是由于儀器出廠時(shí)間較長(zhǎng)、部件老化引起的；相同儀器在近20 a內(nèi)所確定的周期并不相同，分析認(rèn)為是相對(duì)重力儀解算的周期項(xiàng)的振幅和相位大小與采用的觀測(cè)數(shù)據(jù)及儀器觀測(cè)時(shí)所處的環(huán)境有關(guān)。

表4 重力儀周期項(xiàng)計(jì)算結(jié)果

考慮周期項(xiàng)前后重力基準(zhǔn)網(wǎng)平差重力值的變化見表5?？梢钥闯觯紤]周期項(xiàng)后，對(duì)2000國家重力基本網(wǎng)的平差重力值最大影響可達(dá)6.22 μGal，對(duì)2020最新國家重力基準(zhǔn)網(wǎng)的影響可達(dá)4.73 μGal。由此可見，周期項(xiàng)對(duì)重力基準(zhǔn)網(wǎng)影響的數(shù)量級(jí)達(dá)到μGal，不可忽視。

表5 考慮周期項(xiàng)前后重力值變化的統(tǒng)計(jì)

4 結(jié) 語

本文基于國家重力基準(zhǔn)網(wǎng)觀測(cè)數(shù)據(jù)，采用正則化方法及周期項(xiàng)取舍標(biāo)準(zhǔn)對(duì)G796、G922兩臺(tái)相對(duì)重力儀的周期項(xiàng)進(jìn)行篩選，并通過平差對(duì)周期項(xiàng)進(jìn)行精確解算，得出結(jié)論如下：1)使用正則化方法結(jié)合周期項(xiàng)取舍標(biāo)準(zhǔn)，可以有效篩選出帶有齒輪系統(tǒng)及量測(cè)螺旋系統(tǒng)的相對(duì)重力儀的周期項(xiàng)，并且通過實(shí)驗(yàn)確定，正則化方法懲罰因子λ取80最優(yōu)，正則化參數(shù)α的選取對(duì)周期項(xiàng)幾乎沒有影響；2)周期誤差對(duì)2000國家重力基本網(wǎng)和2020最新國家重力基準(zhǔn)網(wǎng)的重力值平差結(jié)果影響最大分別可達(dá)6.22 μGal和4.73 μGal，由此可見，周期項(xiàng)對(duì)重力基準(zhǔn)網(wǎng)影響的數(shù)量級(jí)達(dá)到μGal；3)相對(duì)重力儀的周期會(huì)隨著觀測(cè)環(huán)境或者時(shí)間變化，最終解算結(jié)果不一致，這是由于周期函數(shù)模型與實(shí)際周期變化有差異，模型只能近似而不能真正模擬周期的變化。