歐陽明達(dá)，周 巍，張敏利

(1.地理信息工程國家重點實驗室，陜西 西安 710054;2.西安測繪總站，陜西 西安 710054)

地形改正在大地水準(zhǔn)面計算、重力勘探、地殼密度結(jié)構(gòu)研究等方面有重要應(yīng)用[1-2]?，F(xiàn)有大量文獻(xiàn)對地形改正計算提出了一系列解決方案，且通常將地形分為內(nèi)區(qū)和外區(qū)。傳統(tǒng)的平面計算方案，外區(qū)被視為無限布格平板，內(nèi)區(qū)被近似為一有限范圍且以計算點為中心的平面[3-4]。FFT方法的出現(xiàn)提高了地形改正值的計算效率，但需要提供規(guī)則化格網(wǎng)形式的地形高輸入數(shù)據(jù)[5-7]。隨著衛(wèi)星重力觀測技術(shù)的發(fā)展，以及全球性大尺度海量觀測數(shù)據(jù)的獲取，學(xué)者們開始了對球坐標(biāo)系下重力場元的研究，該方法將地表近似為貼合的球面，外區(qū)為球面布格層，其量級約為平面無限布格層的兩倍，內(nèi)區(qū)則采用球面積分方法進(jìn)行計算[8-10]。

本文重點關(guān)注局部重力地形改正，即內(nèi)區(qū)地形改正值的獲取，給出了平面積分方法、平面FFT方法、球面積分方法的計算公式，以西部某山區(qū)為例開展試算，對其結(jié)果進(jìn)行了比較，并得出了一些有益的結(jié)論。

1 平面積分

根據(jù)定義，不規(guī)則地形起伏相對于計算點P(x,y,z)產(chǎn)生的引力位為

(1)

式中，G為牛頓引力常數(shù)，ρ為地殼密度，σ代表積分區(qū)域，x、y、z、h為流動點的坐標(biāo)及高程，xp、yp、zp、hp為計算點的坐標(biāo)及高程，式(1)的負(fù)垂直方向分量即為重力的地形改正：

(2)

假設(shè)DTM以N×M的規(guī)則格網(wǎng)數(shù)據(jù)表示，每個網(wǎng)格內(nèi)的地形用一個質(zhì)量均勻分布而高程不變的棱柱體代替，則得到用質(zhì)量棱柱地形模型表示的地形改正式：

(3)

令

(4)

(5)

上式右側(cè)可改寫為[1+x]-1/2的形式，其中，x=[(hp-h)/l]2，級數(shù)展開為

(6)

從式(6)可知，地形改正計算時，流動點與計算點連線的坡度不能超過45°，很明顯，這種情況非常少見，即便在山區(qū)也與大多數(shù)地形地貌嚴(yán)重不符。將式(6)引入式(5)，則積分式可改寫為

(7)

即

δg≈C1+C2+C3.

(8)

本文取其1階近似，得到

(9)

對其進(jìn)行FFT變換，改寫成譜形式：

(10)

采用FFT方法計算時，在計算點處，核函數(shù)出現(xiàn)奇異，通常采用引入核函數(shù)項增加常數(shù)因子(α)的方法解決這一問題，即使用新的水平距離函數(shù)代替原來的函數(shù)，具體實現(xiàn)過程參考文獻(xiàn)[6]。

2 球面積分

平面積分忽視了地球曲率影響，對較遠(yuǎn)區(qū)來說存在誤差。將積分區(qū)域視為球近似，則重力的局部地形改正可表示為

(11)

式中，λ′、φ′、r′為流動點的球坐標(biāo)，L-1為牛頓積分核函數(shù)(即流動點到計算點空間距離L的倒數(shù))，且

(12)

cosψ=cosφcosφ′+sinφsinφ′cos(λ-λ′).

(13)

r(3t2-1)ln|r′-rt+L|+C.

(14)

其中，t=cosψ，C為積分常數(shù)。

計算地形改正的奇異積分問題，采用“梯度法”處理[11]。以計算點為中心，取一足夠小的球冠，將其視為平面，半徑為Δφ0=Δλ0，表達(dá)式為

(15)

其中，

(16)

[h(φp,λp+Δλ0)-h(φp,λp-Δλ0)].

(17)

3 算例與結(jié)果分析

3.1 數(shù)據(jù)準(zhǔn)備與方案設(shè)計

