（22）

利用一维FFT及卷积定理进而得到

（23）

（24）

式中表示逆FFT变换。

对于一维球面FFT算法而言，在每一纬度圈进行FFT计算，而在经度方向上进行求和运算，此时与核函数可看做两个离散的实数序列，序列长度是同一纬度圈被计算点的格网数n。此时如果直接进行FFT计算，则如前所述是不准确的，为了解决这一问题，首先对核函数进行分析，对于同一纬度圈而言，在解析公式计算中核函数实际组成了一个n维的对称阵或反对称阵,而此时序列组成一个观测矩阵。

（25）

而对于一个卷积计算，如下式

（26）

上式中A、a表示两个序列，其中（）

如果由序列a组成了一个循环矩阵C，见下式

（27）

则循环矩阵的特点可知，

（28）

上式标明，如果核函数组成的矩阵是一个循环矩阵，则可以直接利用一维球面FFT算法获得与原解析表达式同样精度的结果。对于全球海洋扰动重力反演而言，将观测得到的海面高数据或垂线偏差数据补充至全球区域，则可满足上述条件进而可利用式（28）进行计算。

对于局部区域而言，由于观测数据覆盖地理区域有限，核函数不能构成循环矩阵，则需要将核函数序列进行改造，增加n-2个数，然后变成新的序列，具体如下。

（29）

同时也增加n-2个零值，形成新的序列，则

   （30）

其中，分别是、组成的矩阵，而重要的是的前n×n个元素正好构成矩阵，因此的前n个序列既是我们所要求的。

5 全球及局部海洋扰动重力反演试验

首先进行全球海洋区域扰动重力反演，使用数据是由EGM2008地球重力场模型[17]（Holmes et al，2008）生成2.5′分辨率的大地水准面高数据，经度范围是0°-360°，纬度范围是90°S-90°N。使用的常参数如下表。

表1 计算用常参数

Tab.1 Constant parameter in computation

地心引力常数GM	3986004.418´108m3s-2
长半轴a	6378137.0 m
正常重力	9.8m/s2
地球平均半径R	6371000.0m

利用大地水准面数据按照式（4）、式（28）反演获得全球海洋区域扰动重力，见图1。在2.5′分辨率大地水准面高数据基础上，生成全球2.5′分辨率的垂线偏差，而后按照式（7）、式（28）反演获得全球海洋区域扰动重力，见图2。 

 

图1 由大地水准面高反演的全球海洋区域扰动重力

Fig.1 Global ocean disturbing gravity derived from geoid height 

 

 

图2 由垂线偏差反演的全球海洋区域扰动重力

Fig.2 Global ocean disturbing gravity derived from vertical deflection 

 

其次进行局部区域扰动重力计算，使用的数据是EGM2008地球重力场模型生成的2.5′分辨率大地水准面及垂线偏差数据，数据范围是东经110°-东经120°，纬度范围是北纬10°-北纬20°。采用的方法同上，在实际计算时采用了式（30）的改化方法，最终得到的局部区域海洋扰动重力见图3、图4。 

 

图3 由大地水准面高反演的局部海洋区域扰动重力

Fig.3 Local ocean disturbing gravity derived from geoid height 

 

 

图4 由垂线偏差反演获得的局部海洋区域扰动重力

Fig.4 Local ocean disturbing gravity derived from vertical deflection 

 

通过比较分析，采用本文提出的精确快速算法可以得到与解析算法同样的结果，且计算速度提高20倍左右。为了对两种反演方法进行对比，以局部区域反演为例对由大地水准面和垂线偏差分别计算的扰动重力进行比较，结果见下表：

表2 两种方法获得的局部区域扰动重力比较结果（mGal）

Tab.2  The comparison results derived from two methods(mGal)

比较范围	最大差值	最小差值	差值均值	标准差
10°×10°	14.66	-17.41	-0.43	0.79

从图1-图4以及表2可以看出，由大地水准面和垂线偏差分别反演获得的扰动重力其差异在0.8mGal以内，这说明两种反演方法在实际计算中是基本一致的。

6 结语

本文对海洋扰动重力反演理论及实现方法进行了深入研究，得到了以下结论。

<!--[if !supportLists]-->1） <!--[endif]-->从经典边值问题理论及球谐函数理论出发推导获得了由大地水准面高以及垂线偏差计算扰动重力的解析计算公式，解决了利用卫星测高数据反演海洋扰动重力的理论问题，也为扰动重力数据的来源及第二大地边值问题的应用提供了解决途径。

