球坐標系下基于全球DEM 數據的重力遠二區高精度地改算法

張 偉 廖國忠 張秋冬 李 華

1 中國地質調查局成都地質調查中心，成都市一環路北三段2號，610082

2 成都理工大學地球物理學院，成都市二仙橋東三1號，610059

隨著重力傳感器制造工藝的不斷發展，μGal級的重力測量工作已成為可能。文獻［1－2］對影響重力地改精度的主要原因進行探討，認為誤差精度取決于數字地形高程與真實地形地貌的逼近程度。數字地形高程的網格距越小，高程值越精確，其地形改正值的計算精度就越高［3－4］。而對于球坐標系下重力異常的計算，文獻［5，6］從不同角度推導了相同的解析式。本文在球坐標系下重力遠區地形改正理論模型的基礎上，提出使用公開的全球高分辨率DEM 數據計算遠二區地改。實例表明，該方法能夠大大提高遠二區地改的計算精度，達到μGal量級。

1 球坐標系下的遠二區地改理論

1.1 球坐標下的遠區地改計算公式

球坐標系下的球體理論模型如圖1所示，地球球心O為原點，N、S點分別為地球北、南極，P、N、S在同一個過球心O的平面上，點M在空間曲面上，以OP的延長線方向為Z軸，OX軸垂直于PNS平面，OY軸垂直于XOZ平面，X、Y、Z軸滿足右手坐標系。α1、α2、α3可看成以P為北極的“經度”，β1、β2為其“緯度”，空間曲面上的點M可表示為Μ（α，β，h）。在球極坐標系下，隨地形起伏的質量體Ω可描述為地殼扇形塊h）dv，其底面是以地球平均半徑R所作的球面，頂面是地形的自由起伏面。當積分區域為Ω＝｛（α，β，h）｜α1≤α≤α2，β1≤β≤β2，0≤h≤H｝時，質量體Ω相對于測點P的地形改正值為：

圖1 球坐標系下的重力遠二區地改模型Fig.1 Far zone II’s terrain correction model in spherical coordinates

式中，R為地球平均半徑，取大地水準面；z、H為P、M點的海拔高程；G為引力常數；ρ為地殼平均密度。

式（1）中的積分函數對變量α而言是常量積分，當α∈［0，2π］時，ΔgΩ變為環帶（β1，β2）相對于測點P的 重力地改值Δg［β1，β2］。將式（1）中的變量h設為常量，取h＝，為積分區域的平均高程，此時式（1）可簡化為［6］：

再令

于是，

1.2 地殼扇形“質量體”的高程

計算遠二區的地形改正值時，采用的高程資料是全國5′×5′高程節點網，以每個地殼扇形“質量體”落入5′×5′高程數據個數相加求和取平均值作為該扇形體的平均高程。但由于這類數據網格距較大，式（2）中與真實地貌間必然存在較大誤差。而全球ASTER GDEM 數據的分辨率能達到3″×3″，從而大大提高遠二區地改的精度。ASTER GDEM 和SRTM3 DEM 數據基于WGS－84坐標系，該坐標系下的重力遠區地改模型如圖2所示。當β∈［β1，β2］，α∈［0，2π］時，M點的“軌跡”是球面上以OP為軸的環帶，可將式（5）離散化成式（6），即將環帶等分為M個網格，M值越大，環帶上的網格劃分越精細，對應式（2）中的越逼近真實地形。當M→∞時，式（6）便變為理論解析式（2）。

圖2 WGS－84坐標下的重力遠區地改模型Fig.2 Far zone II’s terrain correction model in WGS－84coordinates

M為某環上的網格數，而P點遠區地改值可表示為二重離散積分公式。式（7）可由計算機編程實現：

N為20～166.7km間的環帶數。

代入式（7）前的一個關鍵問題是將WGS－84坐標系下的經緯度坐標（θ，Φ）轉換為球坐標下的坐標（α，β）。根據球面三角形公式［7］，兩坐標系之間的轉換關系可由式（8）和式（9）表示。如圖2所示，α的取值以PN方向為起始0°，令方位角按相對于球心O的順時針方向增加。當M點位于平面POY前方時，即θ0－θ＞0，α的取值范圍為180°～360°；當M點位于平面POY后方時，即θ0－θ≤0，α的取值范圍為0～180°。

2 計算方法與實現

目前，基本覆蓋全球的DEM 數據主要有：ASTER GDEM V2數據、SRTM3 V4.1數據、GTOPO30數據以及我國的資源三號測繪衛星立體影像數據等。其中，ASTER GDEM V2數據的分辨率最高（1″×1″，約30m×30m），SRTM3次之（3″×3″，約90m×90m）。兩種數據歷經幾次版本升級后，質量日趨精良，絕對高度誤差小于16m，相對高度誤差小于10m，滿足重力遠區地改的高程精度要求。

式（7）中，M、N越大，其環帶和環帶上的網格剖分越精細，計算結果越逼近解析值，但計算量也越大。區域重力調查規范推薦的分環和網格參數如圖3和表1所示。

表1 區域重力調查規范中遠二區計算的環帶和方位數（R＝6 371.025km）Tab.1 The subdivision of zone and azimuth in regional gravity survey specification（R＝6 371.025km）

圖3 區域重力調查規范中環帶及其網格劃分模型Fig.3 Zone subdivision model and its mesh in regional gravity survey specification

