海洋重力场反演频域形式的探讨

翟振和， 李　冬

1.西安测绘研究所，陕西 西安，710054；2.地理信息工程国家重点实验室，陕西 西安，710054;3.中国洛阳电子装备试验中心，河南 洛阳，471003



海洋重力场反演频域形式的探讨

翟振和1,2， 李冬3

1.西安测绘研究所，陕西 西安，710054；2.地理信息工程国家重点实验室，陕西 西安，710054;3.中国洛阳电子装备试验中心，河南 洛阳，471003

本文针对海洋重力场反演的频域方法的严密性问题，从理论和实践两个角度进行了分析和计算。理论分析表明，Sandwell提出的海洋重力场反演的频域方法本质上是求解海洋扰动重力而非重力异常，该方法在理论推导过程中有多项不严密之处，其主要原因是在推导过程中，将扰动重力的径向导数进行简化，结果忽视了扰动重力的三维空间特征。通过在南海区域与DTU10重力数据比较，表明纯粹基于频域方法的海洋扰动重力精度较之于解析方法约低6mGal，进一步说明了频域方法本身存在的问题。考虑到反演的严密性，建议以解析理论公式结合快速计算策略作为反演海洋扰动重力的首选方法。

海洋重力场反演；频域；扰动重力

1　引　言

海洋重力场反演(包括重力异常和扰动重力)是卫星测高技术领域的重要研究内容，对于地球重力场的精化和重力场模型的构建具有重要意义。Molodensky(1960)给出了反演重力异常的逆VeningMeinesz公式[1]；Smith、Rapp提出了利用最小二乘配置法推估海洋重力异常的理论和方法[2]；黄金维(1998)提出了利用垂线偏差反演重力异常的新的核函数和解算方法，该核函数与Molodensky提出的核函数本质上是等价的[3]。基于上述理论成果，国内外学者利用多代卫星测高数据获得了全球及局部区域的重力异常数据[4-9]。翟振和等(2015)推导获得了海面高及垂线偏差反演扰动重力的解析公式，且理论分析表明，扰动重力不仅在重力场逼近方面有一定优势，同时在海底地形反演中也比重力异常更加严密[10-11]。与上述基于经典物理大地测量理论获得的成果不同，Haxby等(1983)、Sandwell(1997)从频域中的拉普拉斯方程变换导出了利用垂线偏差反演重力异常的频域计算公式[12]，文献[13]引用该公式并将其变换为沿轨迹垂线偏差计算沿轨迹重力异常公式。上述频域公式在形式上与其它反演理论模型存在较大差异，为了探寻这种差异究竟是表面现象还是客观存在，本文将从理论分析和实际计算两个方面进行深入分析与比较。

2　海洋扰动重力反演的频域形式

为了完整描述，本文首先按照Sandwell(1997)的文献对海洋扰动重力反演的频域形式进行重新梳理。考虑到扰动位T是调和函数，因此，T满足Laplace方程，见下式：

(1)

在空间直角坐标系下，扰动重力δg及垂线偏差的东西分量η、南北分量ζ可表示如下。

(2)

(3)

(4)

γ表示一点的正常重力，将以上三式代入式(1)得到：

(5)

考虑到傅里叶变换在求导数时的性质，对式(5)进行二维傅里叶变换后可得到：

(6)

u、v分别表示x、y方向上的空间频率，F2表示二维傅里叶变换。

接着，Sandwell并未直接引用，给出如下公式：

F2(δg(z))=F2(δg(z=0))exp(-2π|k|z)

(7)

将式(7)对z求导(并令z=0)得到

(8)

(9)

上式即利用垂线偏差反演扰动重力的频域形式，通过逆傅里叶变换后可得到扰动重力的空域形式。

3　理论分析

本节首先从理论上对上述推导过程进行探讨分析。从公式上看，式(1)是严密的，式(2)～(5)是将扰动位径向导数用空间直角坐标系中的z方向代替，在局部区域这种近似是可以接受的。对式(5)进行二维傅里叶变换没有问题，但要注意，扰动重力、垂线偏差是一个三维变量，而式(6)只是在(x,y)方向上的二维变换。

