馬 健，魏子卿

1. 地理信息工程國家重點實驗室，陜西 西安 710054； 2. 信息工程大學(xué)地理空間信息學(xué)院，河南 鄭州 450001； 3. 西安測繪研究所，陜西 西安 710054

應(yīng)用Stokes/Hotine邊值理論解算(似)大地水準面時要求邊界面外不存在地形質(zhì)量[1-2]，對此可通過Helmert第二壓縮法將邊界面外的地形壓縮成覆蓋在邊界面上的質(zhì)量薄層，由此產(chǎn)生了Stokes-Helmert大地邊值理論[3]并逐步應(yīng)用到實踐中[4-7]。按照Helmert第二壓縮法進行地形壓縮時，地形質(zhì)量的移動對重力產(chǎn)生直接影響，對大地水準面起伏/高程異常產(chǎn)生間接影響。在基于Helmert第二壓縮法的邊值問題中，直接影響將重力數(shù)據(jù)從真實空間轉(zhuǎn)化到Helmert空間，而間接影響將Stokes/Hotine積分得到的Helmert高程異常/大地水準面高從Helmert空間恢復(fù)到真實空間。由于直接、間接影響直接參與(似)大地水準面解算，因此提高直接、間接影響的計算精度對于生成高精度的(似)大地水準面有著重要意義。

目前國外對地形直接、間接影響算法的研究相對較多，其中，文獻[8—9]給出了地形直接、間接影響的球諧表達式，文獻[10—11]提出綜合嚴格積分方法和球諧方法計算地形影響的算法。國內(nèi)在解算大地水準面中地形的研究主要為局部地形改正[12-16]、層間改正算法[17]等方面，對Helmert壓縮理論中地形直接、間接影響研究相對較少。隨著國內(nèi)專家學(xué)者對基于Helmert第二壓縮法的邊值解算的應(yīng)用與研究[4，18]的深化，對地形直接、間接影響展開研究是非常有必要的。

為了兼顧計算精度和計算效率，近區(qū)地形影響通常采用積分算法計算，遠區(qū)則使用Molodensky截斷理論或FFT技術(shù)結(jié)合高程球諧模型計算[19-21]。本文主要研究近區(qū)地形影響的計算。目前使用的地形直接、間接影響算法仍是二重積分形式，該算法在中央?yún)^(qū)(與計算點位于同一地心向徑的積分點所在的地形格網(wǎng)，也可稱為中央網(wǎng)格)存在奇異性。實踐中通常采用圓柱近似或借助球面布格層等方法單獨計算中央?yún)^(qū)地形影響[22]，增大了計算的復(fù)雜性。

目前地形質(zhì)量產(chǎn)生的引力位和引力的棱柱模型算法已經(jīng)比較成熟[23-24]，而局部地形改正的棱柱算法在實踐中已經(jīng)獲得廣泛應(yīng)用[14]，但尚未出現(xiàn)關(guān)于地形直接、間接影響的棱柱模型算法的研究。根據(jù)定義，地形改正是計算點周圍地形起伏對計算點產(chǎn)生的引力作用，而地形直接、間接影響是邊界面外地形壓縮到邊界面上產(chǎn)生的引力和引力位，其本質(zhì)均是地形質(zhì)量(及其壓縮形式的變體)產(chǎn)生的引力位的函數(shù)[25]，因此理論上應(yīng)可得到直接、間接影響的棱柱模型算法。鑒于此，本文通過一系列推導(dǎo)給出了近區(qū)地形直接、間接影響的棱柱模型公式。為了避免棱柱模型算法本身存在平面近似誤差，可進一步根據(jù)文獻[24]提出的方法方便地轉(zhuǎn)化為顧及地球曲率的棱柱模型算法。地形影響的棱柱模型算法不僅提高了地形影響的計算精度，而且在中央?yún)^(qū)不存在奇異性，無需單獨計算中央?yún)^(qū)地形影響，從而使地形影響傳統(tǒng)的“中央?yún)^(qū)影響+近區(qū)影響+遠區(qū)影響”計算模式轉(zhuǎn)變?yōu)椤敖鼌^(qū)影響+遠區(qū)影響”計算模式。

1 地形影響的傳統(tǒng)積分公式

使用Helmert第二壓縮法進行地形壓縮時，需依據(jù)一定的壓縮準則。本文采用質(zhì)量守恒的壓縮準則，即在球近似下壓縮面密度σ與積分網(wǎng)格的地形高度h的關(guān)系為[26]

(1)

式中，ρ為地形密度；R為地球平均半徑。將大地水準面外真實地形產(chǎn)生的引力位表示為Vt，壓縮地形產(chǎn)生的引力位表示為Vc，其算法分別為

(2)

直接影響的積分算法公式為[8,19,27]

