翟振和,孫中苗

西安測繪研究所,陜西西安710054

渤海灣多源重力數(shù)據(jù)的自適應(yīng)融合處理

翟振和,孫中苗

西安測繪研究所,陜西西安710054

在最小二乘逐步配置基礎(chǔ)上,提出多源重力數(shù)據(jù)自適應(yīng)融合處理方法,構(gòu)建基于傳統(tǒng)逐步配置的自適應(yīng)融合以及基于遞推配置的自適應(yīng)融合兩種融合模式。分別利用兩種融合模式對渤海灣陸海交界區(qū)域的航空重力數(shù)據(jù)、衛(wèi)星測高反演重力數(shù)據(jù)以及陸地重力數(shù)據(jù)進(jìn)行融合處理,其中基于遞推配置的自適應(yīng)融合模式取得了較優(yōu)的融合效果。通過兩個區(qū)域船測重力數(shù)據(jù)的外部檢核表明,融合數(shù)據(jù)的精度分別達(dá)到3.7×10-5m/s2和2.1×10-5m/s2。

多源重力數(shù)據(jù);自適應(yīng)融合;最小二乘逐步配置;融合模式

1 引 言

隨著高精尖測量技術(shù)的突飛猛進(jìn)以及相關(guān)科學(xué)理論的發(fā)展,目前在地球重力場測量領(lǐng)域已經(jīng)形成了陸、海、空、天等全方位的測量體系,重力測量數(shù)據(jù)日益豐富和多源化。為了準(zhǔn)確、可靠地解算確定地球重力場,就必須將這些具有不同頻譜、不同高度、不同誤差特性的重力數(shù)據(jù)有效融合起來,在融合過程中消除數(shù)據(jù)間的矛盾和差異從而獲得可靠、統(tǒng)一的重力基礎(chǔ)數(shù)據(jù)。在融合方法上,多數(shù)學(xué)者采用最小二乘配置(簡稱配置法)進(jìn)行融合處理[1-4]。一方面這是因為配置法具有聯(lián)合多種數(shù)據(jù)類型的天然優(yōu)勢,另一方面配置法可以提供融合結(jié)果的內(nèi)部精度估計。還有學(xué)者利用輸入-輸出系統(tǒng)理論在頻域中對重力數(shù)據(jù)進(jìn)行融合處理,這種融合實際上等效于配置法的頻域形式[5-7]。此外,也有學(xué)者如Kern、郝燕玲等利用迭代法的思想對多種重力數(shù)據(jù)進(jìn)行融合[8-9],這種融合的本質(zhì)實際上是用局部重力數(shù)據(jù)不斷地改進(jìn)原有的位模型或大地水準(zhǔn)面。在融合的具體模式上,楊元喜、曾安敏提出了大地測量數(shù)據(jù)融合的兩種模式即基于觀測信息的融合以及基于觀測信息平差結(jié)果的融合[10]。

在運(yùn)用配置法進(jìn)行融合時,多數(shù)學(xué)者都是在固定的協(xié)方差模型及模型參數(shù)下進(jìn)行解算,即所有觀測信號都采用同一個先驗信息,而后按照這個先驗信息完成融合的功能。針對觀測值協(xié)方差與觀測信號協(xié)方差比例不協(xié)調(diào)的問題,楊元喜、張菊清等提出了基于方差分量估計的擬合推估方法[11],其主要內(nèi)容就是調(diào)整對應(yīng)于觀測噪聲和觀測值協(xié)方差的單位權(quán)方差。鑒于方差分量估計的條件限制,從另一方面考慮,即便是同一協(xié)方差模型,不同的重力數(shù)據(jù)確定的模型參數(shù)也會有差異,從而導(dǎo)致重力數(shù)據(jù)的自協(xié)方差及互協(xié)方差不盡相同。這種差異很可能就包含多源重力數(shù)據(jù)的各種誤差(如系統(tǒng)偏差),當(dāng)各種數(shù)據(jù)的觀測噪聲相對固定時,若能夠?qū)Χ嘣粗亓?shù)據(jù)建立的不同協(xié)方差函數(shù)進(jìn)行深入分析,就有可能找出有效融合這些數(shù)據(jù)的處理方法。

