重力局部地形改正量的計算方法比較

歐陽明達，周 巍，張敏利

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

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

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

1 平面積分

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

(1)

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

(2)

假設DTM以N×M的規則格網數據表示，每個網格內的地形用一個質量均勻分布而高程不變的棱柱體代替，則得到用質量棱柱地形模型表示的地形改正式：

(3)

令

(4)

(5)

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

(6)

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

(7)

即

δg≈C1+C2+C3.

(8)

本文取其1階近似，得到

(9)

對其進行FFT變換，改寫成譜形式：

(10)

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

2 球面積分

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

(11)

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

(12)

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

(13)

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

(14)

其中，t=cosψ，C為積分常數。

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

(15)

其中，

(16)

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

(17)

3 算例與結果分析

3.1 數據準備與方案設計

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

圖1 地形模型

3.2 結果分析與比較

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

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

圖2 局部重力地形改正值

圖3 局部地形改正值結果比較

表1 局部地形改正值結果統計信息 mGal

表2 局部地形改正值比較結果統計信息 mGal

4 結束語

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