程一，李桐林，周帥

吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，長春 130026

0 引言

重力梯度測量相比重力測量具有更高的空間分辨率，可以提供更豐富的重力場信息，基于航空飛行平臺的重力梯度測量技術(shù)可快速、高效地完成面積性重力梯度數(shù)據(jù)采集工作，在礦產(chǎn)資源勘查、軍事目標(biāo)探測、慣性導(dǎo)航等諸多科學(xué)領(lǐng)域具有廣泛的應(yīng)用(Barnes and Barraud，2012；Cevallos et al.，2013；Dai et al.，2016；Rahimi and Naeeni，2017；馬國慶等，2020).而航空高動態(tài)條件下重力梯度測量相比傳統(tǒng)地面測量受到嚴(yán)重的載體平臺運(yùn)動的影響，航空飛行試驗(yàn)結(jié)果顯示動態(tài)測量環(huán)境中載體的殘余線運(yùn)動和角運(yùn)動是重力梯度動態(tài)測量噪聲中的主要誤差源(Roberts et al.，2016)，其中垂向運(yùn)動線加速度引起的上千厄缶(E)的重力梯度誤差將完全湮沒地質(zhì)體引起的重力梯度信號(楊曄和李達(dá)，2017).所以針對航空測量環(huán)境外部干擾源，應(yīng)明確各種外部干擾源的誤差傳遞機(jī)理，在此基礎(chǔ)上建立有效的動態(tài)測量誤差補(bǔ)償方法是重力梯度動態(tài)測量數(shù)據(jù)處理中的關(guān)鍵.

在航空重力梯度測量誤差補(bǔ)償方面，國外商用航空重力梯度儀FALCONTM、Air-FTG的生產(chǎn)報告描述了重力梯度儀實(shí)際測量中的數(shù)據(jù)處理流程，其中，事后誤差補(bǔ)償(Post-Mission Compensation，PMC)技術(shù)作為其中一項(xiàng)關(guān)鍵技術(shù)，由于商業(yè)保密性只進(jìn)行了簡單文字描述：“一系列復(fù)雜的算法對采集的數(shù)據(jù)進(jìn)行補(bǔ)償系數(shù)計算，而后利用它對采集的重力梯度數(shù)據(jù)進(jìn)行誤差補(bǔ)償”，但公開文獻(xiàn)未見詳細(xì)的實(shí)現(xiàn)算法細(xì)節(jié).由國外公開資料可以推斷出航空重力梯度動態(tài)測量過程中的誤差干擾是由確定的干擾源引起的，且有一定的規(guī)律，而不是隨機(jī)分布的，可采用數(shù)學(xué)方法通過補(bǔ)償飛行數(shù)據(jù)建立誤差補(bǔ)償模型實(shí)現(xiàn)測線數(shù)據(jù)的誤差補(bǔ)償.

在動態(tài)測量噪聲事后誤差補(bǔ)償方法建立之前需要明確誤差傳遞機(jī)理和模型，進(jìn)而建立有效的方法進(jìn)行誤差補(bǔ)償.王樹甫等(2000)給出了標(biāo)度因數(shù)不一致和加速度計安裝誤差的誤差傳遞模型，并提出了標(biāo)度因數(shù)反饋調(diào)整方法.羅嗣成(2007)給出了部分離心梯度和加速度計非線性系數(shù)的誤差模型并計算了各類誤差的臨界值.李海兵和蔡體菁(2010)給出了一種飛行載體自梯度誤差的模型并進(jìn)行了仿真試驗(yàn).韋宏瑋(2013)、袁園等(2014)基于數(shù)學(xué)仿真方法分析了單一系統(tǒng)誤差因素對儀器輸出的定量影響關(guān)系.Yu和Cai(2019)提出了一種考慮多種非理想因素的建模方法，并基于簡化后的誤差傳遞模型采用線性回歸進(jìn)行動態(tài)誤差補(bǔ)償.航空重力梯度儀動態(tài)測量誤差是由標(biāo)度因數(shù)不一致、加速度計位置安裝誤差、加速度計非線性系數(shù)、加速度計輸入軸安裝誤差等非理想因素與載體運(yùn)動參數(shù)相互耦合的結(jié)果，在誤差傳遞模型建立方面應(yīng)考慮多項(xiàng)因素間的相互作用.此外，各項(xiàng)非理想因素還隨著溫度、濕度、氣壓等外部環(huán)境因素變化且影響機(jī)制未知，基于傳統(tǒng)線性回歸分析的誤差補(bǔ)償方法很難有效地補(bǔ)償重力梯度儀輸出的動態(tài)測量誤差.