當(dāng)配置法分步完成時就形成了逐步配置,逐步配置最大的特點是考慮到了不同觀測值的相關(guān)關(guān)系,當(dāng)不考慮觀測值的相關(guān)性及未測信號時,逐步配置則轉(zhuǎn)化為序貫平差。若將觀測值逐個加入便形成了序列配置。文獻(xiàn)[12]對兩步配置進(jìn)行了詳細(xì)的描述和驗證,在此基礎(chǔ)上文獻(xiàn)[13]對逐步配置的建立和應(yīng)用效果進(jìn)行了深入的研究分析。研究表明,逐步配置在一定觀測數(shù)量條件下可以大幅提高解算速度,更加重要的是,若將不同重力數(shù)據(jù)視為不同的觀測值,則各種重力數(shù)據(jù)對應(yīng)的協(xié)方差、互協(xié)方差將清晰地表示在逐步配置中。因此本文將利用逐步配置的優(yōu)勢,結(jié)合兩種融合模式,以渤海灣航空重力數(shù)據(jù)、衛(wèi)星測高數(shù)據(jù)、陸地重力數(shù)據(jù)的融合處理為例建立多源重力數(shù)據(jù)的自適應(yīng)融合方法。

2 最小二乘配置的遞推公式

3 基于遞推配置的自適應(yīng)融合

理論上可以找到一個滿足重力場需要且合理的協(xié)方差數(shù)學(xué)模型如 Tscherning的全球重力異常階方差模型以及 Forsberg的局部擾動位協(xié)方差模型等,但將這個數(shù)學(xué)模型運(yùn)用到局部區(qū)域時,還必須根據(jù)實測數(shù)據(jù)對模型中的參數(shù)進(jìn)行確定。在此情況下,多源重力數(shù)據(jù)中可能包含的誤差信息也必然會反映在協(xié)方差模型參數(shù)中,具體表現(xiàn)為協(xié)方差模型參數(shù)的不同。為此本文提出新的融合思路即首先在協(xié)方差模型參數(shù)中引入自適應(yīng)因子,而后利用這個自適應(yīng)因子對各種數(shù)據(jù)對應(yīng)的協(xié)方差模型參數(shù)進(jìn)行調(diào)整,直到得到最優(yōu)的結(jié)果,最終達(dá)到有效融合的目的,該融合思路本文稱之為自適應(yīng)融合。

首先討論第一種融合模式即基于遞推配置的融合。融合的具體思路就是從所有觀測值出發(fā),事先將觀測值分組,而后在遞推配置的框架內(nèi),對各種數(shù)據(jù)對應(yīng)的協(xié)方差函數(shù)加入自適應(yīng)因子,在求解過程中,一次性完成求解。為了討論方便,這里以兩種重力數(shù)據(jù)的融合為例對基于觀測信號的自適應(yīng)融合進(jìn)行分析。

假設(shè)需要融合的兩類重力數(shù)據(jù)分別為Δg1、Δg2,模型參數(shù)假設(shè)為ai,i=1,2,…,n。首先利用兩類數(shù)據(jù)分別確定各自對應(yīng)的參數(shù)得到 a1i、a2i。運(yùn)用遞推配置的通用公式得到系統(tǒng)參數(shù)的估值,為

未測信號的估值為

式(9)和式(10)中具體符號含義見前面敘述,從上式看出,兩步遞推配置中需要計算的協(xié)方差函數(shù)有CS1、CS2、C21、C12、C22、C11。

經(jīng)過一定的試驗分析,確定以下計算方案。兩種數(shù)據(jù)的自協(xié)方差分別采用各自數(shù)據(jù)確定的參數(shù)計算,即Δg1的自協(xié)方差C11利用參數(shù)a1i計算,Δg2的自協(xié)方差C22利用參數(shù) a2i計算。對于兩種數(shù)據(jù)互協(xié)方差的計算,首先加入兩個自適應(yīng)因子 P1、P2對模型參數(shù)進(jìn)行調(diào)整,為了保證加入自適應(yīng)因子后的參數(shù)符合實際情況,P1、P2應(yīng)滿足如下關(guān)系

