周波陽(yáng) 羅志才 鐘 波 鄭 凱 魏艷平

1 廣東工業(yè)大學(xué)測(cè)繪工程系，廣州市大學(xué)城外環(huán)西路100號(hào)，510006

2 武漢大學(xué)測(cè)繪學(xué)院，武漢市珞喻路129號(hào)，430079

3 武漢大學(xué)地球空間環(huán)境與大地測(cè)量教育部重點(diǎn)實(shí)驗(yàn)室，武漢市珞喻路129號(hào)，430079

4 中國(guó)石油集團(tuán)工程設(shè)計(jì)有限責(zé)任公司華北分公司，任丘市建設(shè)路，062552

對(duì)機(jī)載應(yīng)用來(lái)說(shuō)，確定載體加速度的3種基本方法是位置求導(dǎo)法、相位差分法、多普勒頻移法［1－4］，這3種方法都需要通過(guò)差分運(yùn)算才能得出載體的加速度。牛頓中心差分器是運(yùn)用最為廣泛的一類(lèi)差分器，其關(guān)鍵在于點(diǎn)數(shù)的選取。文獻(xiàn)［5］分析了利用一階牛頓中心差分器進(jìn)行GPS定速的主要誤差來(lái)源，認(rèn)為一階牛頓中心差分器的最佳點(diǎn)數(shù)應(yīng)使截?cái)嗾`差與觀測(cè)誤差之和最小。這一結(jié)論是在全頻譜范圍內(nèi)進(jìn)行考慮的，適合于航空攝影測(cè)量、GPS／INS組合導(dǎo)航等領(lǐng)域。航空重力測(cè)量最終獲取的是帶限重力信號(hào)，需要將載體加速度從比力觀測(cè)值中分離出來(lái)。由于重力信號(hào)表現(xiàn)為長(zhǎng)波特性，在數(shù)據(jù)處理中需對(duì)原始觀測(cè)值進(jìn)行低通濾波處理，因此對(duì)航空重力測(cè)量而言，傳感器原始觀測(cè)值的高頻部分都可認(rèn)為是噪聲，載體加速度只需在低頻端保持足夠的精度即可。考慮到這一情況，在航空重力測(cè)量的數(shù)據(jù)處理中，在全頻譜范圍內(nèi)考慮差分器的精度可能并不合適。本文將從頻譜的角度分析牛頓中心差分器的點(diǎn)數(shù)對(duì)航空重力測(cè)量中載體加速度的影響，并采用實(shí)測(cè)的機(jī)載GPS數(shù)據(jù)進(jìn)行驗(yàn)證。

1 理想差分器與理想低通差分器

利用載體的位置信息通過(guò)差分方式可以計(jì)算出載體的速度和加速度。為簡(jiǎn)化問(wèn)題，本文僅以垂直方向?yàn)槔?。由微分的定義可知：

對(duì)式（2）取極限，得：

因此，理想差分器的頻率響應(yīng)為：

式中，ω＝2πf為圓頻率。由式（4）可知，理想差分器的頻率特性在0～π內(nèi)是線(xiàn)性增長(zhǎng)的。和信號(hào)相比，噪聲一般在高頻端，而差分過(guò)程會(huì)放大噪聲，特別是高頻噪聲。這一問(wèn)題的解決辦法是采用低通差分器，理想低通差分濾波器的頻率響應(yīng)為：

式中，απ為理想低通差分器的截止頻率，頻率歸一化后則稱(chēng)α為截止頻率。

2 牛頓中心差分器

理想低通差分器在物理上無(wú)法實(shí)現(xiàn)，人們通常采用頻率響應(yīng)函數(shù)H（ejω）來(lái)逼近H′d（ejω）。此時(shí)，H（ejω）一般具有以下形式：

相應(yīng)的轉(zhuǎn)移函數(shù)為：

若差分器的輸入為x（n），輸出為y（n），則有：

表1 一階牛頓中心差分器的系數(shù)Tab.1 Coefficients of one－order Newton central differentiator

依據(jù)式（4）、（10），可給出理想差分器與牛頓中心差分器的幅頻響應(yīng)，如圖1所示。從圖中可以看出，牛頓中心差分器是對(duì)理想差分器的逼近，點(diǎn)數(shù)越多，對(duì)理想差分器的逼近效果越好，但同時(shí)差分器的低通效果在理論上也越差。文獻(xiàn)［7］指出，低通差分器對(duì)理想差分器的逼近也可采用截止頻率α來(lái)描述：απ為曲線(xiàn)｜H（ejω）｜與直線(xiàn)的交點(diǎn)的橫坐標(biāo)，如圖1所示。表2給出了不同點(diǎn)數(shù)牛頓中心差分器的截止頻率αi（i＝3，5，7，9，11），可知α3＜α5＜α7＜α9＜α11。當(dāng)信號(hào)的有效頻段位于低頻段，且信號(hào)的截止頻率α?α3時(shí)，從理論上可以得出兩條結(jié)論：1）為避免高頻噪聲的影響，應(yīng)對(duì)差分結(jié)果進(jìn)行低通濾波；2）由于不同點(diǎn)數(shù)的牛頓中心差分器的幅頻響應(yīng)在低頻端完全吻合，低通濾波后不同點(diǎn)數(shù)的牛頓中心差分器的輸出結(jié)果之間的差異應(yīng)比較小。

圖1 理想差分器及牛頓中心差分器的幅頻特性Fig.1 Magnitude－frequency characteristic of ideal differentiator and Newton central differentiators

表2 一階牛頓中心差分器的截止頻率Tab.2 Cutoff frequencies of one－order Newton central differentiators

3 數(shù)值實(shí)驗(yàn)與分析

3.1 靜態(tài)實(shí)驗(yàn)

選取美國(guó)某地區(qū)2個(gè)IGS跟蹤站的GPS觀測(cè)數(shù)據(jù)，測(cè)站編號(hào)分別為A和B，基線(xiàn)長(zhǎng)度為30 km，數(shù)據(jù)采樣間隔為1s，觀測(cè)時(shí)間為2h。設(shè)測(cè)站A的坐標(biāo)固定，將B的坐標(biāo)觀測(cè)值與A的固定坐標(biāo)相減，可得到時(shí)間序列→ZAB，對(duì)→ZAB作二次牛頓中心差分后可得到測(cè)站B相對(duì)于測(cè)站A的加速度。由于差分過(guò)程會(huì)放大噪聲，必須對(duì)差分結(jié)果作低通濾波處理。本文選用文獻(xiàn)［9］給出的迭代高斯低通濾波器對(duì)加速度進(jìn)行濾波，若載體的飛行速度為155 m／s，所需重力數(shù)據(jù)的半波長(zhǎng)分辨率為10km，則高斯濾波器的長(zhǎng)度為65s，對(duì)應(yīng)的截止頻率為α＝0.015Hz。在短時(shí)間內(nèi)可認(rèn)為測(cè)站B相對(duì)加速度的真值為零，因此低通濾波后所得加速度的標(biāo)準(zhǔn)差即可反映數(shù)值處理的絕對(duì)精度，表3給出了相關(guān)的統(tǒng)計(jì)信息。從表中可以看出，3點(diǎn)牛頓中心差分器輸出結(jié)果的精度最優(yōu)，但隨著差分器點(diǎn)數(shù)的增加，輸出加速度的精度也隨之降低，同時(shí)也可認(rèn)為采用牛頓中心差分器結(jié)合低通濾波器的處理方法是比較可靠的，其絕對(duì)精度優(yōu)于1.8mGal。

表3 基準(zhǔn)站B 的加速度絕對(duì)誤差統(tǒng)計(jì)信息Tab.3 Statistics of absolute acceleration error of station B

3.2 動(dòng)態(tài)實(shí)驗(yàn)

數(shù)據(jù)來(lái)源于美國(guó)某地區(qū)的某次航空重力測(cè)量結(jié)果。測(cè)區(qū)位于路易斯安那州，航線(xiàn)高度約為10 668m，飛機(jī)的地面速度約為155m／s，數(shù)據(jù)采樣率為1 Hz，觀測(cè)時(shí)間為6.6h。采用武漢大學(xué)研制的精密單點(diǎn)定位軟件TriP V2.0［10］解算飛機(jī)的位置信息，飛行軌跡如圖2所示。航空重力測(cè)量作業(yè)的理想情況是飛機(jī)能保持直線(xiàn)平穩(wěn)水平飛行，因此我們選取圖2中的淺色部分作為測(cè)試數(shù)據(jù)，該段數(shù)據(jù)包括3 200個(gè)歷元。

圖2 航跡圖Fig.2 Fight trajectory

對(duì)飛機(jī)z方向的位置信息作二次牛頓中心差分可得到加速度，直接差分所得加速度的頻譜包含所有頻段內(nèi)的信息。由于觀測(cè)噪聲不可避免，而差分過(guò)程會(huì)對(duì)噪聲進(jìn)行放大，如果采用理想差分器進(jìn)行數(shù)據(jù)處理，理論上在全頻段內(nèi)噪聲的頻率越高，對(duì)噪聲放大得越嚴(yán)重。但牛頓中心差分器只是對(duì)理想差分器進(jìn)行逼近，以3點(diǎn)牛頓中心差分為例，由表2 可知其截止頻率為α3＝0.603Hz。由圖3可知，在［0，α3］內(nèi)，噪聲頻率越高放大效果越明顯，甚至遠(yuǎn)遠(yuǎn)超過(guò)了信號(hào)本身的大??；而在［α3，1］范圍內(nèi)差分器對(duì)噪聲的放大程度有所減緩，說(shuō)明差分器也具有一部分低通濾波的效果。上述結(jié)果也與圖1中牛頓中心差分器的幅頻特性一致。

已有航空重力測(cè)量觀測(cè)值的頻譜分析表明，重力信號(hào)表現(xiàn)為長(zhǎng)波特性，在高頻端的能量幾乎可以忽略不計(jì)，但圖3表明，直接差分后得到的加速度的能量主要集中在高頻端，因此需對(duì)直接差分得到的加速度作低通濾波處理。低通濾波器為靜態(tài)實(shí)驗(yàn)中提及的迭代高斯低通濾波器，其對(duì)應(yīng)的截止頻率為α0＝0.015 Hz，易知α0?α3。圖3給出了對(duì)位置信息進(jìn)行兩次3點(diǎn)牛頓中心差分后所得加速度在低通濾波前后的PSD。從圖中可以看出，低通濾波后高頻噪聲已得到有效消除。

圖3 低通濾波前后加速度的PSDFig.3 PSD of acceleration before and after low－pass filtering

由于沒(méi)有真值作為比較，這里僅分別將5、7、9、11點(diǎn)牛頓中心差分器與3點(diǎn)牛頓中心差分器的輸出結(jié)果進(jìn)行比較。表4給出了低通濾波前后不同點(diǎn)數(shù)的牛頓中心差分器輸出結(jié)果的互差統(tǒng)計(jì)信息，其中（3，5）表示5點(diǎn)牛頓中心差分器與3點(diǎn)牛頓中心差分器輸出結(jié)果的較差，依此類(lèi)推。由表4可知，低通濾波前5、7、9、11點(diǎn)牛頓中心差分器與3點(diǎn)牛頓中心差分器的輸出結(jié)果較差的標(biāo)準(zhǔn)偏差在102mGal級(jí)，而濾波后較差的標(biāo)準(zhǔn)偏差已經(jīng)降低至2.1mGal以?xún)?nèi)。一方面，不同點(diǎn)數(shù)的牛頓中心差分器輸出結(jié)果之間差異較小，這是因?yàn)楦咚篂V波器的通帶［0，α0］位于所有牛頓中心差分的通帶［0，αi］（i＝3，5，7，9，11）內(nèi)，且圖1表明，不同點(diǎn)數(shù)的牛頓中心差分器的幅頻響應(yīng)在［0，α0］內(nèi)完全一致，有差異的部分則位于高斯低通濾波器的阻帶范圍內(nèi)；另一方面，因不同點(diǎn)數(shù)的牛頓中心差分器的截止頻率不一樣，在高頻段對(duì)噪聲的放大情況也不一致，同時(shí)低通濾波器不能濾除所有的高頻噪聲，這是不同點(diǎn)數(shù)的牛頓中心差分器的輸出結(jié)果之間存在較差的主要原因。此外，因3點(diǎn)牛頓中心差分器的截止頻率最小，對(duì)高頻噪聲的放大效果要弱于其他點(diǎn)數(shù)的牛頓中心差分器，同時(shí)牛頓3點(diǎn)中心差分器計(jì)算方便，使用簡(jiǎn)單，邊界效應(yīng)小。因此，在實(shí)際應(yīng)用中建議使用牛頓3點(diǎn)中心差分器。

4 結(jié) 語(yǔ)

航空重力測(cè)量的有效頻段位于低頻端，采用牛頓中心差分器得到載體的加速度后還需對(duì)加速度作低通濾波處理。本文分析了不同點(diǎn)數(shù)的牛頓中心差分器的幅頻響應(yīng)，從理論和數(shù)值實(shí)驗(yàn)兩方面驗(yàn)證了當(dāng)?shù)屯V波器的通帶位于所有牛頓中心差分器的通帶內(nèi)時(shí)，不同點(diǎn)數(shù)的牛頓中心差分器所得結(jié)果經(jīng)低通濾波后差異較小，建議航空重力測(cè)量數(shù)據(jù)處理時(shí)采用3點(diǎn)牛頓中心差分器。衛(wèi)星重力測(cè)量中高低衛(wèi)－衛(wèi)跟蹤數(shù)據(jù)差分衛(wèi)星加速度的過(guò)程也可采用類(lèi)似的方法進(jìn)行分析。

表4 低通濾波前后不同點(diǎn)數(shù)的牛頓中心差分器的輸出結(jié)果比較Tab.4 Result comparison of different point Newton central differentiator before and after low－pass filtering

［1］Ryan S，Lachapelle G，Gannon ME.DGPS Kinematic Carrier Phase Signal Simulation Analysis in the Velocity Domain［C］.ION，Kansas，1997

［2］Bruton A M，Gleeie G L，Schwarz K P.Differentiation for High－Precision GPS Velocity and Acceleration Determination［J］.GPS Solution，1999，2（4）：7－21

［3］肖云.利用GPS確定航空重力測(cè)量載體運(yùn)動(dòng)狀態(tài)的理論與方法［D］.武漢：武漢測(cè)繪科技大學(xué)，2000（Xiao Yun.Research on the Theories and Methods of Determination of Moving－Base Status for Airborne Gravimetry by Using GPS［D］.Wuhan：Wuhan Technical University of Surveying and Mapping，2000）

［4］肖云，夏哲仁.利用相位率和多普勒觀測(cè)值精確確定載體的加速度［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2003，28（5）：581－584（Xiao Yun，Xia Zheren.Comparison between Phase－Rate and Doppler to Determine Velocity［J］.Geomatics and Information Science of Wuhan University，2003，28（5）：581－584）

［5］王潛心，徐天河，龔佑興.利用中心差分法進(jìn)行GPS定速時(shí)最佳點(diǎn)數(shù)的選?。跩］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2012，37（3）：265－268（Wang Qianxin，Xu Tianhe，Gong Youxing.Optimal Points of Using Central Differences Method for GPS Velocity Determination［J］.Geomatics and Information Science of Wuhan University，2012，37（3）：265－268）

［6］歐陽(yáng)永忠.海空重力測(cè)量數(shù)據(jù)處理關(guān)鍵技術(shù)研究［D］.武漢：武漢大學(xué)，2013（Ouyang Yongzhong.Study on Key Technologies of Data Processing for Air－Sea Gravity Survey［D］.Wuhan：Wuhan University，2013）

［7］胡廣書(shū).數(shù)字信號(hào)處理理論、算法與實(shí)現(xiàn)［M］.北京：清華大學(xué)出版社，2012（Hu Guangshu.Theory，Algorithm and Implementation of Digital Signal Processing［M］.Beijing：Tsinghua University Press，2012）

［8］Cynar S J.Using Guassian Elimination for Computation of the Central Difference Coefficients［J］.Acm Signum Newsletter，1987，22（2）：12－19

［9］Hwang C，Hsiao Y S，Shih H C.Data Reduction in Scalar Airborne Gravimetry：Theory，Software and Case Study in Taiwan［J］.Computers ＆Geosciences，2006，32：1 573－1 584

［10］張小紅，劉經(jīng)南，F(xiàn)orsberg R.基于精密單點(diǎn)定位技術(shù)的航空測(cè)量應(yīng)用實(shí)踐［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2006，31（1）：19－22（Zhang Xiaohong，Liu Jingnan，Rene Forsberg.Application of Precise Point Positioning in Airborne Survey［J］.Geomatics and Information Science of Wuhan University，2006，31（1）：19－22）