基于數(shù)據(jù)驅(qū)動的深度學(xué)習(xí)方法可直接從數(shù)據(jù)中學(xué)習(xí)隱含的規(guī)律和特征，針對外部各種干擾源引起的動態(tài)測量誤差，深度學(xué)習(xí)為該問題的解決提供了一種有效的策略.深度學(xué)習(xí)是機(jī)器學(xué)習(xí)的分支，通過構(gòu)建深層非線性網(wǎng)絡(luò)結(jié)構(gòu)和訓(xùn)練樣本集，自主學(xué)習(xí)隱含在輸入數(shù)據(jù)和輸出數(shù)據(jù)間的復(fù)雜特征信息和非線性函數(shù)關(guān)系.得益于計算機(jī)硬件的發(fā)展，尤其是圖形處理器(GPU)性能的提升，從2010年以來，深度學(xué)習(xí)成為了計算機(jī)科學(xué)的熱門研究領(lǐng)域，其中代表性的深度學(xué)習(xí)神經(jīng)網(wǎng)絡(luò)有深度卷積神經(jīng)網(wǎng)絡(luò)、深度置信網(wǎng)絡(luò)、循環(huán)神經(jīng)網(wǎng)絡(luò)等，深度學(xué)習(xí)網(wǎng)絡(luò)通過神經(jīng)元之間稀疏相連，極大地減少了需要訓(xùn)練的權(quán)值數(shù)和網(wǎng)絡(luò)計算量，在圖像去噪、語音識別、醫(yī)療診斷、地球物理數(shù)據(jù)處理等領(lǐng)域得到了廣泛應(yīng)用(Kononenko，2001；Zhang et al.，2017；Wu et al.，2019).同時，深度學(xué)習(xí)在其他領(lǐng)域中的誤差傳遞模型建立方面得到了一定的應(yīng)用.組合加工中心機(jī)床中的側(cè)隙誤差會導(dǎo)致零件幾何形狀的一系列變化，進(jìn)而影響設(shè)備的整體性能，但由于運(yùn)動副之間的機(jī)械磨損的不確定性，采用傳統(tǒng)解析誤差傳遞模型很難進(jìn)行側(cè)隙誤差的估計，Li等(2020)采用數(shù)據(jù)驅(qū)動的深度學(xué)習(xí)進(jìn)行側(cè)隙誤差的預(yù)測，深度置信網(wǎng)絡(luò)方法實(shí)現(xiàn)了隱式誤差傳遞模型的建立.另外，在渦輪機(jī)床的熱誤差補(bǔ)償方面，深度學(xué)習(xí)在復(fù)雜誤差傳遞模型建立中得到良好的應(yīng)用效果(Abdulshahed et al.，2015；Ma et al.，2017).而深度學(xué)習(xí)在地球物理測量數(shù)據(jù)動態(tài)噪聲壓制方面，Williams(1993)針對航空磁測過程中飛機(jī)姿態(tài)變化引起的動態(tài)噪聲干擾，基于神經(jīng)網(wǎng)絡(luò)算法進(jìn)行了高精度動態(tài)噪聲壓制；Yu等(2019)、Zhu等(2019)針對地震勘探數(shù)據(jù)中的時間相關(guān)噪聲基于深度學(xué)習(xí)神經(jīng)網(wǎng)絡(luò)建立干擾抑制方法，有效地提升了地震數(shù)據(jù)預(yù)處理的精度；許滔滔等(2020)采用循環(huán)神經(jīng)網(wǎng)絡(luò)深度學(xué)習(xí)方法對大地電磁測量數(shù)據(jù)中包含的時間相關(guān)外界人文設(shè)施電磁噪聲干擾進(jìn)行壓制，在實(shí)際數(shù)據(jù)應(yīng)用中得到良好的處理效果.

本文綜合考慮各項(xiàng)非理想因素與外部運(yùn)動參數(shù)間的耦合關(guān)系，基于旋轉(zhuǎn)式重力梯度儀的測量原理推導(dǎo)并建立了重力梯度儀誤差傳遞模型，在此基礎(chǔ)上建立重力梯度儀動態(tài)測量的仿真數(shù)據(jù)，采用基于深度學(xué)習(xí)網(wǎng)絡(luò)的事后誤差補(bǔ)償方法實(shí)現(xiàn)了重力梯度儀動態(tài)測量誤差的補(bǔ)償，利用數(shù)據(jù)驅(qū)動的誤差補(bǔ)償方法有效地壓制了動態(tài)測量噪聲，并在實(shí)際飛行測量數(shù)據(jù)的處理中進(jìn)一步驗(yàn)證了本文方法的有效性和實(shí)用性.

1 重力梯度儀動態(tài)測量誤差傳遞模型

1.1 理想情況下旋轉(zhuǎn)式重力梯度儀測量方程

圖1為旋轉(zhuǎn)加速度計式重力梯度儀(Gravity Gradient Instrument，GGI)工作原理圖.加速度計檢驗(yàn)質(zhì)量質(zhì)心與圓盤中心距離為R，加速度計i(i=1,2,3,4)檢驗(yàn)質(zhì)量質(zhì)心與原點(diǎn)的連線與X正半軸夾角為θi，圓盤以角速度ω逆時針方向旋轉(zhuǎn).圖1中測量坐標(biāo)系不隨圓盤旋轉(zhuǎn)，且選擇以O(shè)為原點(diǎn)的慣性坐標(biāo)系作為慣性參考系.