利用 P1、P2對模型參數(shù)進(jìn)行調(diào)整得到模型新參數(shù)

獲得Ai后,計算兩種數(shù)據(jù)的互協(xié)方差 C21、 C12。未測信號與已測信號的互協(xié)方差 CS1、CS2也按照調(diào)整后的模型參數(shù)Ai計算。

在實際應(yīng)用中,模型新參數(shù)Ai需要根據(jù) P1、P2的循環(huán)調(diào)整來最優(yōu)確定(如 P1取值為0.1,則P2取值為0.9,若取循環(huán)間隔為0.05,則有20組自適應(yīng)因子)。考慮到融合結(jié)果的內(nèi)部精度不一定完全反映融合結(jié)果的真實精度,同時在融合區(qū)域一般有少量基準(zhǔn)數(shù)據(jù)可以利用,因此本文將以一定外部基準(zhǔn)數(shù)據(jù)為循環(huán)截止的評定標(biāo)準(zhǔn)。兩種重力數(shù)據(jù)融合處理的流程圖見圖1。

圖1 兩種重力數(shù)據(jù)融合處理流程圖Fig.1 The flow chart of two kind of gravity data fusion

對于多種重力數(shù)據(jù)的融合,可以按照兩種重力數(shù)據(jù)的融合進(jìn)行參考分析。此時不僅要在協(xié)方差模型參數(shù)中引入自適應(yīng)因子,還要對各種協(xié)方差函數(shù)的具體計算進(jìn)行詳細(xì)討論。以三種重力數(shù)據(jù)為例進(jìn)行分析,將三種重力數(shù)據(jù)分成三組觀測值按照三步遞推配置解算,需要計算的協(xié)方差函數(shù)有 CS1、CS2、CS3、C11、C22、C33、C21、C12、C31、C32、C13、C23。對于三種數(shù)據(jù)確定的三組協(xié)方差模型參數(shù) a1i、a2i、a3i,首先引入三個自適應(yīng)因子 P1、P2、P3對模型參數(shù)進(jìn)行調(diào)整,為了保證加入自適應(yīng)因子后的參數(shù)符合實際情況,P1、P2、P3同樣應(yīng)滿足如下關(guān)系

而后利用這些自適應(yīng)因子對三種數(shù)據(jù)對應(yīng)的參數(shù)進(jìn)行調(diào)整得到模型新參數(shù)

最后利用加入自適應(yīng)因子的參數(shù)Ai按照兩種重力數(shù)據(jù)融合的思路對各種數(shù)據(jù)的互協(xié)方差函數(shù)進(jìn)行調(diào)整從而完成三種數(shù)據(jù)的融合處理。

4 基于逐步配置的自適應(yīng)融合

按照基于逐步配置的自適應(yīng)融合思路處理多種重力數(shù)據(jù)的融合問題時,首先將所有數(shù)據(jù)按照類型分為若干組的觀測值,而后將其中兩組觀測值按照上一節(jié)的融合方案進(jìn)行融合,再對融合結(jié)果與新的觀測值進(jìn)行融合,以此類推,完成整個數(shù)據(jù)的融合。在每次融合時,都在兩步配置下通過自適應(yīng)因子進(jìn)行調(diào)整以達(dá)到最優(yōu)。在具體計算過程中,由于融合順序以及融合方式的不同使得基于逐步配置的融合具有多種形式。以三種重力數(shù)據(jù)為例進(jìn)行分析,假設(shè)有三種重力數(shù)據(jù)Δg1、Δg2、Δg3,其具體融合形式有多種,選取典型的三種見表1。

表1 基于逐步配置的自適應(yīng)融合的三種形式Tab.1 Three fusion modes based on the stepwise collocation

5 算例分析

