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

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

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

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

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

一、背景

大地水準(zhǔn)面或似大地水準(zhǔn)面是獲取地理空間信息的高程基準(zhǔn)面。高精度局部大地水準(zhǔn)面將為測(cè)繪學(xué)、地球物理、地球動(dòng)力學(xué)及海洋學(xué)等相關(guān)學(xué)科的發(fā)展和應(yīng)用提供重要的基礎(chǔ)地球空間信息。我國似大地水準(zhǔn)面的精化工作雖然取得了階段性的進(jìn)展，但是與歐美國家相比還有較大的差距，同時(shí)也無法完全滿足我國目前測(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ù)估來確定重力場的方法稱為最小二乘配置法。在重力場逼近中，最小二乘配置法的最大優(yōu)點(diǎn)是可以同時(shí)利用多種不同類型的重力觀測(cè)值來確定地球外部重力場，運(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ù)處理類型不同。地球重力場實(shí)際上是一個(gè)各向異性場。由于理論上的假設(shè)，最小二乘配置法只能計(jì)算各向同性的觀測(cè)數(shù)據(jù)，在計(jì)算前需要對(duì)數(shù)據(jù)進(jìn)行中心化，而頻域輸入輸出法不存在那樣的假設(shè)，因此可以計(jì)算各向異性和各向同性的數(shù)據(jù)。而寧津生等分別對(duì)各向異性和各向同性數(shù)據(jù)恢復(fù)局部重力場，結(jié)果表明：利用各向異性的功率譜密度函數(shù)作為先驗(yàn)信息所確定的局部重力場精度要優(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ù)為我國某地區(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)面的要求。