圖1 旋轉(zhuǎn)加速度計式重力梯度儀工作原理圖Fig.1 Working principle diagram of rotating accelerometer gravity gradiometer

基于旋轉(zhuǎn)加速度計式重力梯度儀工作原理，通過計算加計位置矢量Ri、測量坐標(biāo)系下加計處加速度ai、加計輸入軸方向加速度aIi、加計輸出Ui便可得到梯度儀輸出Eout.設(shè)加速度計i的位置矢量Ri為(Xi，Yi，Zi)T(i=1,2,3,4)，則Ri的各分量表達(dá)式為：

(1)

設(shè)原點(diǎn)O處的加速度a0為(aX0，aY0，aZ0)T，加速度計i處加速度ai為(aXi，aYi，aZi)T，對ai在原點(diǎn)O處進(jìn)行泰勒一階展開：

ai=a0+?！i，

(2)

其中，Γ為重力梯度張量矩陣.加速度計i輸入軸方向單位向量eIi為(sinθi,-cosθi,0).ai在輸入軸方向的投影為加速度計敏感到的加速度，加速度計i輸入軸方向加速度aIi為：

aIi=ai·eIi=aXisinθi-aYicosθi

-RΓXYcos2θi.

(3)

加速度計也會敏感到其輸入軸方向的旋轉(zhuǎn)加速度.一般認(rèn)為GGI圓盤旋轉(zhuǎn)速度恒定，故該旋轉(zhuǎn)加速度為測量坐標(biāo)系相對慣性系統(tǒng)繞Z軸變速旋轉(zhuǎn)的分量.設(shè)加速度計檢驗(yàn)質(zhì)量質(zhì)心位置處該旋轉(zhuǎn)角加速度為aωZ，則aIi變?yōu)椋?/p>

-RΓXYcos2θi-Raω Z.

(4)

將各加速度計輸入軸方向的加速度乘以其標(biāo)度因數(shù)，即可得到各加速度計的輸出.理想情況下，各加速度計標(biāo)度因數(shù)一致，設(shè)為k，則加速度計i的輸出Ui為：

×sin2θi-kRΓXYcos2θi-kRaω Z.

(5)

設(shè)t=0時刻加速度計1檢驗(yàn)質(zhì)量質(zhì)心位于x軸正半軸上，則θi=ωt+(i-1)×π/2.通過相對加速度計的輸出求和可去除共模加速度aX0、aY0的影響，即：

U1+U3=kR(ΓXX-ΓYY)sin2ωt-2kRΓXYcos2ωt

-2kRaω Z，

(6)

U2+U4=kR(-ΓXX+ΓYY)sin2ωt

+2kRΓXYcos2ωt-2kRaω Z.

(7)

對兩組相對加速度計的輸出和做差運(yùn)算，可去除旋轉(zhuǎn)角加速度aω Z的影響，重力梯度儀理想情況下的測量方程為：

Eout=(U1+U3)-(U2+U4)

=2kR(ΓXX-ΓYY)sin2ωt-4kRΓXYcos2ωt.(8)

通過對重力梯度儀輸出信號的解調(diào)處理便可獲取重力梯度分量ΓXX-ΓYY、ΓXY.

1.2 旋轉(zhuǎn)式重力梯度儀誤差傳遞模型

重力梯度儀測量誤差為系統(tǒng)各種非理想因素與動態(tài)測量外部干擾源(線運(yùn)動、角運(yùn)動等)相互耦合的結(jié)果，明確各種因素對重力梯度儀輸出的影響機(jī)制是進(jìn)行高精度事后誤差補(bǔ)償?shù)年P(guān)鍵.在理想情況下旋轉(zhuǎn)式重力梯度儀測量方程的基礎(chǔ)上，圖2為考慮各種非理想因素耦合下重力梯度儀誤差傳遞模型.本文將介紹各種非理想因素并給出其對重力梯度儀輸出的影響機(jī)制.類似上一節(jié)通過逐步推導(dǎo)各種非理想因素影響下的Ri、ai、實(shí)際加計坐標(biāo)系下加計處加速度a′i和加計輸出Ui便可完成綜合各種非理想因素的誤差傳遞模型的構(gòu)建.綜合考慮各類非理想因素時該方程過于復(fù)雜，文中不給出最終的測量誤差傳遞表達(dá)式.

圖2 旋轉(zhuǎn)加速度計式重力梯度儀誤差傳遞模型Fig.2 Error transfer mechanism of rotating accelerometer gravity gradiometer