選取渤海灣某區(qū)域的航空重力數(shù)據(jù)(如圖2直線所示區(qū)域)、衛(wèi)星測高反演的重力數(shù)據(jù)以及陸地重力數(shù)據(jù)進(jìn)行融合處理,所有數(shù)據(jù)都已歸算為5′分辨率的格網(wǎng)數(shù)據(jù),地形數(shù)據(jù)采用由 Smith. W.H.F、D.T.Sandwell構(gòu)建的基于海平面的全球2′地形數(shù)據(jù)[14]。融合前,分別選取近海(圖2中A表示)、遠(yuǎn)海(圖2中B表示)和陸地(圖2中C表示)三個代表性區(qū)域進(jìn)行精度統(tǒng)計,統(tǒng)計時,區(qū)域A、B以該區(qū)域的船測重力數(shù)據(jù)為基準(zhǔn)進(jìn)行比較,區(qū)域 C以陸地數(shù)據(jù)為基準(zhǔn)進(jìn)行比較,航空重力數(shù)據(jù)先延拓至比較區(qū)域再進(jìn)行比較。具體比較結(jié)果見表2。

表2 三種重力數(shù)據(jù)融合前的精度統(tǒng)計Tab.2 The accuracy statistic of three kinds of gravity data before fusion /(10-5m/s2)

通過三個區(qū)域的對比分析可以發(fā)現(xiàn)各種重力數(shù)據(jù)之間差異較大,具體表現(xiàn)在精度差別較大,且存在一定的系統(tǒng)偏差,如航空重力數(shù)據(jù)未經(jīng)嚴(yán)密處理可能有4×10-5m/s2的偏差,而測高數(shù)據(jù)可能有2× 10-5m/s2的偏差。值得注意的是衛(wèi)星測高數(shù)據(jù)在靠近大陸部分精度與遠(yuǎn)海區(qū)域相比明顯下降,這也說明了由于邊界效應(yīng)、近岸海域復(fù)雜的環(huán)境等因素使得近海的測高數(shù)據(jù)的質(zhì)量和可靠性并不理想。

具體融合時,考慮到重力數(shù)據(jù)分布于不同高度,因此選取近似三維的重力異常協(xié)方差模型[15],其協(xié)方差模型如下

該模型在Forsberg局部擾動位模型基礎(chǔ)上推導(dǎo)而來,模型構(gòu)建過程中顧及了重力數(shù)據(jù)隨高度產(chǎn)生的變化,能夠滿足近地空間不同高度重力數(shù)據(jù)的融合要求。模型中共有3個參數(shù),即高、低頻衰減因子以及比例系數(shù),三種重力數(shù)據(jù)對應(yīng)協(xié)方差模型參數(shù)值見表3。

表3 三種重力數(shù)據(jù)對應(yīng)的協(xié)方差模型參數(shù)值Tab.3 The covariance model parameter of three kinds of gravity data

融合的范圍包括A、B、C三個區(qū)域在內(nèi)的陸海交界區(qū)域,最終形成2°×3°范圍的融合區(qū)域,如圖2中大方框所示。融合過程中,選取的基準(zhǔn)數(shù)據(jù)為該區(qū)域1°×30′的船測重力數(shù)據(jù)(共計72個5′格網(wǎng)數(shù)據(jù),如圖2中區(qū)域D所示)。在兩種融合模式下,得到融合區(qū)域的5′平均重力異常,等值線圖見圖3,其中基于逐步配置的融合選取了兩種形式。最終融合效果通過區(qū)域 A和區(qū)域B的船測重力數(shù)據(jù)進(jìn)行評定,統(tǒng)計結(jié)果見表4。

圖3 渤海灣區(qū)域融合數(shù)據(jù)等值線圖Fig.3 The isoline figure of fusion data in Bohai gulf

表4 區(qū)域A和區(qū)域B在兩種融合模式下的自適應(yīng)融合結(jié)果Tab.4 The adaptive fusion results under two modes in area Band areaA /(10-5m/s2)

通過表4的融合結(jié)果可以看出,基于遞推配置的自適應(yīng)融合結(jié)果總體較優(yōu),其中A區(qū)精度優(yōu)于3.7×10-5m/s2,B區(qū)優(yōu)于2.1×10-5m/s2,且融合后無明顯的系統(tǒng)偏差存在。在計算過程中,衛(wèi)星測高數(shù)據(jù)、航空重力數(shù)據(jù)、陸地重力數(shù)據(jù)在區(qū)域A和區(qū)域B的三個最優(yōu)自適應(yīng)因子完全相同,分別為 P1=0.5、P2=0.4、P3=0.1。

