路媛琦，章傳銀，蔣 濤

(1. 青島市勘察測繪研究院，山東 青島 266032； 2. 中國測繪科學(xué)研究院，北京 100039； 3. 山東科技大學(xué)，山東 青島 266590)

大地水準面是一個與靜止的平均海水面重合并延伸到大陸內(nèi)部的、包圍整個地球的、封閉的重力位水準面，可以作為定義全球或區(qū)域高程系統(tǒng)的基準面[1]。確定大地水準面是建立高程基準統(tǒng)一的有效途徑，并可結(jié)合GNSS實現(xiàn)大地高到正高的轉(zhuǎn)換，實現(xiàn)衛(wèi)星三維空間定位。

全球大地水準面的確定通常采用球諧展開模型，但確定局部大地水準面所利用的局部重力數(shù)據(jù)無法滿足球諧函數(shù)的正交性，球諧展開模型并不適用于表達區(qū)域大地水準面。區(qū)域大地水準面確定的常用方法一種是基于Stokes或Molodensky理論求解大地測量邊值問題，我國第二代似大地水準面模型(CQG2000)即是采用此類傳統(tǒng)方法構(gòu)建。另一種被許多學(xué)者應(yīng)用的解析法為球冠諧分析；李建成等首次提出將球冠諧分析應(yīng)用于區(qū)域重力場建模，證明了球冠諧展開表達區(qū)域重力場的可行性。球冠諧分析在大區(qū)域范圍內(nèi)建立大地水準面是一種理想方法，但應(yīng)用于較小范圍的大地水準面確定時，存在數(shù)值計算上的困難，精度較差。

確定局部大地水準面的譜方法除球冠諧分析外，還有一種方法是矩諧分析。矩諧分析在地磁學(xué)研究中應(yīng)用廣泛；Alldredge最早使用矩諧分析用于區(qū)域地磁場建模。1992年邊少鋒[1]將矩諧分析引入地球重力場逼近，由于矩諧分析的基函數(shù)是三角函數(shù)，球冠諧分析的基函數(shù)是非整階締合Legendre函數(shù)，展開到較高階次時前者的計算遠比后者穩(wěn)定、快捷。

本文利用EGM-2008重力位模型，加入標準差為2 mGal的高斯白噪聲，模擬陸地重力及飛行高度為4 km的航空重力值；并在研究區(qū)域內(nèi)建立直角坐標系，求解地球引力位的Laplace方程，得到由矩諧系數(shù)表達的重力場參數(shù)表達式；將模擬重力數(shù)值代入矩諧展開式中求得矩諧系數(shù)、殘余重力擾動及殘余大地水準面；最后恢復(fù)大地水準面并與真實值進行對比，以此說明矩諧分析確定大地水準面的可靠性、準確性。

1 矩諧分析

已知地球外部空間某一區(qū)域或地球表面的重力觀測值(重力異?；蛑亓_動)，求定該區(qū)域內(nèi)地球重力場的矩諧系數(shù)，這一過程叫作矩諧分析[2]。

1.1 矩諧展開模型

矩諧分析(rectangular harmonic analysis，RHA)的基本原理是對地球某一地區(qū)構(gòu)建區(qū)域大小的矩形面，以此近似球面；以矩形面的中心點為原點建立局部直角坐標系(如圖1所示)；并在其中求解Laplace方程，在選定一定的截斷水平后，引力場和引力位是包含有限個待定系數(shù)的已知函數(shù)，可以用所研究區(qū)域內(nèi)一組觀測點上的引力場值確定這些系數(shù)，并用這些系數(shù)表示不同波長的重力場信息，便于求解大地水準面[3-6]。

地球引力位滿足Laplace方程[6-7]

(1)

Laplace方程在直角坐標系中表示為

(2)

采用分離變量法，引力位可寫為

V(x，y，z)=f(x)g(y)h(z)

(3)

將式(3)代入式(2)得局部直角坐標系中Laplace方程的解[3]為

(4)

式中，Cnm為矩諧系數(shù)，即引力位系數(shù)；ψ(x)，ψ(y)為基函數(shù)，可分別表示為

(5)

其中

(6)

τnm定義為

(7)

式中，Dx、Dy分別表示矩形計算區(qū)域在東西(x)和南北(y)方向上的距離。

圖1 矩諧分析使用的局部直角坐標系

實用上矩諧級數(shù)Cnm不可能展開到無窮階數(shù)，應(yīng)截斷至某個最高階數(shù)N和次數(shù)M，則扣除參考橢球的正常引位之后，擾動位可展開為如下有限階次的矩諧級數(shù)