由于機(jī)械制造工藝的限制，加速度計實(shí)際位置與理想位置存在偏差，即存在加速度計位置安裝誤差，將對位置矢量Ri產(chǎn)生影響.加速度計位置安裝誤差包括質(zhì)心垂向誤差zi，質(zhì)心徑向誤差Ti，質(zhì)心切向誤差ri三個部分，分別表示加速度計i檢驗(yàn)質(zhì)量質(zhì)心實(shí)際位置與理想位置的位移在垂直于GGI圓盤方向、圓盤平面內(nèi)垂直于圓盤切線方向、圓盤平面內(nèi)平行于圓盤切線方向的投影，分別以垂直圓盤向上、由圓盤中心向外、順時針方向?yàn)檎较?Ri各分量表達(dá)式將變?yōu)椋?/p>

(9)

由于穩(wěn)定平臺無法完全隔離載體平臺的角運(yùn)動，將會有殘余的角運(yùn)動影響測量坐標(biāo)系下加計處加速度ai.殘余角運(yùn)動通過產(chǎn)生的離心加速度和旋轉(zhuǎn)加速度影響ai，ai將變?yōu)椋?/p>

ai=a0+?！i-aCi+aRi=a0+?！i

-ω×(ω×Ri)+aω×Ri，

(10)

其中，aCi為加速度計處的離心加速度矢量；aRi為加速度計處的旋轉(zhuǎn)加速度矢量；ω為旋轉(zhuǎn)角速度矢量；aω為旋轉(zhuǎn)角加速度矢量.

設(shè)以加速度計輸入軸、輸出軸、擺軸為坐標(biāo)軸的坐標(biāo)系為加計坐標(biāo)系.其中，輸出軸方向理想情況下在GGI圓盤平面內(nèi)，擺軸方向理想情況下與測量坐標(biāo)系Z軸同向.

由于加速度計安裝在GGI圓盤上時，其輸入軸方向往往與理想方向有一定的偏差，即存在加速度計輸入軸安裝誤差，影響實(shí)際加計坐標(biāo)系下加計處加速度a′i.輸入軸安裝誤差包括加速度計i輸入軸理想方向繞其輸入軸旋轉(zhuǎn)角度φIi、繞輸出軸旋轉(zhuǎn)角度φOi和繞擺軸旋轉(zhuǎn)角度φPi.通過坐標(biāo)變換可得到a′i：

a′i=C1i·C2i·C3i·a′i0=C1i·C2i·C3i·Cmai·ai，

(11)

其中，C1i、C2i、C3i為加速度計i理想加計坐標(biāo)系到實(shí)際加計坐標(biāo)系的坐標(biāo)變換矩陣；a′i0為理想加計坐標(biāo)系下加計處加速度矢量；Cmai為測量坐標(biāo)系到理想加計坐標(biāo)系的坐標(biāo)變換矩陣.即：

(12)

其中，aIi、aOi、aPi為a′i各分量；aIi0、aOi0、aPi0為a′i0各分量.a′i0各分量為：

aIi0=aXisinθi-aYicosθi，

(13)

aOi0=aXicosθi+aYisinθi，

(14)

aPi0=aZi.

(15)

由于制造工藝的精度問題，加速度計的標(biāo)度因數(shù)存在不一致性，且在實(shí)際測量中，加速度計并不能只線性敏感其輸入軸方向的加速度，輸出還與一些非線性系數(shù)有關(guān)，以上兩種非理想因素都會影響加速度計輸出Ui.標(biāo)度因數(shù)不一致時，加速度計間的差分運(yùn)算將無法再完全消除平動、旋轉(zhuǎn)產(chǎn)生的共模量.加速度計非線性系數(shù)則會使得輸入軸、輸出軸和擺軸三方向的加速度及其高階項(xiàng)、耦合項(xiàng)都影響加速度計的輸出.Ui將變?yōu)椋?/p>

+KIPiaIiaPi+KI OiaIiaO i+KOPiaO iaPi)，

(16)

其中，ki為標(biāo)度因數(shù)；K0i為零位偏置；KIPi、KI O i、KOPi為交叉耦合系數(shù)；KIIi、KO Oi、KPPi為二階系數(shù).

最后，通過加速度計輸出間的差分可以得到重力梯度儀輸出Eout，完成誤差傳遞模型的構(gòu)建：

Eout=(U1+U3)-(U2+U4).

(17)

2 事后誤差補(bǔ)償方法