M、N和網格剖分方案確定后，即可根據α、β計算出某一網格4個節點在WGS－84坐標系下的經緯度值。根據球面三角形邊余弦定理：

聯立式（9）和（10），可求得網格4個節點的經緯度坐標。當網格較小時，可將4個節點視為同平面，通過平面插值方法，利用待插值點鄰近四周的DEM 數據得到節點的高程值。將4 個節點的高程值取平均，便可得到該網格的平均高程。代入式（7），得到測點在球坐標系下的遠區地形改正值。由于式（9）存在除數為0時的奇異值，因此在程序編寫時需要進行特殊處理，當cosΦ＝0時，Φ＝90°，M點處于北極或者南極，θ＝0°；當cosΦ0＝0時，Φ0＝90°，P點處于北極或者南極，θ＝α。

3 算例分析

3.1 理論模型

選用的理論模型為高程＝1 000m 的有限球殼，且計算點高程z＝1 000m。將、z、β1（β1＝20km／R）、β2（β2＝166.7km／R）代入式（2），可求得有限球殼模型的理論地改值為－3.749 57mGal。

在有限球殼模型下，將不同的環數和方位數代入式（5），計算值與理論值之間的差值如表2所示。可知，在不考慮高程誤差時，式（5）的計算精度主要取決于環數，與方位數無關，環數越大，其計算精度越高，當N＞4 000時，式（5）的計算精度達到μGal級。ASTER GDEM V2 數據網格30m×30 m，其環數最大取值可到5 500，滿足μGal級遠二區地改要求。

表2 有限球殼理論模型下不同環數及方位的試算結果對比Tab.2 Comparative results within different subdivisions in limited spherical shell model

3.2 誤差精度

為進一步分析高程誤差對本文計算方法的影響，在上述有限球殼模型（＝1 000m，z＝1 000 m）基礎上，在參與計算的每個高程網格節點上加入高斯隨機噪聲進行試算，高斯噪聲的均值為0，標準差為20，其噪聲水平與ASTER GDEM V2 DEM 的數據質量相當，計算結果如表3所示。從表3可知，在該噪聲水平下，加入噪聲和沒有加入噪聲時的計算結果差值在μGal量級以下，可以忽略。因此，ASTER GDEM V2和SRTM3數據的質量滿足于μGal級遠二區地改的工作要求。如表3第10行所示，當N＝5 500、M＝35 000時，有噪聲與無噪聲時的計算結果差值小于μGal數量級。

表3 有限球殼理論模型在高斯噪聲條件下的試算結果對比Tab.3 Comparative results within different subdivisions under Gauss noise environment

［1］李東漢.試論區域重力測量遠區地改的精度［J］.物探與化探，1984，8（2）：99－104（Li Donghan.A Discussion on the Accuracy of Far Sector Topographic Correction in Regional Gravity Survey［J］.Geophysical ＆ Geochemical Exploration，1984，8（2）：99－104）

［2］張俊，張寶松，邸兵葉，等.高程數據網格間距對重力中區地形改正精度的影響［J］.物探與化探，2014，38（1）：157－161（Zhang Jun，Zhang Baosong，Di Bingye，et al.The Effect of the Grid Spacing of Elevation on the Accurcy of Median Region Terrain Correction of Gravity［J］.Geophysical ＆Geochemical Exploration，2014，38（1）：157－161）

［3］趙海濤，張兵，左正立，等.中國及周邊區域ASTER GDEM 與SRTM DEM 高程對比分析及互補修復［J］.測繪科學，2012，37（1）：8－11（Zhao Haitao，Zhang Bing，Zuo Zhengli，et al.Elevation Comparison and Complementary Repair of ASTER GDEM and SRTM DEM in China and Surrounding Areas［J］.Science of Surveying and Mapping，2012，37（1）：8－11）

［4］杜小平，郭華東，范湘濤，等.基于ICESat／GLAS數據的中國典型區域SRTM 與ASTER GDEM 高程精度評價［J］.地球科學—中國地質大學學報，2013，38（4）：887－897（Du Xiaoping，Guo Huadong，Fan Xiangtao，et al.Verical Accuracy Assesment of SRTM and ASTER GDEM over Typical Regions of China Using ICESat／GLAS［J］.Earth Science－Journal of China University of Geosciences，2013，38（4）：887－897）

［5］杜勁松，陳超，梁青，等.球冠體積分的重力異常正演方法及其與Tesseriod單元體泰勒級數展開方法的比較［J］.測繪學報，2012，41（3）：339－346（Du Jinsong，Chen Chao，Liang Qing，et al.Gravity Anomaly Calculation Based on Volume Integral in Spherical Cap and Comparision with the Tesserio－Taylor Series Expansion Approach［J］.Acta Geodaetica et Cartographica Sinca，2012，41（3）：339－346）

［6］安玉林，張明華，黃金明，等.純球坐標系內各項重力校正值計算方案和計算過程［J］.物探與化探，2010，34（6）：697－705（An Yulin，Zhang Minghua，Huang Jinming，et al.The Computation Scheme and Computation Process for Gravity Correction Values Within the Pure Spherical Coordinate System［J］.Geophysical ＆Geochemical Exploration，2010，34（6）：697－705）

［7］周雪娟，褚善東.球面三角形公式及其應用［J］.浙江國際海運職業技術學院學報，2008，4（1）：59－62（Zhou Xuejuan，Chu Shandong.Formula of Spherical Triangle and Its Application［J］.Journal of Zhejiang International Maritime College，2008，4（1）：59－62）