(8)

由Molodensky邊值條件得重力異常與擾動位滿足如下近似關(guān)系

(9)

于是聯(lián)合式(8)和式(9)可推出重力異常的矩諧展開式為

(10)

式中，r為地心距離。

同理，重力擾動、大地水準面差距的矩諧展開式分別為

(11)

(12)

式中，γ0為正常重力。

1.2 截斷階數(shù)

在實際應(yīng)用中，無論球諧級數(shù)還是矩諧級數(shù)都不可能展開到無窮多次，因此對級數(shù)表達式選擇適當(dāng)?shù)慕財嚯A數(shù)非常重要，過低的截斷階數(shù)會丟失有價值的信息，而過高的截斷階數(shù)會導(dǎo)致計算量增加，從而使結(jié)果不穩(wěn)定，出現(xiàn)所謂的“龍格”現(xiàn)象[8-9]。因此需要根據(jù)數(shù)據(jù)的數(shù)量和精度要求[10]，確定模型的截斷階數(shù)Nmax。最優(yōu)截斷階數(shù)的確定仍舊困難，通常得到最高截斷階數(shù)后不斷進行嘗試性的計算比較，以獲取最優(yōu)截斷階數(shù)。

最高截斷階數(shù)公式[8]為

最佳截斷階數(shù)公式[11]為

1.3 邊界效應(yīng)

對于區(qū)域大地水準面而言，其研究范圍的局部性限制了重力場信號的完全周期性；但矩諧展開模型采用的周期性函數(shù)Fourier級數(shù)成立的前提是假定待求信號在計算區(qū)域內(nèi)為周期函數(shù)[12]。式(8)在計算區(qū)域內(nèi)收斂于真實地球擾動位，在計算區(qū)域邊界處近似等于兩邊界擾動位的平均值，因此在邊界處會產(chǎn)生Gibbs振蕩現(xiàn)象。因此可擴展所求區(qū)域，以此降低Gibbs震蕩現(xiàn)象的影響，使得計算區(qū)域的重力場信號具有周期性[13-15]。

舉例說明，若研究范圍大小的矩形長、寬分別為為Dx、Dy，數(shù)據(jù)點范圍長、寬同為dx、dy；可將計算區(qū)域的平面范圍擴展為Dx+dx、Dy+dy，此時待估重力場信號在區(qū)間D+d上滿足周期性條件，此時式(7)為

當(dāng)擴展參數(shù)增大到一定大小時可以有效降低邊緣效應(yīng)，但過大的擴展參數(shù)會使求解過程存在病態(tài)，同樣影響帶球擾動位系數(shù)的估算精度[3]。

2 數(shù)值計算與分析

2.1 基于地面重力數(shù)據(jù)的算例與分析

為驗證矩諧分析基于陸地重力試驗數(shù)據(jù)的可靠性，設(shè)計如下試驗：獲取緯度為[29°，34°]、經(jīng)度為[111°，116°]圍成的5°×5°區(qū)域，Dx=504.37 km，Dy=601.34 km，格網(wǎng)間隔為2.5′×2.5′，共計14 400個計算點。利用EGM2008重力位模型的2～2190階系數(shù)計算得到重力異常的模擬觀測值，并加入數(shù)學(xué)期望為0、標準差為2 mGal的高斯白噪聲；移去基于GOCO05S模型的2～200階系數(shù)計算的參考重力異常，得到殘余重力異常；根據(jù)殘余重力異常的矩諧展開式得到矩諧系數(shù)；根據(jù)大地水準面差距的矩諧表達式求得殘余大地水準面，恢復(fù)參考大地水準面得到最終的大地水準面；最后利用EGM2008模型2～2190階位系數(shù)計算的大地水準面高作為真實值，比較基于矩諧分析的結(jié)果與真實值間差距，說明矩諧分析的可靠性與精度。

根據(jù)上文所述可知，矩諧分析受周期延拓邊界效應(yīng)的制約，并成為影響其精度的關(guān)鍵因素。本試驗利用擴展參數(shù)法以減小邊界效應(yīng)的影響，圖2中虛線范圍為3°×3°中心區(qū)域。由于在不同截斷階數(shù)下大地水準面誤差變化規(guī)律相同，因此以截斷階數(shù)30為例進行邊界效應(yīng)的研究，并繪制誤差分布圖。

選取擴展參數(shù)為截斷階數(shù)30、擴展參數(shù)為50 km 繪制大地水準面對比圖，如圖3、圖4所示。

2.2 基于航空重力數(shù)據(jù)的算例與分析