重力梯度儀動態(tài)測量噪聲輸出是由角運(yùn)動、線運(yùn)動等外部干擾源與系統(tǒng)多種非理想因素相互耦合的結(jié)果，與外部各干擾源之間具有非線性、復(fù)雜的表達(dá)關(guān)系，并且外部環(huán)境參數(shù)(溫度、濕度、氣壓、磁場等)對加速度計標(biāo)度因數(shù)、安裝誤差等系統(tǒng)誤差項(xiàng)的影響無解析的定量關(guān)系，很難基于傳統(tǒng)線性回歸的方式進(jìn)行誤差補(bǔ)償.針對外部干擾源與重力梯度儀動態(tài)測量噪聲輸出間的復(fù)雜關(guān)系，深層次特征學(xué)習(xí)的深度學(xué)習(xí)方法為該問題的解決提供了一種有效的策略，通過構(gòu)建深層非線性網(wǎng)絡(luò)結(jié)構(gòu)和大規(guī)模訓(xùn)練樣本集，可以自主學(xué)習(xí)隱含在輸入數(shù)據(jù)和輸出數(shù)據(jù)間的復(fù)雜特征信息和非線性函數(shù)關(guān)系.本文選擇航空重力梯度動態(tài)測量過程中獲取的三方向加速度、角速度、角加速度、加計相位以及重力梯度儀輸出噪聲數(shù)據(jù)建立樣本集，通過搭建深度學(xué)習(xí)神經(jīng)網(wǎng)絡(luò)建立動態(tài)測量誤差與運(yùn)動參數(shù)之間的關(guān)系，進(jìn)而實(shí)現(xiàn)重力梯度儀動態(tài)測量數(shù)據(jù)的事后誤差補(bǔ)償.

本文選擇多層感知機(jī)網(wǎng)絡(luò)進(jìn)行事后誤差補(bǔ)償方法的建立，包括輸入層、隱藏層和輸出層，其中輸入層負(fù)責(zé)接收、輸入外來數(shù)據(jù)；隱藏層負(fù)責(zé)計算、加工輸入層接收的數(shù)據(jù)，并逐層傳遞向輸出層，隱藏層的存在使得神經(jīng)網(wǎng)絡(luò)具備了表述復(fù)雜問題的能力；輸出層負(fù)責(zé)接收、計算前一層傳遞來的數(shù)據(jù)并輸出最后的結(jié)果.如圖3為建立的旋轉(zhuǎn)式重力梯度儀誤差補(bǔ)償?shù)亩鄬痈兄獧C(jī)網(wǎng)絡(luò)模型，其中M-P感知機(jī)是多層感知機(jī)網(wǎng)絡(luò)最基本的模型，M-P感知機(jī)的輸出值y為：

y(x)=f(WTφ(x)+b)，

(18)

其中，W為權(quán)重參數(shù)向量；φ為任意基本函數(shù)；x為輸入數(shù)據(jù)向量；b為偏置參數(shù)；f為激活函數(shù).激活函數(shù)的作用是給線性函數(shù)加權(quán)結(jié)構(gòu)引入非線性因素，使得神經(jīng)網(wǎng)絡(luò)能夠擬合非線性函數(shù).在基于深度學(xué)習(xí)的事后誤差補(bǔ)償方法的網(wǎng)絡(luò)模型建立中，本文設(shè)計的神經(jīng)網(wǎng)絡(luò)模型包括一個輸入層、四個隱藏層和一個輸出層，各層之間采用全連接模式進(jìn)行連接，各層的節(jié)點(diǎn)數(shù)分別為128、64、32、16、8、1，并選擇“Adam”優(yōu)化算法進(jìn)行深度學(xué)習(xí)網(wǎng)絡(luò)模型的訓(xùn)練(圖3).具體的實(shí)現(xiàn)流程如下：

圖3 動態(tài)測量誤差補(bǔ)償深度學(xué)習(xí)網(wǎng)絡(luò)模型Fig.3 Deep learning network for dynamic measurement error compensation

(1)將重力梯度儀的運(yùn)動參數(shù)：角速度(ωX、ωY、ωZ)、線加速度(aX0、aY0、aZ0)、角加速度(aωX、aωY、aωZ)及加速度計實(shí)時相位θ1作為輸入數(shù)據(jù)，梯度儀輸出噪聲信號Eout作為輸出數(shù)據(jù)建立數(shù)據(jù)樣本集，并劃分為訓(xùn)練集和測試集.

(2)搭建圖3所示的神經(jīng)網(wǎng)絡(luò)模型，將輸入映射到輸出.

(3)配置學(xué)習(xí)過程：選擇優(yōu)化器、損失函數(shù)及監(jiān)控指標(biāo)等參數(shù).

(4)將數(shù)據(jù)集進(jìn)行均值方差歸一化，保存歸一化參數(shù)，迭代訓(xùn)練神經(jīng)網(wǎng)絡(luò).

(5)選擇預(yù)計進(jìn)行誤差補(bǔ)償?shù)奶荻葍x輸出對應(yīng)的運(yùn)動參數(shù)及加速度計實(shí)時相位數(shù)據(jù)，讀取保存的歸一化參數(shù)進(jìn)行歸一化，輸入到訓(xùn)練得到的神經(jīng)網(wǎng)絡(luò)中，將輸出進(jìn)行反歸一化得到預(yù)測的梯度儀輸出噪聲信號，進(jìn)而完成動態(tài)測量數(shù)據(jù)事后誤差補(bǔ)償.

