葉周潤(rùn) 柳林濤 梁星輝

1 中國(guó)科學(xué)院大學(xué)，北京市石景山區(qū)玉泉路甲19號(hào)，100049

2 中國(guó)科學(xué)院測(cè)量與地球物理研究所大地測(cè)量與地球動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢市徐東大街340號(hào)，430077

地形和均衡影響的計(jì)算是非常重要的數(shù)據(jù)處理技術(shù)［1－5］。在忽略離心力的微小影響假設(shè)下，地形改正通常通過計(jì)算正演模型引力位的相關(guān)導(dǎo)數(shù)來實(shí)現(xiàn)。一般而言，在局部范圍內(nèi)常選用直角坐標(biāo)系下的棱柱體模型，其引力位及其相關(guān)導(dǎo)數(shù)都有嚴(yán)格的理論表達(dá)式，同時(shí)對(duì)應(yīng)的譜域形式通常采用FFT技術(shù)［6－9］。由于地球形狀接近球體，球坐標(biāo)下的近似模型可以盡量避免平面誤差從而比較滿足研究實(shí)際。在空間域，卡爾斯魯厄理工學(xué)院的研究團(tuán)隊(duì)發(fā)展了一套基于Tessroid 單元體的泰勒展開方法［10－12］。在三維頻譜域，常采用基于單層位技術(shù)的球諧分析方法［13－14］。對(duì)于均衡影響計(jì)算，理論上是壓力平衡原理下地形的補(bǔ)償影響，在計(jì)算方法上和地形正演沒有實(shí)質(zhì)差異［11］。

重力梯度是重力位的二階導(dǎo)數(shù)，能從多角度反映異常體的中心點(diǎn)和邊界等細(xì)節(jié)信息［15］。最新的GOCE模型改進(jìn)了165～220階的地球重力場(chǎng)細(xì)節(jié)，可以捕捉更多的地球內(nèi)部信號(hào)［16］。本文計(jì)算了中國(guó)及其鄰區(qū)地形和均衡的重力梯度影響，并給出均衡重力梯度異常分布圖。實(shí)驗(yàn)中原始重力梯度數(shù)據(jù)來自GOCO03S位系數(shù)模型，且保留了高階GOCE中的高頻數(shù)據(jù)［17］。在地形改正部分，原始數(shù)據(jù)來自DTM2006高程模型［18－19］。在均衡模型選取方面，以Airy－Heiskanen均衡模型起主導(dǎo)作用［3］。

1 計(jì)算方法介紹

1.1 由GOCE模型計(jì)算擾動(dòng)重力梯度值

解算完成的GOCE 模型最終以球諧位系數(shù)形式發(fā)布。如果要得到全球或局部地區(qū)的擾動(dòng)重力梯度值，則需要通過球諧綜合方法。其擾動(dòng)位合成的數(shù)學(xué)基本表達(dá)式如下：

這里，r、θ和λ分別表示計(jì)算點(diǎn)的球心距離、緯度和經(jīng)度；GM代表地球總質(zhì)量和萬有引力常數(shù)的乘積；R是地球平均半徑和表示擾動(dòng)位系數(shù)；是完全正則化的締合勒讓德函數(shù)；Nmin和Nmax分別為擾動(dòng)位球諧表達(dá)的起始和最大階數(shù)，考慮到GOCE模型的最大階次，Nmin和Nmax取值為2 和250。

因?yàn)樗杏?jì)算量是擾動(dòng)位在緯度、經(jīng)度及地心方向的二階導(dǎo)數(shù)，則［17］：

具體表達(dá)式如下：

式中的相關(guān)表達(dá)式詳見文獻(xiàn)［9］。

1.2 球諧譜正演擾動(dòng)重力梯度

在大地水準(zhǔn)面為球面的假設(shè)下，地形單層位3次項(xiàng)球諧譜引力位展開式在形式上和式（1）一致［13］，但為了表達(dá)方便，我們定義：

在具體計(jì)算中，式（1）中的球諧系數(shù)由下式代替：

其中，

這里，i∈｛1 ，2，3｝ ；表示地球平均密度，取值5 500 kg／m3；h2和h1分別表示積分體的高程上下限，具體數(shù)值以球面（大地水準(zhǔn)面）為基準(zhǔn)，球面以上為正，反之為負(fù)（圖1）。這里注意：1）球諧階次n起始項(xiàng)為0?？紤]到地形信號(hào)和和重力模型分辨率的匹配，n最大項(xiàng)截至250。2）從地形擾動(dòng)位計(jì)算各重力梯度分量，可直接采用式（3）。

圖1 地形示意圖Fig.1 Schematic view of topography

1.3 Airy－Heiskanen均衡模型