由于航空重力測量獲取的是飛行高度處的重力擾動值，但所需值為大地水準面上的擾動值，而將飛行高度處的數(shù)值向下延拓至大地水準面的過程中，噪聲會被放大，重力場信息因此失真，這是一個病態(tài)求解的過程。但航空重力測量對于獲取高頻重力場信息有不可替代的作用，因此如何能在延拓過程中將重力場信息高度還原，獲得穩(wěn)定的重力場解是航空重力數(shù)據(jù)處理的難點問題。

為驗證矩諧分析的全面性、可靠性，本文基于航空重力數(shù)據(jù)設(shè)計如下試驗。獲取緯度為[29°，34°]、經(jīng)度為[111°，116°]圍成的5°×5°區(qū)域，Dx=504.37 km，Dy=601.34 km，格網(wǎng)間隔為2.5′×2.5′，共計14 400個計算點。利用EGM2008重力位模型的2～2190階系數(shù)分別計算飛行高度為4 km的重力擾動的模擬觀測值，并加入數(shù)學(xué)期望為0、標準差為2 mGal的高斯白噪聲；移去基于GOCO05S模型的2～200階系數(shù)計算的參考重力擾動，得到殘余重力擾動；根據(jù)殘余重力擾動的矩諧展開式得到矩諧系數(shù)；根據(jù)大地水準面差距的矩諧表達式和重力擾動的矩諧展開式求得殘余大地水準面；最后恢復(fù)參考大地水準面及參考重力擾動得到最終的大地水準面；最后利用EGM2008模型2～2190階位系數(shù)計算的大地水準面作為真實值，比較基于矩諧分析的計算結(jié)果與真實值間的差距。

圖2 截斷階數(shù)30時取不同擴展參數(shù)矩諧分析所得大地水準面的誤差

圖3 真實大地水準面

圖4 由重力異常值所得大地水準面

與重力異常處理方法相同，航空重力值也選用擴展參數(shù)法以減小邊界效應(yīng)的影響。由于不同截斷階數(shù)時的大地水準面誤差變化規(guī)律相同，因此選取截斷階數(shù)30為例繪制誤差分布圖，如圖5所示。

圖5 截斷階數(shù)30時取不同擴展參數(shù)矩諧分析所得大地水準面的誤差

選取擴展參數(shù)為截斷階數(shù)30、擴展參數(shù)為50 km 繪制大地水準面對比圖，如圖6、圖7所示。

圖6 真實大地水準面值

圖7 由重力擾動值所得大地水準面

3 分析與結(jié)論

如圖2、圖5所示，當(dāng)截斷階數(shù)固定時，針對每一截斷階數(shù)，分別選取0、30、50、70 km的擴展參數(shù)計算其大地水準面的誤差分布。由于各截斷階數(shù)的變化規(guī)律一致，因此選取截斷階數(shù)30為案例進行分析。表1、表2分別給出了不同研究范圍內(nèi)取不同擴展參數(shù)時矩諧分析所得的大地水準面的誤差統(tǒng)計信息。

表15°×5°范圍時取不同擴展參數(shù)時矩諧分析所得大地水準面的誤差統(tǒng)計m

范圍擴展參數(shù)Δ最小值最大值平均值均方差5°×5°0-0.3380.311-0.0150.06730-0.3590.299-0.0150.06750-0.5460.273-0.0150.06370-0.5540.269-0.0140.063

對于全圖5°×5°區(qū)域范圍，在4個角點處誤差出現(xiàn)極大值，邊界效應(yīng)明顯，對于本研究區(qū)域尤以左上角部分最為嚴重，誤差最大值達到3 dm；在擴展參數(shù)為0、30 km時均方誤差基本不變；當(dāng)擴展參數(shù)不斷增大時，整體誤差逐步降低，當(dāng)擴展參數(shù)為50 km 時，均方誤差產(chǎn)生小型跳變，最大值降低到2.7 cm，均方誤差為6.3 cm；繼續(xù)增大擴展參數(shù)，達到70 km時，并不能減小均方誤差，且會增長計算時間，影響計算效率。

表23°×3°范圍時取不同擴展參數(shù)時矩諧分析所得大地水準面的誤差統(tǒng)計m

范圍擴展參數(shù)Δ最小值最大值平均值均方差3°×3°0-0.1570.1190.0150.04830-0.1510.1160.0150.04750-0.1510.1150.0150.04270-0.1510.1150.0140.042