多層感知機(jī)網(wǎng)絡(luò)在算法上實(shí)現(xiàn)可分為三步：參數(shù)初始化，正向傳播計算輸出結(jié)果，反向傳播更新網(wǎng)絡(luò)參數(shù).正向傳播即是將數(shù)據(jù)輸入輸入層，經(jīng)由各隱藏層逐層計算、加工，之后推進(jìn)到輸出層輸出結(jié)果的算法.神經(jīng)網(wǎng)絡(luò)反向傳播算法的作用是更新神經(jīng)網(wǎng)絡(luò)內(nèi)的參數(shù)，通過定義一個損失函數(shù)來衡量神經(jīng)網(wǎng)絡(luò)真實(shí)輸出值與期望輸出值(標(biāo)簽值)的差距，其后對損失函數(shù)求導(dǎo)并結(jié)合鏈?zhǔn)椒▌t迭代計算網(wǎng)絡(luò)中每一層參數(shù)的梯度值，最后在當(dāng)前參數(shù)中減去人為設(shè)定的學(xué)習(xí)率與該梯度值的積即可完成參數(shù)的更新.損失函數(shù)選擇L2范數(shù)定義的方式進(jìn)行計算：

(19)

其中，yi表示已知的重力梯度儀輸出誤差；yip表示預(yù)測的重力梯度儀輸出誤差.

3 基于誤差傳遞模型的仿真數(shù)據(jù)試驗(yàn)

基于建立的誤差傳遞模型，按照重力梯度儀實(shí)際工作時狀態(tài)對各種非理想因素進(jìn)行設(shè)置，具體取值如表1所示.此外，加速度計檢驗(yàn)質(zhì)量質(zhì)心到GGI圓盤中心的距離R為0.2 m，重力梯度儀采樣頻率為800 Hz，GGI圓盤旋轉(zhuǎn)頻率f為0.25 Hz，重力梯度儀三方向加速度aX0、aY0、aZ0設(shè)置為由-0.3～0.3 m·s-2內(nèi)隨機(jī)確定的加速度和在0.01～1 Hz均勻分布的100個振幅在0.006～0.015 m·s-2內(nèi)隨機(jī)確定、相位隨機(jī)確定的分量組成.重力梯度儀三方向角速度ωX、ωY、ωZ設(shè)置為由-5～5arcsec·s-1內(nèi)隨機(jī)確定的角速度和在0.01 Hz到1 Hz均勻分布的100個振幅在0.2～0.5arcsec·s-1內(nèi)隨機(jī)確定、相位隨機(jī)確定的分量組成.上述隨機(jī)數(shù)均服從連續(xù)均勻分布，確定后各量不隨時間變化.重力梯度儀三方向旋轉(zhuǎn)角加速度aω X、aω Y、aω Z通過ωX、ωY、ωZ的一階差分獲得.解調(diào)低通濾波窗選擇為矩形窗，截止頻率設(shè)為0.075 Hz.

表1 重力梯度儀系統(tǒng)非理想因素參數(shù)Table 1 Non-ideal factor parameters of gravity gradiometer system

通過對重力梯度儀系統(tǒng)非理想因素和外部運(yùn)動參數(shù)的設(shè)計完成重力梯度儀測量輸出數(shù)據(jù)，并將仿真數(shù)據(jù)集隨機(jī)劃分為訓(xùn)練集和測試集兩部分，訓(xùn)練集用于深度學(xué)習(xí)網(wǎng)絡(luò)的訓(xùn)練，測試集用于建立深度學(xué)習(xí)網(wǎng)絡(luò)的泛化性驗(yàn)證，采用建立的深度學(xué)習(xí)網(wǎng)絡(luò)模型進(jìn)行5輪的迭代計算，訓(xùn)練集和測試集每次迭代的歸一化損失函數(shù)如圖4所示，經(jīng)過前兩次迭代，損失函數(shù)有了較為明顯的下降，其后趨于平穩(wěn).

圖4 訓(xùn)練集與測試集的損失函數(shù)值Fig.4 Loss function values of training set and testing set

