馬 健，魏子卿

1. 地理信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710054； 2. 信息工程大學(xué)地理空間信息學(xué)院，河南 鄭州 450001； 3. 西安測(cè)繪研究所，陜西 西安 710054

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

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

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

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

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

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

(1)

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

(2)

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

(3)

式中，μ=Gρ，計(jì)算時(shí)ρ通常取地殼平均密度。

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

(4)

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

δζ=(Vt-Vc)/γ

(5)

當(dāng)計(jì)算點(diǎn)位于邊界面上時(shí)(即hP=0)，得到地形壓縮對(duì)大地水準(zhǔn)面高的間接影響δN為[9,19,27]

δN=(Vt-Vc)/γ0

(6)

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

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

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

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

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

2.1 直接影響

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

(7)

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

(8)

已知積分方程式

(9)

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

(10)

進(jìn)一步運(yùn)用積分式

(11)

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

(12)

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

(13)

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

(14)

2.2 間接影響

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

(15)

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

(16)

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

(17)

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

(18)

從而可得出

(19)

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

(20)

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

(21)

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

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

(22)

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

(23)

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

(24)

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

(25)

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

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

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

3 試 驗(yàn)

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

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

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

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

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兩種算法與顧及曲率的棱柱模型計(jì)算的直接影響差異

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

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

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

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

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

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

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

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′分辨率的直接影響差異較為明顯，說(shuō)明直接影響在4320～21 600頻率波段的信息含量較多。表1、表4中1°半徑的3種分辨率的間接影響差異很小，說(shuō)明間接影響在超過(guò)4320階的頻率區(qū)間的信息較少。為了說(shuō)明30″×30″、3″×3″分辨率下傳統(tǒng)積分算法計(jì)算的地形影響的誤差情況，將積分算法與顧及曲率的棱柱算法計(jì)算的地形影響差異統(tǒng)計(jì)于表5。

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

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，說(shuō)明分辨率為3″×3″時(shí)傳統(tǒng)積分算法計(jì)算的直接影響誤差基本可以忽略。但對(duì)于間接影響，試驗(yàn)區(qū)30″×30″分辨率的誤差為0.16 cm，而3″×3″分辨率的誤差為0.34 cm，此時(shí)相應(yīng)的最大誤差分別為0.68 cm、0.96 cm(絕對(duì)值)，因此當(dāng)構(gòu)建高精度的(似)大地水準(zhǔn)面時(shí)應(yīng)使用顧及曲率的棱柱模型算法計(jì)算間接影響。

4 結(jié) 論

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

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