式(7)在文献[2]中被描述为由傅里叶变换平移定理得到，但傅里叶平移定理的具体形式如下：

Ft(t+t0)=Ft(t)exp(iλt)

(10)

因此，式(7)显然不是从平移定理得到，而很可能是基于重力数据在空间呈指数衰减的简化公式得到。实际上，从经典理论分析，扰动重力在空域是一个复杂的函数关系，其具体结论可详见文献[14]中解析延拓解中重力异常导数的推导。式(9)的推导实际上忽略了一个过程，即首先使式(6)在z=0条件下得到新的公式如下：

(11)

式(11)中，垂线偏差是在z=0条件下的垂线偏差，而后将式(8)和式(11)比较得到式(9)。

通过以上分析可以看出，式(7)是推导过程中较大的一个近似。此外，基于二维傅立叶变换的推导过程虽没有错误，但其最终的推论式(9)却并不能反映扰动重力三维空间的特征。式(9)表明，地面扰动重力的频域形式只与垂线偏差东西和北南分量及正常重力值相关；而从经典理论推导获得的解析公式中，扰动重力则与地理区域、垂线偏差以及垂线偏差之间的方位角、球面角距相关，如文献[2]给出的公式。

(12)

式中，

(13)

通过以上对比分析，可以看到Sandwell推导扰动重力频域形式过程存在以下几个问题。首先，扰动重力频域公式只与垂线偏差观测数值有关，而与观测数据的地理位置以及观测数据的空间分布毫无关系，这与解析公式明显存在不同；其次，扰动重力频域公式是在一个特殊条件下得到的特定解，若空域中没有得到垂线偏差反演扰动重力的解析公式，则频域公式可作为一种手段以弥补反演的可行性，但在解析公式日趋完善时，频域公式则明显表现不够严密；最后，虽然频域方法易于使用快速傅里叶变换算法，但在空域中一样可以利用快速傅里叶变换进行求解，而且更为严密，这部分资料可见参考文献[15]。

4　计算分析

为了进一步验证频域公式与经典解析公式的差别，本文选用丹麦科技大学发布的海面高数据(DTU10)作为计算用的大地水准面数据，在此基础上分别生成垂线偏差东西和北南分量，而后按照两种反演方法进行计算分析。计算区域范围是(0°N～5°N,110°E～115°E)，格网分辨率为1′，计算结果见图1和图2。选取丹麦科技大学发布的1′分辨率DTU10重力数据作为评价标准，统计结果见表1。

图1　空域解析反演扰动重力

图2　频域反演扰动重力

表1两种方法计算结果的差异见下表(mGal)

类型最大值最小值差值均值标准差解析法70.77-157.93-32.3321.22频域法19.75-232.77-49.5527.25两者互差120.11-61.9417.2128.71

从表1可以看出，两种方法计算结果差异较大，且频域方法精度比解析法约低6mGal。为了说明本文所用解析法的正确性，将数据范围扩大至(-20°S～35°N,90°E～135°E)，而计算区域范围保持不变，计算结果(见图3)与DTU10模型进行比较，结果见下表。

表2解析法计算结果与DTU10模型的差异(mGal)

最大值最小值差值均值标准差19.74-48.72-2.973.52

图3　积分范围扩大后的解析方法结果

从表2可以看出，在仅仅使用平均海面高模型数据的条件下，解析法结果已接近于DTU10模型结果，如果使用沿轨海面高数据并扣除海面地形影响，结果则更加吻合。

5　结　论

本文从理论和实践两方面分析了垂线偏差反演扰动重力(径向)的频域公式，通过分析及实际计算表明，Sandwell实际提出的是径向扰动重力反演频域公式，并且在理论推导过程中进行了若干近似并只考虑二维平面的特殊条件，而忽略了扰动重力的三维空间特征。相比较而言，基于快速傅里叶变换的解析算法仍是目前精确反演海洋重力场的主要方法。