<!--[if !supportLists]-->2） <!--[endif]-->提出了海洋扰动重力的精确快速算法。对于全球海洋扰动重力计算，通过将观测数据补充至全球区域可在每一纬度圈进行FFT计算，而后在经度方向上进行求和运算，这种处理方式保证了FFT计算后的结果与原解析算法是一致的。对于局部海洋扰动重力计算，通过对核函数和观测数据补充一定量数据从而使得核函数构成了循环矩阵，进而利用FFT计算后保证了结果与原解析算法是一致的。利用本方法，海洋扰动重力的计算速度相比较解析方法提高约20倍，且避免了海洋重力场反演中存在的混叠、边缘效应问题，计算精度等同于原解析算法，更加重要的是对观测数据的序列长度没有硬性要求，使得该方法的应用更加灵活。

<!--[if !supportLists]-->3） <!--[endif]-->利用本文提出的反演理论和方法计算获得了全球及局部海洋区域的扰动重力数据，并对两种方法进行了比较分析。本文提出的精确快速算法也可以应用于物理大地测量中大部分积分问题的求解如求解大地水准面、垂线偏差、地形改正等。

 

参考文献

 

[1]管泽霖,管铮,黄谟涛等.局部重力场逼近理论和方法[M]. 1997.北京：测绘出版社

[2]PETR V ,ZHANG Changyou ,LARS E S. A Comparison of Stokes and Hotine’s Approaches to Geoid Computation[J]. Manuscripta Geodaetica , 1992 . 17 : 29-35.

[3]Sten Johan ClaessensSolutions to the Ellipsoidal Boundary Value Problems for Gravity Field Modelling[D]. .2006. Curtin University of Technology 

[4]吴晓平,李姗姗,张传定.扰动重力边值问题与实际数据处理的研究[J].武汉大学学报·信息科学版, 2003.28:73-76

[5]李斐,陈武, 岳建利. GPS在物理大地测量中的应用及GPS 边值问题[J].测绘学报, 2003.32(3):198-201

[6] David T.Sandwell,Walter H.F.Smith.marine gravity anamoly from Geosat and Ers-1 satellite altimetry[J].Journal of Geophysical Research.102,B5:10039-10054

[7] C.Hwang. Inverse Vening Meinesz formula and deflection-geoid formula: applications to the predictions of gravity and geoid

over the South China[J]. Sea journal of geodesy (1998)72:304-312

[8]Colombo O L. Numerical Methods for Harmonic Analysis on the Sphere[R]. OSU Report No. 310. 1981.Ohio: The Ohio State University

[9]Forsberg R,Sideris M G. Geoid Computations by the Multi_banding Spherical FFT Approach[J]. Manuscripta Geodaetica, 1993.18(2):82-90

[10]Haagmans R,de Min E,Gelderen M V. Fast Evaluation of Convolution Integrals on the Sphere Using 1D FFT, and a Comparison with Existing Methods for Stokes’Integral[J].Manuscripta Geodaetica, 1993.18(5):227-241

[11]K.P.Schwarz,M.G.Sideris,R.Forsberg. The use of FFT techniques in physical geodesy[J].Geophys.J. Int. 1990. (100) : 485-514.

[12]李建成,陈俊勇,宁津生等.地球重力场逼近理论与中国2000似大地水准面的确定[M]. 2003.武汉:武汉大学出版社 

[13]黄谟涛.利用卫星测高数据反演海洋重力异常研究[J]. 2001.测绘学报, 30(2):179-183.

[14]李建成.卫星测高在大地测量学中的应用及进展[J].测绘科学, 2006.31(6):19-21.

[15]黄谟涛,翟国军,管铮.海洋重力场测定及其应用[M]. 2005.北京：测绘出版社

[16]王冰,申卫昌,田来科等.快速傅立叶变换Cooley-Tukey算法补零问题[J]. 西北大学学报(自然科学版), 2004. 34(1):31-34

[17]Holmes,Steve C. Kenyon, and John K. Factor. 2008.An Earth Gravitational Model to Degree 2160: EGM2008.

EGU General Assembly 2008 Vienna, Austria, April 13-18

上一篇：全球IGS站周年非线性运动振幅分布规律研究

下一篇：数据空白填充对于抑制向下延拓的边界效应…