局部大地水準(zhǔn)面確定方法研究

北京中鐵第五勘察設(shè)計(jì)院集團(tuán)有限公司 崔立魯

局部大地水準(zhǔn)面確定方法研究

北京中鐵第五勘察設(shè)計(jì)院集團(tuán)有限公司 崔立魯

頻域輸入輸出法和最小二乘配置法是融合多種重力場(chǎng)信息精化局部大地水準(zhǔn)面的兩種方法。本文，筆者主要介紹了這兩種融合方法的基本原理，推導(dǎo)了相關(guān)的實(shí)用計(jì)算公式，并在理論上比較了兩者之間的異同。以融合GPS水準(zhǔn)和重力異常數(shù)據(jù)精化局部大地水準(zhǔn)面為例，通過實(shí)測(cè)數(shù)據(jù)計(jì)算分析了兩種方法在精化局部重力場(chǎng)方面的精度情況。

一、背景

大地水準(zhǔn)面或似大地水準(zhǔn)面是獲取地理空間信息的高程基準(zhǔn)面。高精度局部大地水準(zhǔn)面將為測(cè)繪學(xué)、地球物理、地球動(dòng)力學(xué)及海洋學(xué)等相關(guān)學(xué)科的發(fā)展和應(yīng)用提供重要的基礎(chǔ)地球空間信息。我國(guó)似大地水準(zhǔn)面的精化工作雖然取得了階段性的進(jìn)展，但是與歐美國(guó)家相比還有較大的差距，同時(shí)也無法完全滿足我國(guó)目前測(cè)繪生產(chǎn)的需要。

二、計(jì)算方法

1.頻域輸入輸出法。根據(jù)輸入輸出端口的不同，頻域輸入輸出法可以分為多種模型，常見的由單輸入單輸出、雙輸入單輸出、多輸入單輸出和多輸入多輸出等模型。本文，筆者以雙輸入單輸出模型為例，討論頻域輸入輸出法的原理和計(jì)算模型。雙輸入單輸出模式（Double-input Single-output System）指輸入信號(hào)為兩種不同類型的信號(hào)，而輸出信號(hào)為一種類型的信號(hào)。雙輸入單輸出系統(tǒng)（有噪聲）模型如圖1所示。

在考慮噪聲的情況下，圖1中，n1，n2分別為兩種輸入信號(hào)的噪聲，h′，Δg′分別表示帶有誤差的輸入信號(hào)，即大地水準(zhǔn)面高和重力異常。

式（1）、（2）、（3）中，F(xiàn)表示傅立葉變換，F(xiàn){y}=Y，F(xiàn){b1}+B1，F(xiàn){h′}=H′，F(xiàn){b2}=B2，F(xiàn){Δg′}=ΔG， F{e}=E。

將（3）式改寫為：

式（5）中，E*表示為E的共軛量。其中Pee表示功率譜密度（即PSD），計(jì)算公式為：

式（6）中，Cxy和Pxy分別表示信號(hào)x，y之間的協(xié)方差和功率譜密度。

在Pee=m in準(zhǔn)則下，可求得系數(shù)B1，B2。

若已知相關(guān)變量的功率譜密度，則聯(lián)合（3）、（7）、（8）式即可計(jì)算該系統(tǒng)的輸出信號(hào)。

由公式（7）、（8）式可知，如需計(jì)算輸出信號(hào)必須知道輸入信號(hào)的功率譜密度，本文計(jì)算功率譜密度的公式如下。

式（9）、（10）中，f1（i, j）和f2（k＋i,l＋j）表示2個(gè)不同格網(wǎng)觀測(cè)值，x， y分別表示格網(wǎng)數(shù)據(jù)的橫縱方向和分別為觀測(cè)值f1（i, j）和f2（k＋i,l＋j）的數(shù)學(xué)期望，M和N分別為數(shù)據(jù)點(diǎn)沿著方向上的個(gè)數(shù)，k和l分別為沿著x和y方向上的距離值。根據(jù)式（9）、（10）可計(jì)算出相應(yīng)的功率譜密度。

2.最小二乘配置法。在物理大地測(cè)量學(xué)中，把用最小二乘預(yù)估來確定重力場(chǎng)的方法稱為最小二乘配置法。在重力場(chǎng)逼近中，最小二乘配置法的最大優(yōu)點(diǎn)是可以同時(shí)利用多種不同類型的重力觀測(cè)值來確定地球外部重力場(chǎng)，運(yùn)用該方法的關(guān)鍵是要獲得一個(gè)符合實(shí)際情況的關(guān)于觀測(cè)值的協(xié)方差函數(shù)。

根據(jù)Moritz H.的理論推導(dǎo)，最小二乘預(yù)估公式為：

式（11）中，s表示推估點(diǎn)的值，l為觀測(cè)值的數(shù)據(jù)，Csl為推估值與觀測(cè)值的互相關(guān)協(xié)方差，Cll為觀測(cè)值之間的協(xié)方差。

在重力數(shù)據(jù)融合中，當(dāng)觀測(cè)值為大地水準(zhǔn)面高和重力異常，推估值為任意值時(shí)，其相應(yīng)公式為：

式（14）中，Cyh和CyΔg分別表示推估值與大地水準(zhǔn)面高和重力異常的互協(xié)方差函數(shù)，Chh，ChΔg，CΔgh和CΔgΔg分別表示大地水準(zhǔn)面高與重力異常的自協(xié)方差函數(shù)及其互協(xié)方差函數(shù)，Dhh、DΔgΔg分別表示大地水準(zhǔn)面高與重力異常的噪聲自協(xié)方差函數(shù)，其中ChΔg=ChΔgT。 N、 Δg′分別表示推估信號(hào)大地水準(zhǔn)面高和重力異常，而Cee（N）和Cee（Δg′）分別表示大地水準(zhǔn)面高和重力異常推估信號(hào)噪聲協(xié)方差矩陣。