目前，在科技论文出版中会出现这样一种现象，一个科学理论、模型在某个特定环境、特定条件下产生，且在产生过程中采用了一些近似、假设，而以上种种限制却并不会明显地在其成果(如论文、报告)中体现，此时，论文中的成果就往往会给后来者带来误导。在这一点上，钱学森曾指出“科学文章中的错误必须及时阐明，以免后来工作者误用而误事”，他本人也以身作则改正了发表论文中的错误之处[16]。在Sandwell提出频域解算方法之后，鲜有文献对其质疑，在关于海洋重力场反演的诸多文献中，也只有文献[1]中提到“垂线偏差反演扰动重力的频域公式有一定近似性、不甚严密”。我们追求科学的本质是无限接近地认知这个世界，如果不能严密、谨慎地考虑问题，那么我们的认知则无法接近事实。

[1]Molodensky M S, Eremeev V F,Yrukina M J.Methods for Study of the External Gravitational Field and Figure of the Earth Works of Central Research Lnstitute of Geodesy Aerial Photography and Cartography[M].No131.Moscow:Geodetic Literature Press,1960.

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

[3]C. Hwang.Inverse Vening Meinesz Formula and Deflection-Geoid Formula:Applications to the Predictions of Gravity and Geoid Over the South China Sea[J]. Journal of Geodesy,1998(72):304-312.

[4]许厚泽,王海瑛,陆洋等.利用卫星测高数据推求中国近海及邻域大地水准面起伏与重力异常研究[J].地球物理学报，1999，42(4)：465-468.

[5]翟国君, 黄谟涛,谢锡君等.卫星测高数据处理的理论与方法[M].北京:测绘出版社，2000.

[6]李建成,宁津生,陈俊勇等．中国海域大地水准面和重力异常的确定[J]．测绘学报，2003，32(2)：114-119.

[7]王虎彪,王勇,陆洋.联合多种测高数据确定中国海及其邻域1.5′重力异常[J]. 武汉大学学报·信息科学版，2008(33):1292-1295.

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

[9]黄谟涛,王瑞,翟国君.多代卫星测高数据联合平差及重力场反演[J]. 武汉大学学报·信息科学版,2007,32(11):988-993.

[10]翟振和.海洋测高卫星数据处理理论与应用方法研究[D].郑州：信息工程大学，2015.

[11]翟振和,孙中苗,王兴涛.全球及局部海洋扰动重力反演的快速解析方法[J].测绘学报，2015,44(8):827-830.

[12]David T.Sandwell,Walter H.F.Smith.Marine Gravity Anomaly from Geosat and ERS-1 Satellite Altimetry[J].Journal of Geophysical Research,1997(10):039-054.

[13]王海瑛,王广运.卫星测高数据的沿轨迹重力异常反演法及其应用[J]. 测绘学报，2001,30(1):21-27.

[14]翟振和,王兴涛,李迎春.解析延拓高阶解的推导与比较分析[J].武汉大学学报·信息科学版,2015,40(1)：134-139.

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

[16]叶永烈.走近钱学森[M].上海：上海交通大学出版社，2009.

Discussion about the Frequency Domain Form of Ocean Gravity Field Recovery

Zhai Zhenhe1,2, Li Dong3

1.Xi’an Research Institute of Surveying and Mapping, Xi’an 710054， China 2.State Key Laboratory of Geo-information Engineering, Xi’an 710054, China 3.Luoyang Electronic Equipment Test Center of China, Luoyang 471003, China

The rigor of the frequency domain method of ocean gravity field recovery is analyzed and calcultaed from the theoretical and practical aspects. The theoretical analysis shows that the frequency domain method of ocean gravity field recovery put forward by Sandwell is intended to calculate the disturbing gravity rather than gravity anomaly. Besides, several imprecise derivation exist in the theory because it only considers the special conditions in 2D plane but neglects the 3D spatial characteristics of disturbing gravity. The practical calculation shows that there are great differences between the ocean disturbing gravity results calculated only based on the frequency domain with this method, which also demonstrates that the frequency domain method has problems in itself. Considering the rigor of the ocean gravity field recovery, it is suggested to combine the calculation theory formula with rapid calculation strategy to recover the ocean disturbing gravity.

ocean gravity field recovery; frequency domain; disturbing gravity

2015-12-14。

高分专项青年创新基金资助项目(GFZX04060103-7-14)。

翟振和(1980—)，男，助理研究员，主要从事物理大地测量研究。

P223

A