(3)

式中，μ=Gρ，計算時ρ通常取地殼平均密度。

地形壓縮產(chǎn)生的殘余位(真實地形產(chǎn)生的引力位與壓縮地形產(chǎn)生的引力位之差)的積分算法公式為[9,19,22,27]

(4)

根據(jù)Bruns公式，當計算點位于地面上時，得到地形壓縮對高程異常的間接影響δζ為[9,19,27]

δζ=(Vt-Vc)/γ

(5)

當計算點位于邊界面上時(即hP=0)，得到地形壓縮對大地水準面高的間接影響δN為[9,19,27]

δN=(Vt-Vc)/γ0

(6)

式中，γ為地面點正常高處的正常重力；γ0為參考橢球面上的正常重力。

2 近區(qū)地形影響的棱柱模型算法推導(dǎo)

鑒于目前地形影響積分公式還是二重積分的形式，以下對地形影響的棱柱模型算法展開推導(dǎo)。通常地形是以網(wǎng)格形式存儲的，在空間直角坐標系下，任一地形網(wǎng)格Π與計算點P的相對位置如圖1所示。

圖中P為計算點，坐標系原點O為計算點在邊界面(大地水準面或參考橢球面)上的投影，因此有xP=0、yP=0。x軸指向北、y軸指向東，x1、x2分別為棱柱南、北邊界的x坐標值，y1、y2分別為棱柱西、東邊界的y坐標值。hP為計算點高程，h為地形網(wǎng)格(實際計算中代表流動積分點)的高程。

圖1 地形網(wǎng)格與計算點相對位置示意圖Fig.1 The schematic diagram of a topographic grid away from the computing point

2.1 直接影響

根據(jù)直接影響的定義，直角坐標系下直接影響的公式為

(7)

式(7)可化為下列二重積分

(8)

已知積分方程式

(9)

式中，C為不定積分式中的任意常數(shù)。根據(jù)式(9)可得式(8)的一重積分形式為

(10)

進一步運用積分式

(11)

可將式(10)轉(zhuǎn)化為

(12)

若f(x)為變量x的函數(shù)，g(y)為變量y的函數(shù)，根據(jù)式(13)

(13)

可對式(12)進一步簡化得到直接影響的棱柱模型算法公式

(14)

2.2 間接影響

由于間接影響可由殘余位通過Bruns公式得到，因此這里首先對殘余位的棱柱模型算法進行推導(dǎo)。直角坐標系下殘余地形位的基本公式為

(15)

運用式(9)可將式(15)轉(zhuǎn)化為下列二重積分

(16)

根據(jù)式(9)、式(11)，將式(16)轉(zhuǎn)化為一重積分

(17)

由于當ab<1時，反正切函數(shù)具有下列性質(zhì)

(18)

從而可得出

(19)

另外，根據(jù)分部積分法可以得到下列不定積分成立

(20)

根據(jù)式(19)、式(20)及式(13)，由式(17)得到殘余地形位的棱柱模型公式為

(21)

將式(21)分別代入式(5)、式(6)即為地形壓縮對高程異常、大地水準面高的間接影響。

地球是一個曲面，因此地形直接、間接影響的棱柱模型算法將近區(qū)范圍視為平面，存在平面近似誤差，當近區(qū)范圍較大時，需顧及地球曲率的影響。根據(jù)文獻[24]提出的算法，當顧及地球曲率的影響時需對計算點高程進行如下改化

(22)

下面從理論上對傳統(tǒng)積分算法存在的誤差進行說明。由于角距ψ是緯度θ、經(jīng)度λ的函數(shù)，因此可將式(3)、式(4)的直接、間接影響積分算法簡化為下列形式

(23)

當求取二重積分式(23)的單個網(wǎng)格的離散積分結(jié)果時，計算中采用的近似公式為

(24)

式中，Δθ=θ2-θ1；Δλ=λ2-λ1。這一算法實際上是以網(wǎng)格中心點處積分核的大小作為網(wǎng)格積分核的平均值。對于一元積分

(25)

令f(θ)=sinθ，上下界θ2=π/2、θ1=0，則式(25)左端等于1/2，右端等于π/4，顯然是不成立的。當θ1、θ2差值減小時，式(25)兩端的差異(誤差)也將隨之縮小。需要特別說明的是，當f為kθ/cosθ(k為常數(shù)，此時積分核與自變量為線性關(guān)系)時，式(25)是恒成立的。上述對一元積分式(25)的分析同樣適用于二重積分式(24)。