可以明顯看出圖示中心3°×3°區(qū)域在沒有加入擴展參數(shù)時相較于周邊區(qū)域誤差較小，但依然受到長波系統(tǒng)誤差的影響；當(dāng)加入擴展參數(shù)時，隨著擴展參數(shù)的增大，中心3°區(qū)域大地水準面的誤差逐步明顯降低，與5°區(qū)域變化趨勢一致；在擴展參數(shù)增大到50 km時誤差值發(fā)生一個小型跳變，降低到4.2 cm；繼續(xù)增大擴展參數(shù)(Δ=70 km)均方誤差沒有明顯變化，并不能提高大地水準面的精度。

按同樣的方式對航空重力數(shù)據(jù)進行統(tǒng)計，誤差統(tǒng)計見表3、表4。在不同范圍的不同擴展參數(shù)下其變化趨勢與重力異常數(shù)據(jù)一致。

綜上所述，邊界效用確實影響矩諧分析的精度，但采用擴展參數(shù)方法后，精度上升明顯。擴展參數(shù)為50 km時矩諧分析得到的大地水準面的精度最高，擴展參數(shù)過小(如0、30 km)對精度的提高無法達到最佳，當(dāng)擴展參數(shù)過大時(如70 km)，并不能繼續(xù)提高精度，相反會大大增加解算時間以及復(fù)雜程度；因此實際計算中宜選取擴展參數(shù)為50 km進行解算。

表35°×5°范圍時取不同擴展參數(shù)時矩諧分析所得大地水準面的誤差統(tǒng)計m

范圍擴展參數(shù)Δ最小值最大值平均值均方差5°×5°0-0.3340.273-0.0140.06830-0.3540.300-0.0150.06650-0.5460.275-0.0150.06370-0.5490.271-0.0140.062

表43°×3°范圍時取不同擴展參數(shù)時矩諧分析所得大地水準面的誤差統(tǒng)計m

范圍擴展參數(shù)Δ最小值最大值平均值均方差3°×3°0-0.1580.104-0.0220.04630-0.1560.107-0.0200.04450-0.1530.109-0.0180.04370-0.1540.110-0.0170.043

限制矩諧分析精度的另一大難題是截斷階數(shù)的選取。為獲得針對矩諧分析的適當(dāng)?shù)慕財嚯A數(shù)，避免選取截斷階數(shù)較大時產(chǎn)生“龍格”現(xiàn)象，截斷階數(shù)較小時無法獲取有效的數(shù)據(jù)值，本文針對擴展參數(shù)為50的不同截斷階數(shù)的矩諧展開結(jié)果進行分析，誤差分布如圖8所示。數(shù)值統(tǒng)計分析見表5。

表5取不同展開階數(shù)時矩諧分析所得大地水準面的誤差統(tǒng)計m

范圍截斷階數(shù)(N=M)最小值最大值平均值均方差3°×3°20-0.1560.119-0.0190.05925-0.1520.117-0.0160.04830-0.1530.119-0.0180.04335-0.1510.115-0.0160.04340-0.1520.116-0.0160.043

可以看出，當(dāng)截斷階數(shù)由20增加至40時誤差逐步減小，均方誤差值由5.9 cm減小至4.3 cm，誤差減小了27%；其中尤以截斷階數(shù)從20增加至25時變化最為明顯，均方誤差減少1.1 cm；當(dāng)截斷階數(shù)達到30時，誤差變化趨于穩(wěn)定，繼續(xù)增大截斷階數(shù)并不能減小誤差。相較于截斷到35階時，截斷階數(shù)為30階時誤差分布更為均勻，最大值與最小值差距較小。而截斷階數(shù)達到40時雖然均方誤差值不變，但計算時間大大增加，不利于進行大量數(shù)據(jù)點的解算。故在實際計算中可選取截斷階數(shù)30進行矩諧分析。

圖8 選取擴展參數(shù)為50 km時不同截斷階數(shù)下矩諧分析所得大地水準面的誤差

4 結(jié) 語

針對球諧分析及球冠諧分析在求解小區(qū)域大地水準面時的解算難點，本文利用更穩(wěn)當(dāng)?shù)木刂C分析求解大地水準面。為驗證理論的可靠有效性，設(shè)計完成了基于EGM2008、基于陸地重力數(shù)據(jù)和航空重力數(shù)據(jù)的試驗。試驗結(jié)果表明，利用矩諧分析求解大地水準面，基于陸地觀測數(shù)據(jù)精度可達到4 cm，基于航空重力觀測數(shù)據(jù)精度可達到4 cm。相比于應(yīng)用廣泛的stokes理論求解大地測量邊值問題所得結(jié)果，精度提高明顯。因此，矩諧分析是確定區(qū)域大地水準面的理想方法，值得應(yīng)用推廣。