理論上,自適應(yīng)融合的兩種模式在融合結(jié)果上應(yīng)該是完全一致的,但從融合效果上看,兩種模式獲得的融合效果不盡相同,這也再次印證了兩種融合模式在實際應(yīng)用中確實存在一定的差異。從兩種融合的過程可以看出,基于遞推配置的融合將所有觀測信號同時調(diào)整,在解算時以遞推的形式完成計算,這種融合不僅顧及了未測信號與觀測信號之間的相互關(guān)系,同時也充分利用了原始觀測信號中的有用信息,通過自適應(yīng)因子的調(diào)整就有可能提高各種數(shù)據(jù)有用信息的權(quán)重而降低不利信息的權(quán)重。而基于逐步配置的融合事先已經(jīng)將兩種觀測信號進(jìn)行了融合,而后再將其融合結(jié)果與未測信號進(jìn)行融合,這種融合實際上已經(jīng)部分改變了觀測信號的原始信息。因此總體上基于遞推配置的融合模式要優(yōu)于基于逐步配置的融合模式,這與文獻(xiàn)[10]中得到的觀測信息的融合要優(yōu)于觀測信息平差結(jié)果的融合結(jié)論基本一致。

6 結(jié) 語

通過對最小二乘逐步配置的深入分析,建立了兩種自適應(yīng)融合模式即基于遞推配置的自適應(yīng)融合以及基于逐步配置的自適應(yīng)融合,利用這兩種融合模式對渤海灣區(qū)域的航空重力數(shù)據(jù)、陸地重力數(shù)據(jù)以及衛(wèi)星測高數(shù)據(jù)進(jìn)行了融合處理,獲得了該區(qū)域2°×3°范圍的5′格網(wǎng)數(shù)據(jù),在兩個區(qū)域的檢核精度分別優(yōu)于3.7×10-5m/s2和2.1× 10-5m/s2。通過比較分析表明,自適應(yīng)融合方法能夠在一定程度上解決多源重力數(shù)據(jù)的融合問題,其中基于傳統(tǒng)逐步配置的自適應(yīng)融合的優(yōu)勢在于能夠減少大量數(shù)據(jù)融合中求逆的運(yùn)算量,而基于遞推配置的自適應(yīng)融合更能充分利用各種觀測數(shù)據(jù)的有用信息,從而獲得較優(yōu)的融合效果。

[1] OLESEN A V,ANDERSEN O B,TSCHERNING C C. Merging of Airborne Gravity and Gravity Derived from Satellite Altimetry:Test Cases along the Coast of Green-land[J].Stud Geophys Geod,2002,(46):387-394.

[2] STRYKOWSKI G,FORSBERG R.Operational Merging of Satellite,Airborne and Surface Gravity Data by Draping Techniques[C]∥Geodesy on the Move:Gravity,Geoid, Geodynamics and Antarctica.Berlin:SpringerVerlag, 1997:243-248.

[3] BAYOUD F A.Some Investigations on Local Geoid Determination from Airborne Gravity Data[D]. Calgary: University of Calgary,2001:33-35.

[4] HWANG Cheinway,GUO Jinyun,DENG Xiaoli,et al. Coastal Gravity Anomalies from Retracked Geosat/GM Altimetry:Improvement,Limitation and the Roleof Airborne Gravity Data[J].Journal of Geodesy,2006,80 (4):204-216.

[5] SCHWARZ K P,SIDERIS M G,FORSBERG R.The Use of FFT Techniques in Physical Geodesy[J].Geophys J Int, 1990,100:485-514.

[6] SIDERIS M G.On the Use of Heterogeneous Noisy Data in Spectral Gravity Field Modeling Methods[J].Journal of Geodesy,1996,70(8):470-479.

[7] TZIAVOS I N,SIDERIS M G,VERGOS G S,et al.An Overview of Spectral Methods for the Optimal Processing of Satellite Altimetry and Other Data[C]∥Proceedings of the Symposium on 15 Years of Progress in Radar Altimetry. Venice:European Space Agency,2006:135-138.