傳統(tǒng)積分算法是二重積分形式，由于地形影響積分核函數(shù)與幾何位置并非簡單的線性關(guān)系。通過以上分析可知，傳統(tǒng)積分算法計算地形影響時以網(wǎng)格中心點處積分核的大小作為積分單元積分核的平均值，在一定程度上引入了近似誤差。由于棱柱模型公式不存在此項誤差，因而可提高地形影響的計算精度。此外，棱柱模型公式在中央?yún)^(qū)不存在奇異性，實現(xiàn)了中央網(wǎng)格與其他網(wǎng)格算法的統(tǒng)一。

需要說明的是，同地形改正的棱柱模型算法一樣，當x1、x2、y1、y2中任意一個值接近0時，棱柱算法存在奇異性，但當滿足這一條件時傳統(tǒng)積分算法是非奇異的，因此對于滿足這一條件的個別網(wǎng)格，計算地形影響時使用傳統(tǒng)積分算法即可。

3 試 驗

為了定量說明本文推導(dǎo)的棱柱模型算法的精度，本文在我國中部地區(qū)進行了地形影響試驗，區(qū)域范圍為104°E～108°E，32°N～36°N，海拔最高3941 m，最低415 m，平均1562 m。由于顧及地球曲率的棱柱模型算法既避免了未顧及地球曲率的棱柱算法的平面近似誤差，又克服了傳統(tǒng)積分算法以格網(wǎng)中心點處積分核大小作為格網(wǎng)平均值存在的近似誤差，因此該算法計算的地形影響更加準確、可靠。本文試驗中將顧及地球曲率的棱柱模型計算的地形影響視為地形影響的準確值。

利用顧及地球曲率的棱柱算法分別計算了2.5′×2.5′分辨率下0.5°×0.5°、1°×1°、3°×3°共3種積分半徑的直接、間接影響，其中間接影響為地形壓縮對高程異常的間接影響(下同)。

表1中3種積分半徑下直接、間接影響各自的差異較小，說明近區(qū)地形影響主要受計算點附近地形的影響，而積分半徑對其影響不明顯(1 mGal=1×10-5m/s2)。為了說明不顧及地球曲率的棱柱算法的平面近似誤差和傳統(tǒng)積分算法以網(wǎng)格中點積分核大小作為網(wǎng)格積分核平均值的近似誤差的量級，將兩種算法計算的地形影響分別與表1中顧及地球曲率的棱柱模型計算的地形影響作差。表2對直接影響的差異進行了統(tǒng)計。

表1顧及曲率的棱柱模型計算的不同積分半徑的地形影響

Tab.1Topographiceffectscalculatedbyprismmodelsconsideringcurvatureunderdifferentintegralradius

地形影響算法最小值最大值平均值RMS直接影響/(mGal)0.5°×0.5°-64.68061.066-2.06910.7671°×1°-63.96264.690-0.67210.5483°×3°-63.45266.8960.43910.545間接影響/m0.5°×0.5°0.0080.7510.1520.1861°×1°0.0070.7560.1510.1863°×3°0.0010.7410.1410.177

表2兩種算法與顧及曲率的棱柱模型計算的直接影響差異

Tab.2 The direct effect differences between the two algorithms and the prism model considering curvature mGal

由于顧及曲率的棱柱算法計算的地形影響近似作為地形影響的真實值，因此表2中未顧及曲率與顧及曲率的棱柱算法差異反映了未顧及曲率的棱柱算法的平面近似誤差，而積分算法與顧及曲率的棱柱算法差異反映了積分算法以網(wǎng)格中點的積分核大小作為網(wǎng)格平均值的計算模式引入的近似誤差。從表2可以看出，兩種誤差隨積分半徑變化不明顯。3種積分半徑下，未顧及地球曲率的直接影響棱柱算法優(yōu)于直接影響的傳統(tǒng)積分算法：未顧及曲率的棱柱模型計算的直接影響的平面近似誤差很小，僅為0.06 mGal，而傳統(tǒng)積分算法計算的直接影響誤差達3.36 mGal。表3統(tǒng)計了兩種算法計算的間接影響與顧及曲率的棱柱模型計算的間接影響之間的差異。

表3兩種算法與顧及曲率的棱柱模型計算的間接影響差異

Tab.3 The indirect effect differences between the two algorithms and the prism model considering curvature cm

由于直接、間接影響具有不同的頻譜特性，兩種算法的直接、間接影響誤差也呈現(xiàn)不同的特點。表3說明了未顧及曲率的棱柱算法計算的間接影響誤差隨積分半徑的增大而增大：0.5°積分半徑下誤差僅為0.08 cm，比積分算法的間接影響誤差約優(yōu)一個數(shù)量級；1°積分半徑下誤差達到0.17 cm；當積分半徑為3°時，未顧及曲率的棱柱模型算法計算的間接影響誤差明顯增大，已超過積分算法計算的間接影響誤差，說明當間接影響的積分半徑較大時，應(yīng)使用顧及地球曲率的棱柱算法。