仿真數(shù)據(jù)誤差補(bǔ)償結(jié)果如圖5所示，其中，ΓXX-ΓYY梯度分量通過正弦二倍頻解調(diào)得到，ΓXY梯度分量通過余弦二倍頻解調(diào)得到.利用訓(xùn)練集數(shù)據(jù)建立的深度學(xué)習(xí)網(wǎng)絡(luò)模型對測試集數(shù)據(jù)進(jìn)行應(yīng)用，分析事后誤差補(bǔ)償方法的泛化性，通過重力梯度儀輸出數(shù)據(jù)誤差補(bǔ)償圖可以看出，誤差補(bǔ)償后剩余的ΓXX-ΓYY、ΓXY梯度分量噪聲水平下降近兩個數(shù)量級，有效地壓制了動態(tài)測量噪聲.其中，ΓXX-ΓYY梯度分量誤差補(bǔ)償前后的均方根誤差分別為1111 E、15 E，ΓXY梯度分量誤差補(bǔ)償前后的均方根誤差分別為681 E、9 E，ΓXX-ΓYY與ΓXY梯度分量補(bǔ)償后均方根誤差與補(bǔ)償前均方根誤差比值分別為1.4%、1.3%，建立的事后誤差補(bǔ)償方法的噪聲壓制程度在98.5%以上，說明現(xiàn)有基于數(shù)據(jù)驅(qū)動的深度學(xué)習(xí)方法能夠有效地進(jìn)行高動態(tài)測量環(huán)境中的重力梯度動態(tài)測量誤差補(bǔ)償，但補(bǔ)償結(jié)果中仍殘留部分噪聲干擾，這主要與深度學(xué)習(xí)網(wǎng)絡(luò)模型結(jié)構(gòu)設(shè)計、參數(shù)選取等有關(guān)，改善深度學(xué)習(xí)方法有望進(jìn)一步提高事后誤差補(bǔ)償精度.

圖5 仿真數(shù)據(jù)誤差補(bǔ)償結(jié)果圖(a,b)誤差補(bǔ)償前后的重力梯度儀輸出數(shù)據(jù)；(c,d)誤差補(bǔ)償前后的ΓXX-ΓYY梯度分量；(e,f)誤差補(bǔ)償前后的ΓXY梯度分量.Fig.5 Error compensation result of synthetic data(a,b)Output data of gravity gradiometer before and after error compensation;(c,d)ΓXX-ΓYY gravity gradient component before and after error compensation;(e,f)ΓXY gravity gradient component before and after error compensation.

4 實(shí)測航空重力梯度數(shù)據(jù)應(yīng)用

為了驗(yàn)證建立的事后誤差補(bǔ)償在實(shí)際數(shù)據(jù)中的有效性，對航空重力梯度儀實(shí)際飛行數(shù)據(jù)進(jìn)行方法測試.采用深度學(xué)習(xí)策略進(jìn)行動態(tài)測量數(shù)據(jù)事后誤差補(bǔ)償?shù)年P(guān)鍵在于特征全面的訓(xùn)練樣本集的建立，為了獲得僅包含重力梯度動態(tài)測量噪聲的樣本集，在試驗(yàn)區(qū)進(jìn)行高空(3600 m)方形重復(fù)線飛行，如圖6黑色邊框圈定區(qū)域?yàn)閷?shí)際航空飛行路徑位置.本文收集了該區(qū)域的航空重力測量數(shù)據(jù)，采用波數(shù)域?qū)?shù)計算方法進(jìn)行了重力梯度數(shù)據(jù)的轉(zhuǎn)換計算(周帥等，2016)，結(jié)果顯示在高空飛行路徑位置的ΓXX-ΓYY、ΓXY重力梯度分量異常變化在10 E以內(nèi)，可認(rèn)為在該飛行高度獲取的重力梯度儀輸出數(shù)據(jù)不包含由地質(zhì)體引起的重力梯度信號，重力梯度儀的原始輸出均為動態(tài)測量噪聲(圖6).通過載體平臺的慣性測量單元可測得重力梯度儀的加速度、角速度等運(yùn)動參數(shù)，采用加裝的光柵編碼器可測得加速度計的實(shí)時相位，進(jìn)而結(jié)合重力梯度儀輸出數(shù)據(jù)建立深度學(xué)習(xí)樣本集.為了測試建立的深度學(xué)習(xí)網(wǎng)絡(luò)模型的泛化性，以4：1的比例將該實(shí)測數(shù)據(jù)集隨機(jī)劃分為訓(xùn)練集和測試集兩部分，其中訓(xùn)練集用于誤差補(bǔ)償深度學(xué)習(xí)模型的建立，測試集用于模型泛化性驗(yàn)證.

圖6 試驗(yàn)區(qū)換算重力梯度數(shù)據(jù)(a)ΓXX-ΓYY梯度分量；(b)ΓXY梯度分量.Fig.6 The calculated gravity gradient data(a)ΓXX-ΓYY gravity gradient component;(b)ΓXY gravity gradient component.