如圖2，在Airy－Heiskanen（A－H）模型中，地殼漂浮在較重的均質(zhì)巖漿上，較輕的山根補(bǔ)償山體的質(zhì)量過剩，較重的反山根補(bǔ)償海水的質(zhì)量不足。一般設(shè)定地球正常地殼深度D＝30km，補(bǔ)償物質(zhì)和地幔物質(zhì)密度差為600kg／m3。

如果假定大地水準(zhǔn)面為近似球面，根據(jù)等質(zhì)量原則，則（反）山根公式表述如下［11］：這里ρ0和ρw分別為巖石和海水密度，取值2 670 kg／m3和1 025kg／m3。

圖2 Airy－Heiskanen均衡模型Fig.2 Isostatic model of Airy－Heiskanen

2 實(shí)驗(yàn)及結(jié)果分析

實(shí)驗(yàn)區(qū)域?yàn)榫暥?°～59°、經(jīng)度60°～149°的中國(guó)及其鄰近區(qū)域，最終圖形結(jié)果均以1°×1°分辨率展示。實(shí)驗(yàn)中，球諧合成時(shí)的正常重力場(chǎng)模型為GRS80，所有計(jì)算參考高度為GOCE軌道平均高度255km。此外，由于空氣的重力梯度效應(yīng)微小，作忽略處理［14］。

GOCO03S模型恢復(fù)得到的255km 中國(guó)及其鄰近區(qū)域上空的擾動(dòng)重力梯度分布見圖3，其統(tǒng)計(jì)結(jié)果見表1。從圖可見，緯度方向的擾動(dòng)梯度在中國(guó)范圍內(nèi)主要集中在青藏高原及其鄰區(qū)，表現(xiàn)為高原區(qū)域?yàn)樨?fù)、兩側(cè)為正的近似對(duì)稱異常分布。較于前者，經(jīng)度方向的重力梯度異常信息更加豐富，在經(jīng)度90°～110°范圍和太平洋區(qū)域內(nèi)存在顯著的豎向梯度異常帶；在球心方向的垂直梯度圖中，青藏高原和新疆地區(qū)圖形輪廓類似于緯度方向但數(shù)值正負(fù)相反，且在日本和太平洋地區(qū)異常圖形與環(huán)太平地震帶高度吻合。

圖4是基于圖3數(shù)據(jù)經(jīng)過地形改正后的擾動(dòng)重力梯度圖。從圖4相對(duì)于圖3的輪廓變化及表2的極值分布范圍可知，地形產(chǎn)生的擾動(dòng)影響較大，尤其是高山密布區(qū)域的青藏地區(qū)。所以當(dāng)完成地形擾動(dòng)高頻信號(hào)濾波后，會(huì)凸顯地球內(nèi)部構(gòu)造等中低頻信息。圖4中，緯度方向在塔里木盆地、青藏高原和喜馬拉雅山南部附近區(qū)域存在3條明顯的緯度方向梯度帶，同時(shí)在太平洋地區(qū)表現(xiàn)出大片負(fù)異常；相對(duì)前者，圖2（b）中經(jīng)度梯度帶的數(shù)量明顯多于緯度方向，東北－太行－武夷山梯度帶清晰可見，并在經(jīng)度95°左右有一條縱穿中國(guó)大陸的中央梯度帶［2］；垂向?qū)?shù)的異常結(jié)果總體呈現(xiàn)出陸地為正、海洋為負(fù)的圖形分布。

在圖4基礎(chǔ)上經(jīng)Airy－Heiskanen模型改正后的均衡重力梯度異常見圖5，其統(tǒng)計(jì)結(jié)果見表3。均衡異常反映了由地殼運(yùn)動(dòng)導(dǎo)致的靜力平衡偏離程度。因?yàn)榫馐Ш庑枰鄳?yīng)均衡調(diào)整而產(chǎn)生新的構(gòu)造運(yùn)動(dòng)，所以通常重力失衡帶就是活動(dòng)構(gòu)造帶，并且 地 震 多 發(fā) 生 在 附 近 異 常 梯 度 帶 上［3，20，21］。如圖所示，均衡異常圖和GOCE 衛(wèi)星觀測(cè)到的擾動(dòng)異常具有總體相似性。在喜馬拉雅山脈附近地區(qū)，所有重力梯度結(jié)果都顯示該地區(qū)屬于非均衡狀態(tài)。異常梯度帶按照?qǐng)D5（c）顯示，山脈南側(cè)為均衡負(fù)異常而北側(cè)為正異常，可能表示印度洋板塊的碰撞作用導(dǎo)致北部地殼在增厚［22］。此外，在中國(guó)區(qū)域內(nèi)有4個(gè)明顯的均衡異常區(qū)域，從3個(gè)分量圖像及重力梯度的幾何圖形意義［15］可知，這些異?；境尸F(xiàn)點(diǎn)狀分布。對(duì)比文獻(xiàn)［21］的重力均衡異常結(jié)果發(fā)現(xiàn)，在31°N、103°E的汶川地區(qū)附近有兩個(gè)相近的異常區(qū)域，且汶川地區(qū)異常值較大，但經(jīng)過汶川地震能量釋放后該區(qū)均衡異常已相對(duì)減弱。在太平洋區(qū)域，均衡異常圖形顯示的是和環(huán)太平地震帶高度符合的正負(fù)異常相對(duì)的海溝－?；》植肌?/p>

