張 偉 廖國忠 張秋冬 李 華

1 中國地質(zhì)調(diào)查局成都地質(zhì)調(diào)查中心，成都市一環(huán)路北三段2號，610082

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

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

1 球坐標系下的遠二區(qū)地改理論

1.1 球坐標下的遠區(qū)地改計算公式

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

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

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

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

再令

于是，

1.2 地殼扇形“質(zhì)量體”的高程

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

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

M為某環(huán)上的網(wǎng)格數(shù)，而P點遠區(qū)地改值可表示為二重離散積分公式。式（7）可由計算機編程實現(xiàn)：

N為20～166.7km間的環(huán)帶數(shù)。

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

2 計算方法與實現(xiàn)

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

式（7）中，M、N越大，其環(huán)帶和環(huán)帶上的網(wǎng)格剖分越精細，計算結果越逼近解析值，但計算量也越大。區(qū)域重力調(diào)查規(guī)范推薦的分環(huán)和網(wǎng)格參數(shù)如圖3和表1所示。

表1 區(qū)域重力調(diào)查規(guī)范中遠二區(qū)計算的環(huán)帶和方位數(shù)（R＝6 371.025km）Tab.1 The subdivision of zone and azimuth in regional gravity survey specification（R＝6 371.025km）

圖3 區(qū)域重力調(diào)查規(guī)范中環(huán)帶及其網(wǎng)格劃分模型Fig.3 Zone subdivision model and its mesh in regional gravity survey specification

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

聯(lián)立式（9）和（10），可求得網(wǎng)格4個節(jié)點的經(jīng)緯度坐標。當網(wǎng)格較小時，可將4個節(jié)點視為同平面，通過平面插值方法，利用待插值點鄰近四周的DEM 數(shù)據(jù)得到節(jié)點的高程值。將4 個節(jié)點的高程值取平均，便可得到該網(wǎng)格的平均高程。代入式（7），得到測點在球坐標系下的遠區(qū)地形改正值。由于式（9）存在除數(shù)為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。

在有限球殼模型下，將不同的環(huán)數(shù)和方位數(shù)代入式（5），計算值與理論值之間的差值如表2所示?？芍诓豢紤]高程誤差時，式（5）的計算精度主要取決于環(huán)數(shù)，與方位數(shù)無關，環(huán)數(shù)越大，其計算精度越高，當N＞4 000時，式（5）的計算精度達到μGal級。ASTER GDEM V2 數(shù)據(jù)網(wǎng)格30m×30 m，其環(huán)數(shù)最大取值可到5 500，滿足μGal級遠二區(qū)地改要求。

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

3.2 誤差精度

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

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

［1］李東漢.試論區(qū)域重力測量遠區(qū)地改的精度［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］張俊，張寶松，邸兵葉，等.高程數(shù)據(jù)網(wǎng)格間距對重力中區(qū)地形改正精度的影響［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］趙海濤，張兵，左正立，等.中國及周邊區(qū)域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數(shù)據(jù)的中國典型區(qū)域SRTM 與ASTER GDEM 高程精度評價［J］.地球科學—中國地質(zhì)大學學報，2013，38（4）：887－897（Du Xiaoping，Guo Huadong，F(xiàn)an 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單元體泰勒級數(shù)展開方法的比較［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］安玉林，張明華，黃金明，等.純球坐標系內(nèi)各項重力校正值計算方案和計算過程［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］.浙江國際海運職業(yè)技術學院學報，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）