訓(xùn)練共計進(jìn)行10輪迭代，10輪迭代中訓(xùn)練集及測試集的歸一化損失函數(shù)如圖7a所示，訓(xùn)練集和測試集損失函數(shù)都是明顯下降趨勢且維持在較低水平，深度學(xué)習(xí)網(wǎng)絡(luò)模型可以有效地提取出數(shù)據(jù)特征.由圖7b—d重力梯度儀原始輸出數(shù)據(jù)補(bǔ)償效果圖可以看出，該方法能夠有效地對動態(tài)測量噪聲進(jìn)行抑制，提高了重力梯度儀動態(tài)測量精度.對補(bǔ)償前后的重力梯度儀輸出數(shù)據(jù)進(jìn)行解調(diào)，ΓXX-ΓYY梯度分量采用正弦二倍頻進(jìn)行解調(diào)，ΓXY梯度分量采用余弦二倍頻進(jìn)行解調(diào)，獲得如圖8所示補(bǔ)償前后的ΓXX-ΓYY和ΓXY重力梯度分量.補(bǔ)償后兩組梯度分量的噪聲水平相比補(bǔ)償前抑制程度超過一個數(shù)量級.由補(bǔ)償前后的均方根誤差進(jìn)行補(bǔ)償效果的定量評價，ΓXX-ΓYY梯度分量補(bǔ)償前后的均方根誤差分別為2806 E、56 E，ΓXY梯度分量補(bǔ)償前后的均方根誤差分別為1459 E、24 E，ΓXX-ΓYY與ΓXY梯度分量補(bǔ)償后均方根誤差與補(bǔ)償前均方根誤差比值分別為1.9%、1.6%，建立的事后誤差補(bǔ)償方法的噪聲壓制程度在98%以上，結(jié)果顯示基于深度學(xué)習(xí)的事后補(bǔ)償方法能夠有效地進(jìn)行實(shí)際數(shù)據(jù)的噪聲抑制處理.基于實(shí)際數(shù)據(jù)的誤差補(bǔ)償結(jié)果相比仿真數(shù)據(jù)而言具有較低的噪聲壓制能力，主要由于實(shí)際飛行測量過程中溫度、濕度、氣壓等外部環(huán)境對系統(tǒng)非理想因素的影響，間接影響了重力梯度儀輸出，而現(xiàn)有深度學(xué)習(xí)網(wǎng)絡(luò)模型訓(xùn)練僅考慮了外部運(yùn)動參數(shù)，未考慮外部環(huán)境參數(shù).

圖7 重力梯度儀原始輸出數(shù)據(jù)補(bǔ)償結(jié)果圖(a)訓(xùn)練集與測試集的損失函數(shù)值；(b)梯度儀原始輸出數(shù)據(jù)；(c)預(yù)測的梯度儀輸出噪聲；(d)誤差補(bǔ)償后的梯度儀輸出數(shù)據(jù).Fig.7 Compensation result of gravity gradiometer original output data(a)Loss function values of training set and testing set;(b)Original output data of gravity gradiometer;(c)Predicted output noise of gravity gradiometer;(d)Gravity gradiometer output data after error compensation.

圖8 實(shí)測數(shù)據(jù)誤差補(bǔ)償結(jié)果圖(a,b)補(bǔ)償前后的ΓXX-ΓYY梯度分量；(c,d)補(bǔ)償前后的ΓXY梯度分量.Fig.8 Compensation result of real measured data(a,b)ΓXX-ΓYY gravity gradient component before and after error compensation;(c,d)ΓXY gravity gradient component before and after error compensation.

5 結(jié)論

通過建立的重力梯度儀誤差傳遞模型可以看出，動態(tài)測量過程中梯度儀輸出噪聲主要是由外部線運(yùn)動、角運(yùn)動等參數(shù)與系統(tǒng)各項(xiàng)非理想因素耦合激勵的結(jié)果，它們之間具有復(fù)雜、隱式、非線性的關(guān)系，采用傳統(tǒng)線性回歸方法進(jìn)行誤差補(bǔ)償存在一定局限性，而基于數(shù)據(jù)驅(qū)動的深度學(xué)習(xí)方法為該問題求解提供了一種有效地策略.通過仿真數(shù)據(jù)測試顯示本文建立事后誤差補(bǔ)償方法壓制的動態(tài)測量噪聲水平近兩個數(shù)量級，該方法對動態(tài)測量噪聲壓制程度高于98.5%，說明了深度學(xué)習(xí)方法在解決重力梯度儀誤差補(bǔ)償問題的有效性，優(yōu)化深度學(xué)習(xí)網(wǎng)絡(luò)模型結(jié)構(gòu)設(shè)計、參數(shù)選取、優(yōu)化算法選擇等，有望進(jìn)一步提高深度學(xué)習(xí)網(wǎng)絡(luò)模型的誤差補(bǔ)償精度.在實(shí)際航空重力梯度測量數(shù)據(jù)的應(yīng)用中，補(bǔ)償結(jié)果顯示動態(tài)測量噪聲壓制程度高于98%，驗(yàn)證了該方法實(shí)際數(shù)據(jù)處理的實(shí)用性和有效性，為了進(jìn)一步提高實(shí)際航空重力梯度測量精度，應(yīng)考慮加入外部測量環(huán)境參數(shù)(磁場、溫度、濕度、氣壓等)作為樣本集輸入變量，進(jìn)而建立包含特征更加完備的數(shù)據(jù)樣本集.