選擇我國西部山區(qū) 35°～38°N，95°～98°E范圍進(jìn)行計算，高程數(shù)據(jù)來自于SRTM模型，經(jīng)平滑得到分辨率為30″的高精度地形數(shù)據(jù)，如圖1所示，橫坐標(biāo)、縱坐標(biāo)為經(jīng)緯度范圍，數(shù)值單位為m。由于計算需要，地形模型需向外分別擴(kuò)充2°范圍。試驗區(qū)地形高程的最大值為5 592 m，最小值為2 675 m，平均值為3 695 m。

圖1 地形模型

3.2 結(jié)果分析與比較

根據(jù)相關(guān)文獻(xiàn)結(jié)論，選擇積分半徑為100 km[12]。圖2表示出了采用3種方案得到的局部重力地形改正值。表1給出了其相關(guān)統(tǒng)計信息，可以看出，地形改正值結(jié)果多為正值，不同方法具有顯著差異；計算點周邊地形起伏大，其改正值較大；中部盆地受較遠(yuǎn)地區(qū)復(fù)雜地形的影響較小，改正值較小，多位于0 mGal附近；局部地形改正量受山區(qū)、平原、盆地等的地貌影響較小(見圖2)，數(shù)值單位為mGal，盡管南部為山區(qū)地貌，但其上地形起伏較為平緩，因而地形改正量并不大。

圖3給出了3種方案的計算結(jié)果對比，表2給出了其相關(guān)統(tǒng)計信息。從表2和圖3可以看出，不同方案計算結(jié)果存在顯著差異，平面方法與球面方法差異較大，差值結(jié)果的標(biāo)準(zhǔn)差分別為1.59 mGal和2.29 mGal，主要原因是兩種方法對地形體的近似模型不同。平面積分方法和平面FFT方法的差值位于-1～6 mGal，標(biāo)準(zhǔn)差為0.78 mGal，不存在較大差別。中部地區(qū)差異較小，在0 mGal左右，南部和北部山區(qū)計算結(jié)果的差異化較大，說明計算點周邊的復(fù)雜地形對改正量造成顯著影響。后續(xù)將進(jìn)一步研究改進(jìn)、提高精度的空間：一是進(jìn)一步提高輸入的地形模型分辨率和精度；二是盡可能擴(kuò)大積分半徑，對較遠(yuǎn)區(qū)地形影響可以采用粗/細(xì)格網(wǎng)相結(jié)合的方式進(jìn)行計算提高效率，改善精度；三是采用嚴(yán)密計算方式，提高平面FFT方法計算階數(shù)。

圖2 局部重力地形改正值

圖3 局部地形改正值結(jié)果比較

表1 局部地形改正值結(jié)果統(tǒng)計信息 mGal

表2 局部地形改正值比較結(jié)果統(tǒng)計信息 mGal

4 結(jié)束語

本文給出了局部地形改正的平面積分方法、平面FFT方法、球面積分方法的具體計算表達(dá)式，以西部山區(qū)為例，計算了重力局部地形改正值，比較了不同結(jié)果的差異。結(jié)果表明：平面積分方法和平面FFT方法計算結(jié)果接近；球面積分計算結(jié)果與平面積分方法和平面FFT方法具有較大差異。三種方法存在差異的本質(zhì)原因是近似面分別作為球面和平面產(chǎn)生的系統(tǒng)性影響不同，在對大面積地形開展積分時，球面更加能貼合地球形狀，精度較高；局部地形改正值大小主要與計算點近區(qū)地形起伏有關(guān)，若近區(qū)地形起伏較大會對改正值產(chǎn)生很大影響，若近區(qū)地形較為平緩則其改正值多在0mGal附近，重力值幾乎不受影響。計算重力局部地形改正時應(yīng)根據(jù)地形起伏、分辨率，計算效率等要求選擇合適的計算方法。在平原地區(qū)，地形起伏變化不大，無論采用球面積分或是平面積分，近區(qū)地形在垂直方向上沒有較大差異，積分結(jié)果主要受算式中常數(shù)因子影響，這種情況考慮平面方法即可。鑒于FFT方法的速度優(yōu)勢，針對大面積的地形改正計算，此方法既能保證精度和球面積分近似，又能提高計算效率；在山區(qū)，球面積分對平面積分有很強(qiáng)的改善作用，為精度考慮，采用球面積分方法為宜。