圖3 GOCE高度中國(guó)及其鄰區(qū)擾動(dòng)重力梯度圖Fig.3 The map of gravity gradient disturbance in China and its neighboring area at GOCE altitude

圖4 經(jīng)過地形改正后的中國(guó)及其鄰區(qū)擾動(dòng)重力梯度圖Fig.4 The map of gravity gradient disturbance after topography correction in China and its neighboring area

表1 GOCE高度中國(guó)及其鄰區(qū)擾動(dòng)重力梯度結(jié)果統(tǒng)計(jì)Tab.1 The statistics of gravity gradient disturbance in China and its neighboring area at GOCE altitude

表2 經(jīng)過地形改正后中國(guó)及其鄰區(qū)的擾動(dòng)重力梯度結(jié)果統(tǒng)計(jì)Tab.2 The statistics of gravity gradient disturbance after topography correction in China and its neighboring area

圖5 中國(guó)及其鄰區(qū)均衡重力梯度異常圖Fig.5 The isostatic anomaly map of gravity gradient disturbance in China and its neighboring area

表3 中國(guó)及其鄰區(qū)均衡重力梯度異常結(jié)果統(tǒng)計(jì)Tab.3 The isostatic anomaly statistics of gravity gradient in China and its neighboring area

3 結(jié) 語

本文在球坐標(biāo)下采用球諧譜方法計(jì)算地形和均衡影響，并給出相應(yīng)的重力梯度異常結(jié)果。由于地形和均衡公式都基于球面模型，從而有效避免了大空間范圍計(jì)算時(shí)的平面誤差影響。從地形改正結(jié)果初步分析，較好地符合了一些地質(zhì)構(gòu)造的先驗(yàn)結(jié)論。從均衡異常圖的初步解讀，與地殼運(yùn)動(dòng)和地震活動(dòng)帶有較好的一致性。文中所示梯度分量能從多角度提供原始信息用于結(jié)果分析，但文中沒有考慮均衡模型橫向密度不均，重力梯度在地球物理學(xué)中的應(yīng)用解釋尚未完善，需要在以后研究中加強(qiáng)。

［1］álvarez O，Gimenez M，Braitenberg C，et al.GOCE Satellite Derived Gravity and Gravity Gradient Corrected for Topographic Effect in the South Central Andes Region［J］.Geophysical Journal International，2012，190（2）：941－959

［2］馬宗晉，高祥林，宋正范.中國(guó)布格重力異常水平梯度圖的判讀和構(gòu)造解釋［J］.地球物理學(xué)報(bào)，2006，49（1）：106－114（Ma Zongjin，Gao Xianglin，Song Zhengfan.Analysis and Tectonic Interpretation to the Horizontal Gradient Map Calculated from Bouguer Gravity Data in the China Mainland［J］.Chinese Journal Geophysics，2006 49（1）：106－14）

［3］方盛明，馮銳，田長(zhǎng)征，等.亞洲中部地區(qū)均衡重力異常特征及地震活動(dòng)性［J］.地震學(xué)報(bào)，1997，19（6）：650－654（Fang Shengming，F(xiàn)eng Rui，Tian Changzheng，et al.Isostatic Gravity Anomaly and Seismic Activity in Central Asia［J］.Acta Seismologica Sinica，1997，19（6）：650－654）

［4］張赤軍，方劍，馬宗晉.均衡異常與現(xiàn)代構(gòu)造應(yīng)力場(chǎng)的初步研究——以東海及其東緣海為例［J］.地震學(xué)報(bào)，2003，25（1）：32－39（Zhang Chijun，F(xiàn)ang Jian，Ma Zongjin.Preliminary Study of Isostatic Anomaly and Modern Tectonic Stress Field［J］.Acta Seismologica Sinica，2003，25（1）：32－39）

［5］方劍.中國(guó)及鄰區(qū)均衡重力異常及其地球動(dòng)力學(xué)特征［D］.北京：中國(guó)地震局地質(zhì)研究所，2006（Fang Jian.Isostatic Gravity Anomaly and Its Geodynamic Characters in China and Its Adjacent Regions［D］.Institute of Geology，CEA，2006）