為了反映地形影響在不同分辨率下的情況，表4統(tǒng)計了顧及地球曲率的棱柱算法計算的30″×30″與3″×3″分辨率的直接、間接影響，積分半徑為1°。

表4顧及曲率的棱柱算法計算的不同分辨率的地形影響

Tab.4Thetopographiceffectsbytheprismalgorithmsconsideringcurvaturewithdifferentresolutions

地形影響分辨率最小值最大值平均值RMS直接影響/(mGal)30″×30″-54.34290.3658.24014.5183″×3″-53.12892.03210.24415.694間接影響/m30″×30″0.0070.7640.1510.1873″×3″0.0070.7660.1510.187

根據(jù)Nyquist采樣定律，2.5′×2.5′網(wǎng)格含有0～4320階地形信息，30″×30″網(wǎng)格包含0～21 600頻段的地形信息，而3″×3″網(wǎng)格包含0～216 000頻段的地形信息。相同積分半徑(1°)下，表1、表4中30″×30″分辨率的直接影響與2.5′×2.5′分辨率的直接影響差異較為明顯，說明直接影響在4320～21 600頻率波段的信息含量較多。表1、表4中1°半徑的3種分辨率的間接影響差異很小，說明間接影響在超過4320階的頻率區(qū)間的信息較少。為了說明30″×30″、3″×3″分辨率下傳統(tǒng)積分算法計算的地形影響的誤差情況，將積分算法與顧及曲率的棱柱算法計算的地形影響差異統(tǒng)計于表5。

表5傳統(tǒng)積分算法與顧及曲率的棱柱算法計算的地形影響差異

Tab.5Thedifferencesofthetopographiceffectsbetweenthetraditionalintegrationalgorithmsandtheprismalgorithms

地形影響分辨率最小差值最大差值平均差值RMS直接影響/(mGal)30″×30″-0.5685.3043.8374.0403″×3″-0.158-0.015-0.1060.110間接影響/cm30″×30″-0.1870.6860.0480.1613″×3″-0.964-0.060-0.2910.341

從表5可以看出，30″×30″分辨率的直接影響誤差為4.04 mGal，3″×3″誤差僅為0.11 mGal，說明分辨率為3″×3″時傳統(tǒng)積分算法計算的直接影響誤差基本可以忽略。但對于間接影響，試驗區(qū)30″×30″分辨率的誤差為0.16 cm，而3″×3″分辨率的誤差為0.34 cm，此時相應(yīng)的最大誤差分別為0.68 cm、0.96 cm(絕對值)，因此當構(gòu)建高精度的(似)大地水準面時應(yīng)使用顧及曲率的棱柱模型算法計算間接影響。

4 結(jié) 論

使用Helmert第二壓縮法對邊界面外地形進行壓縮時，對重力數(shù)據(jù)和(似)大地水準面分別產(chǎn)生直接、間接影響。近區(qū)地形影響的傳統(tǒng)算法是積分算法，由于傳統(tǒng)積分算法為二重積分形式，計算時以網(wǎng)格中點積分核的大小作為網(wǎng)格積分核平均值存在一定的近似誤差。此外，傳統(tǒng)積分算法在中央?yún)^(qū)存在奇異性，需單獨計算中央網(wǎng)格地形影響，增加了地形影響計算的復(fù)雜度。為此，本文推導(dǎo)了近區(qū)地形直接、間接影響的棱柱模型算法，為進一步避免棱柱模型存在的平面近似誤差，可使用顧及地球曲率的棱柱模型算法計算地形影響。由于棱柱模型算法在中央?yún)^(qū)不存在奇異性，因此地形影響的計算模式由傳統(tǒng)的“中央?yún)^(qū)影響+近區(qū)影響+遠區(qū)影響”轉(zhuǎn)變?yōu)椤敖鼌^(qū)影響+遠區(qū)影響”，實現(xiàn)了算法的簡化。

當試驗區(qū)地形影響取1°積分半徑時，對于2.5′×2.5′分辨率的地形，未顧及曲率的棱柱模型計算的直接影響的平面近似誤差僅為0.06 mGal，而傳統(tǒng)積分算法的誤差達3.36 mGal，此時兩種算法計算的間接影響誤差分別為0.17 cm與0.74 cm；當使用30″×30″與3″×3″分辨率地形時，傳統(tǒng)積分算法計算的直接影響誤差分別為4.04 mGal與0.11 mGal，間接影響誤差分別為0.16 cm與0.34 cm。綜上所述，對于直接影響，未顧及曲率的棱柱模型的平面近似誤差很小，在地形分辨率較高時傳統(tǒng)積分算法的誤差也是可忽略的，但對于間接影響，在高精度要求的大地水準面計算中應(yīng)盡量使用顧及地球曲率的棱柱模型算法。