由最小二乘配置法原理可知，獲得與觀測(cè)值相關(guān)的協(xié)方差函數(shù)是該方法計(jì)算的關(guān)鍵。本文主要采用格網(wǎng)數(shù)據(jù)，其協(xié)方差計(jì)算公式如下。

式（15）、（16）中，f*是觀測(cè)值f與該地區(qū)觀測(cè)值的平均值之差，C（S）橫和 C（S）縱分別表示橫向和縱向協(xié)方差，分別按下式計(jì)算。

由（15）和（17）式得到離散的協(xié)方差值，再由多項(xiàng)式函數(shù)擬合出經(jīng)驗(yàn)協(xié)方差函數(shù)p（x）：

式（18）中， p（x）為函數(shù)值；（p1，p2，··，pn，pn＋1）為函數(shù)的參數(shù)值，即需擬合的參數(shù)；x 為函數(shù)自變量。

3.兩種方法的比較。從兩種方法的基本原理，可以看出頻域輸入輸出法與最小二乘配置法具有一定的相似性，但也存在著一些區(qū)別。下面從四個(gè)不同的角度對(duì)兩種方法進(jìn)行比較。

（1）計(jì)算原理不同。頻域輸入輸出法是基于頻域的數(shù)據(jù)處理方法，其基本原理是根據(jù)不同數(shù)據(jù)的頻譜特性（即功率譜密度）確定其在計(jì)算結(jié)果中所占的比重；而最小二乘配置則是一種空域的數(shù)據(jù)處理方法，該方法主要是根據(jù)不同數(shù)據(jù)的數(shù)學(xué)特性（即協(xié)方差）來確定相應(yīng)的比重。兩者的區(qū)別在于數(shù)學(xué)運(yùn)算法則在頻域和空域內(nèi)是不同的，其中空域的矩陣求逆到了頻域就變成了矩陣相除，這樣就極大地減少了計(jì)算所需要的軟硬件資源，并節(jié)省了大量的時(shí)間，提高了計(jì)算效率。

（2）計(jì)算準(zhǔn)則不同。頻域輸入輸出法和最小二乘配置法的計(jì)算準(zhǔn)則分別是 和 。根據(jù)頻域輸入輸出法和最小二乘配置法原理，可知 和 在頻域和空域內(nèi)是相互對(duì)應(yīng)的，分別表示的是輸出信號(hào)噪聲的功率譜密度和推估值的誤差協(xié)方差。因此，兩種方法的計(jì)算準(zhǔn)則是一致的。

（3）數(shù)據(jù)處理類型不同。地球重力場(chǎng)實(shí)際上是一個(gè)各向異性場(chǎng)。由于理論上的假設(shè)，最小二乘配置法只能計(jì)算各向同性的觀測(cè)數(shù)據(jù)，在計(jì)算前需要對(duì)數(shù)據(jù)進(jìn)行中心化，而頻域輸入輸出法不存在那樣的假設(shè)，因此可以計(jì)算各向異性和各向同性的數(shù)據(jù)。而寧津生等分別對(duì)各向異性和各向同性數(shù)據(jù)恢復(fù)局部重力場(chǎng)，結(jié)果表明：利用各向異性的功率譜密度函數(shù)作為先驗(yàn)信息所確定的局部重力場(chǎng)精度要優(yōu)于由各向同性功率譜密度函數(shù)所得到的精度。

（4）兩種方法在計(jì)算前都需要知道觀測(cè)數(shù)據(jù)的先驗(yàn)信息，即功率譜密度和協(xié)方差。在計(jì)算過程中，將根據(jù)不同觀測(cè)數(shù)據(jù)之間的相互關(guān)系和自身的特性來確定各種類型數(shù)據(jù)在計(jì)算中所占的比重。而功率譜密度和協(xié)方差可以通過快速傅立葉變換（FFT）實(shí)現(xiàn)兩者之間的轉(zhuǎn)換。

三、數(shù)值計(jì)算與分析

實(shí)測(cè)數(shù)據(jù)為我國(guó)某地區(qū)53×101格網(wǎng)重力異常數(shù)據(jù)和65個(gè)GPS水準(zhǔn)數(shù)據(jù)，其緯度范圍為22°26′17.5″N- 22°55′00″N，經(jīng)度為113°41′00″E- 114°39′55″E，格網(wǎng)的橫、縱間隔分別為32.5″和35″。格網(wǎng)重力異常的標(biāo)準(zhǔn)差為1.600毫伽，GPS水準(zhǔn)數(shù)據(jù)的標(biāo)準(zhǔn)差為0.010 m。另有29個(gè)大致均勻分布的GPS水準(zhǔn)點(diǎn)作為外部檢核點(diǎn)。

表1 融合結(jié)果誤差

由表1可知，實(shí)測(cè)數(shù)據(jù)的計(jì)算結(jié)果，融合后的大地水準(zhǔn)面內(nèi)外符合精度均達(dá)到厘米級(jí)，滿足目前大地水準(zhǔn)面的精度要求。從實(shí)測(cè)數(shù)據(jù)計(jì)算結(jié)果可以看出：IOST的大地水準(zhǔn)面高內(nèi)符合精度為±0.020 2 m，而LSC的內(nèi)符合精度為±0.064 4 m。IOST和LSC大地水準(zhǔn)面高的外符合精度分別為±0.037 6 m和±0.042 1 m。融合結(jié)果為重力異常的計(jì)算精度也是如此。由此可知，頻域輸入輸出法的融合精度要略優(yōu)于最小二乘配置法。兩種方法的計(jì)算精度均達(dá)到精化局部大地水準(zhǔn)面的要求。