［6］羅志才，陳永奇，寧津生.地形對(duì)確定高精度局部大地水準(zhǔn)面的影響［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2003，28（3）：340－344（Luo Zhicai，Chen Yongqi，Ning Jinsheng.Terrain Effect in the Determination of Local High Precise Geoid［J］.Geomatics and Information Science of Wuhan University，2003，28（3）：340－344）

［7］黃謨濤，翟國(guó)君.利用FFT 技術(shù)計(jì)算地形改正和間接效應(yīng)［J］.測(cè)繪學(xué)院學(xué)報(bào)，2000，17（4）：242－246（Huang Motao，Zhai Guojun.Terrain Correction and Indirect Effects Using FFT Technology［J］.Journal of Institute of Surveying and Mapping，2000，17（4）：242－246）

［8］Tziavos I N，Sideris M G，F(xiàn)orsberg R，et al.The Effect of the Terrain on Airborne Gravity and Gradiometry［J］.Journal of Geophysical Research：Solid Earth（1978 –2012），1988，93（B8）：9 173－9 186

［9］Zhu L.Gradient Modelling with Gravity and DEM［D］.The Ohio State University，2007

［10］Heck B，Seitz K.A Comparison of the Tesseroid，Prism and Point－mass Approaches for Mass Reductions in Gravity Field Modelling［J］.Journal of Geodesy，2007，81（2）：121－136

［11］Wild F，Heck B.Topographic and Isostatic Reductions for Use in Satellite Gravity Gradiometry［C］.VI Hotine－Marussi Symposium on Theoretical and Computational Geodesy，2008

［12］Grombein T，Seitz K，Heck B.Untersuchungen zur Effizienten Berechnung Topographischer Effekte auf den Gradiententensor am Fallbeispiel der Satellitengradiometriemission GOCE［M］.KIT Scientific Publishing，2010

［13］Tsoulis D.A Comparison between the Airy／Heiskanen and the Pratt／Hayford Isostatic Models for the Computation of Potential Harmonic Coefficients［J］.Journal of Geodesy，2001，74（9）：637－643

［14］Novák P，Grafarend E W.The Effect of Topographical and Atmospheric Masses on Spaceborne Gravimetric and Gradiometric Data［J］.Studia Geophysica et Geodaetica，2006，50（4）：549－582

［15］Saad A H.Understanding Gravity Gradients—A Tutorial［J］.The Leading Edge，2006，25（8）：942－949

［16］Hirt C，Kuhn M，F(xiàn)eatherstone W E，et al.Topographic／Isostatic Evaluation of New－Generation GOCE Gravity Field Models［J］.Journal of Geophysical Research：Solid Earth（1978–2012），2012，117（B5）

［17］Mayer－Gürr T，Rieser D，H?ck E，et al.The New Combined Satellite Only Model GOCO03s［C］.GGHS2012，Venice，Poster，2012

［18］張興福，劉成.綜合EGM2008 模型和SRTM／DTM2006.0 剩余地形模型的GPS 高程轉(zhuǎn)換方法［J］.測(cè)繪學(xué)報(bào)，2012，41（1）：25－32（Zhang Xingfu，Liu Cheng.The Approach of GPS Height Transformation Basedon EGM2008 and SRTM／DTM2006.0Residual Terrain Model［J］.Acta Geodaeticaet Cartographica Sinica，2012，41（1）：25－32）

［19］Pavlis N K，Holmes S A，Kenyon S C，et al.An Earth Gravitational Model to Degree 2160：EGM2008［J］.EGU General Assembly，2008：13－18

［20］蔣福珍，張赤軍，倪志宏，等.青藏高原均衡場(chǎng)與非潮汐重力測(cè)量［J］.地球物理學(xué)報(bào)，1989，32（4）：465－469（Jiang Fuzhen，Zhang Chijun，Ni Zhihong，et al.Isostatic Field and Non－tidal Gravity Measurement in Qinghai－Xizang Plateau［J］.Chinese Journal Geophysics，32（4）：465－469）

［21］王懋基，程振炎.均衡異常與地殼結(jié)構(gòu)［J］.地質(zhì)學(xué)報(bào)，1982，56（1）：52－61（Wang Maoji，Cheng Zhenyan.Isostatic Anomaly and Crustal Structure［J］.Acta Geologica Sinica，1982，56（1）：52－61）

［22］劉燕平，蔡少華.青藏地區(qū)的均衡異常研究［J］.地殼形變與地震，1988（4）：8－12（Liu Yanping，Cai Shaohua.Research of Gravity Anomalies in Qinghai－Tibet Plateau［J］.Crustal Deformation and Earthquake，1988（4）：8－12）