[8] HAO Yanling,CHENG Yi,LIU Fanming,et al.Simulation of Combination Algorithm for Heterogeneous Marine Gravity Data[J].Journal of System Simulation,2007,19 (21):4897-4899.(郝燕玲,成怡,劉繁明,等.融合多類型海洋重力數(shù)據(jù)算法仿真研究[J].系統(tǒng)仿真學(xué)報,2007,19 (21):4897-4899.)

[9] KERN M,SCHWARZ K K P P,SNEEUW N.A Study on the Combination of Satellite,Airborne,and Terrestrial Gravity Data[J].Journal of Geodesy,2003,77(3-4):217-225.

[10] YANG Yuanxi,ZENG Anmin.Fusion Modes of Various Geodetic Observations and Their Analysis[J].Geomatics and Information Science of Wuhan University,2008,33 (8):771-774.(楊元喜,曾安敏.大地測量數(shù)據(jù)融合模式及其分析[J].武漢大學(xué)學(xué)報:信息科學(xué)版,2008,33(8): 771-774.)

[11] YANG Yuanxi,ZHANGJuqing,ZHANG Liang.Variance Component Estimation Based Collocation and Its Application in GIS Error Fitting[J].Acta Geodaetica et Cartographica Sinica,2008,37(2):152-155.(楊元喜,張菊清,張亮.基于方差分量估計的擬合推估及其在GIS誤差糾正的應(yīng)用[J].測繪學(xué)報,2008,37(2):152-155.)

[12] LU Zhonglian.Theory and Method of Earth Gravity Field [M].Beijing:Liberation Army Publishing House,1996: 316-320.(陸仲連.地球重力場理論與方法[M].北京:解放軍出版社,1996:316-320.)

[13] ZHAI Zhenhe,SUN Zhongmiao,XIAO Yun.The Realization and Application of Least Stepwise Estimation[J]. Geomatics and Information Science of Wuhan University, 2009,34(7):761-763.(翟振和,孫中苗,肖云.最小二乘逐步、序列推估的實現(xiàn)和應(yīng)用分析[J].武漢大學(xué)學(xué)報:信息科學(xué)版,2009,34(7):761-763.)

[14] SMITH W H F,SANDWELL D T.Global Sea Floor Topography from Satellite Altimetry and Ship Depth Soundings[J].Science,1997,277(5334):1956-1962.

[15] ZHAI Zhenhe,SUN Zhongmiao.Continuation Model Construction and Application of LocalGravity Field Based on LeastSquare Collocation[J]. Chinese J Geophys,2009,52(7):1700-1706.(翟振和,孫中苗.基于配置法的局部重力場延拓模型構(gòu)建與應(yīng)用分析[J].地球物理學(xué)報,2009,52(7):1700-1706.)

(責(zé)任編輯:叢樹平)

The Adaptive Fusion of Multi-source Gravity Data in Bohai Gulf

ZHAI Zhenhe,SUN Zhongmiao
Xi’an Research Institute of Surveying and Mapping,Xi’an 710054,China

Based on the least square stepwise collocation,the adaptive fusion method is presented which has two fusion modes including the mode based on the traditional stepwise collocation and the mode based on the new recursion collocation.The airborne gravity data,satellite altimetry gravity data and land gravity data are fused by using the two fusion modes.The computation results show that the fusion based on recursion collocation has better accuracy and reliability than the other modes.According to the extra-verification with ship borne gravity data in two areas,the fusion accuracy is about 3.7×10-5m/s2and 2.1×10-5m/s2respectively.

multi-source gravity data;adaptive fusion;least square stepwise collocation;fusion mode

ZHAI Zhenhe(1980—),male,master, engineer,majors in gravity datafusion,geodetic boundary value problem,etc.

E-mail:zhaizhenhe1980@163.com

1001-1595(2010)05-0444-06

P228

A

2009-07-01

2010-03-09

翟振和(1980—),男,碩士,工程師,研究方向為重力數(shù)據(jù)融合處理、大地測量邊值問題等。