路媛琦，章傳銀，蔣 濤

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

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

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

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

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

1 矩諧分析

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

1.1 矩諧展開(kāi)模型

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

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

(1)

Laplace方程在直角坐標(biāo)系中表示為

(2)

采用分離變量法，引力位可寫(xiě)為

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

(3)

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

(4)

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

(5)

其中

(6)

τnm定義為

(7)

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

圖1 矩諧分析使用的局部直角坐標(biāo)系

實(shí)用上矩諧級(jí)數(shù)Cnm不可能展開(kāi)到無(wú)窮階數(shù)，應(yīng)截?cái)嘀聊硞€(gè)最高階數(shù)N和次數(shù)M，則扣除參考橢球的正常引位之后，擾動(dòng)位可展開(kāi)為如下有限階次的矩諧級(jí)數(shù)

(8)

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

(9)

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

(10)

式中，r為地心距離。

同理，重力擾動(dòng)、大地水準(zhǔn)面差距的矩諧展開(kāi)式分別為

(11)

(12)

式中，γ0為正常重力。

1.2 截?cái)嚯A數(shù)

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

最高截?cái)嚯A數(shù)公式[8]為

最佳截?cái)嚯A數(shù)公式[11]為

1.3 邊界效應(yīng)

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

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

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

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

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

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

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

選取擴(kuò)展參數(shù)為截?cái)嚯A數(shù)30、擴(kuò)展參數(shù)為50 km 繪制大地水準(zhǔn)面對(duì)比圖，如圖3、圖4所示。

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

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

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

圖2 截?cái)嚯A數(shù)30時(shí)取不同擴(kuò)展參數(shù)矩諧分析所得大地水準(zhǔn)面的誤差

圖3 真實(shí)大地水準(zhǔn)面

圖4 由重力異常值所得大地水準(zhǔn)面

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

圖5 截?cái)嚯A數(shù)30時(shí)取不同擴(kuò)展參數(shù)矩諧分析所得大地水準(zhǔn)面的誤差

選取擴(kuò)展參數(shù)為截?cái)嚯A數(shù)30、擴(kuò)展參數(shù)為50 km 繪制大地水準(zhǔn)面對(duì)比圖，如圖6、圖7所示。

圖6 真實(shí)大地水準(zhǔn)面值

圖7 由重力擾動(dòng)值所得大地水準(zhǔn)面

3 分析與結(jié)論

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

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

范圍擴(kuò)展參數(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

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

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

范圍擴(kuò)展參數(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ū)域在沒(méi)有加入擴(kuò)展參數(shù)時(shí)相較于周邊區(qū)域誤差較小，但依然受到長(zhǎng)波系統(tǒng)誤差的影響；當(dāng)加入擴(kuò)展參數(shù)時(shí)，隨著擴(kuò)展參數(shù)的增大，中心3°區(qū)域大地水準(zhǔn)面的誤差逐步明顯降低，與5°區(qū)域變化趨勢(shì)一致；在擴(kuò)展參數(shù)增大到50 km時(shí)誤差值發(fā)生一個(gè)小型跳變，降低到4.2 cm；繼續(xù)增大擴(kuò)展參數(shù)(Δ=70 km)均方誤差沒(méi)有明顯變化，并不能提高大地水準(zhǔn)面的精度。

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

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

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

范圍擴(kuò)展參數(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í)取不同擴(kuò)展參數(shù)時(shí)矩諧分析所得大地水準(zhǔn)面的誤差統(tǒng)計(jì)m

范圍擴(kuò)展參數(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

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

表5取不同展開(kāi)階數(shù)時(shí)矩諧分析所得大地水準(zhǔn)面的誤差統(tǒng)計(jì)m

范圍截?cái)嚯A數(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)截?cái)嚯A數(shù)由20增加至40時(shí)誤差逐步減小，均方誤差值由5.9 cm減小至4.3 cm，誤差減小了27%；其中尤以截?cái)嚯A數(shù)從20增加至25時(shí)變化最為明顯，均方誤差減少1.1 cm；當(dāng)截?cái)嚯A數(shù)達(dá)到30時(shí)，誤差變化趨于穩(wěn)定，繼續(xù)增大截?cái)嚯A數(shù)并不能減小誤差。相較于截?cái)嗟?5階時(shí)，截?cái)嚯A數(shù)為30階時(shí)誤差分布更為均勻，最大值與最小值差距較小。而截?cái)嚯A數(shù)達(dá)到40時(shí)雖然均方誤差值不變，但計(jì)算時(shí)間大大增加，不利于進(jìn)行大量數(shù)據(jù)點(diǎn)的解算。故在實(shí)際計(jì)算中可選取截?cái)嚯A數(shù)30進(jìn)行矩諧分析。

圖8 選取擴(kuò)展參數(shù)為50 km時(shí)不同截?cái)嚯A數(shù)下矩諧分析所得大地水準(zhǔn)面的誤差

4 結(jié) 語(yǔ)

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