2010, Vol. 53
2. 61363部队, 西安 710054
2. 61363 Troops, Xi'an 710054, China
高频重力场信息是制约精密重力场模型构建的关键[1, 2].要达到中国地区5′×5′网格平均重力异常精度±3×10-5m·s-2，部分地区1′×1′网格平均重力异常精度±5×10-5m·s-2，除引进新型测量技术外，利用精细高程数据和有限重力测量点值构制全国5′×5′和1′×1′网格平均重力异常技术也至关重要[3].
传统5′×5′网格平均重力异常主要是利用布格异常间接推估而得出的.众所周知，由于去掉了布格层物质，在山区布格异常的变化幅度比空间异常的变化幅度要小得多[4, 5]，因此运用布格异常按间接内插法推求空间点的布格异常值可以提高精度.但这一方法没能充分利用地形信息，在地形激烈起伏的大山区，当实测点密度较小时，也能产生不能容忍的误差.产生这种现象的主要原因是由于层间改正仅消去了物质层的引力，它把点周围认为是平坦的，且把物质层的密度看作是相同的，因而造成了误差.为了减小这部分误差，在运用间接内插法时要把地形改正[6]考虑进去，以提高网格平均重力异常的分辨率及推值精度.本文在传统空间异常和布格异常的基础上，以相对平滑的完全布格异常和完全空间异常作为研究对象，统计计算它们的特征参数，这些参数是进行重力异常的推值、提高插值精度以及应用最小二乘配置技术的关键之一[7~10].
2 重力异常的统计特征参数重力异常的协方差函数是进行重力异常统计分析的基本依据，如果协方差函数正确地反映了本区域重力异常变化规律和相关性，那么由此可以导出比较好的结果，反之就会歪曲结果.
2.1 局部重力异常协方差模型经验协方差函数应该是相关函数.设l为距离，从数学角度分析，由相关函数的定义，协方差函数应当满足：（1）当点的距离为0时，相关性为正，即C（0）=C0＞0；（2）协方差函数应为偶函数，即C（-l）=C（l）；（3）当点的距离为0时相关性最强，即C（0）≥C（l），当且仅当l=0时等号成立.同时满足这样三个条件的经验协方差函数为强相关协方差函数.一般研究中都选择“强相关协方差函数”作为局部重力异常经验协方差函数.本文选择局部经验空间调和协方差公式[2, 11]作为统计的经验协方差函数，如下式所示：
|
(1) |
调和扩展到Z＞0的外部空间：
|
(2) |
式中C0为重力异常方差，B是模型参数，为待求参数，b=1/B.
局部重力异常协方差模型可以用三个基本参数来表征[11~13]：重力异常方差C0、相关长度S以及曲率参数χ.其中C0的大小反映了本区域内重力异常变化的幅度；相关长度S是当协方差值为1/2C0时的l值，体现了本区域内重力异常的相关程度；曲率参数χ表明了本区域重力异常变化的激烈程度，它是一个与协方差函数曲线在l=0处的曲率τ0有关的量，即χ=τ0/C0·S2.
当选定了协方差函数模型（1）式时，我们可以得出：
|
(3) |
|
(4) |
以面积内任一点的重力异常代表该面积平均重力异常所产生的中误差称之为代表误差[11].若用
|
(5) |
式中N为网格内实测重力点数，Δgi为网格内第i个重力点的空间异常，
|
(6) |
代入（5）式，代表误差又可表示为
|
(7) |
根据亨特代表误差经验公式[4]
|
(8) |
因此代表误差系数c可表示为
|
(9) |
统计计算时，统计区域大小大致选为大于10km×10km，小于50km×50km.由（7）式可看出，只有当网格内重力点数N较大时算出的代表误差值才可靠.
3 统计计算分析 3.1 重力异常协方差函数模型统计计算在统计分析前，我们均将统计区域按所需分辨率处理成了矩形网格的形式，假定横向网格数N个，纵向网格数M个.
3.1.1 分析统计过程（1）计算归中后的异常值Δg i，j*将选择的统计区域内的所有格网Δgi，j值取中数，作为区域平均异常值，减去这个平均值以得到归中后的新异常值Δg i，j*，即
|
(10) |
（2）计算重力异常方差C0
用上面计算的归中后的异常值Δg i，j*作为计算协方差值的起始观测数据，方差C0的计算公式为：
|
(11) |
（3）统计实际的协方差
若以q表示正整数1、2、3、4、……，相应距离的协方差值可以分别按纵横两方向计算，计算公式为：
|
(12) |
|
(13) |
把所有具有系统距离l值的纵横方向算得的协方差值取中数作为C（lq）的统计数值：
|
(14) |
式中，C（lq）φ为纬向协方差值，C（lq）λ为经向协方差值.
（4）拟合协方差
利用已知的众多距离l和与l相对应的C（lq），通过恒等变换，把方程（1）进行线性化，依据最小二乘原理解得协方差系数B，进而求得相关长度S.
3.1.2 实验计算本次初步实验工作以某区大约1°×1°范围的12569个离散点进行统计计算.在试验区内，实测重力点的分布很不均匀，距离完全相等或接近相等的对点很少存在，区域中部存在大量的重复观测，因此在数据统计前先将重复观测数据进行了整理，认为距离相差300 m的点为重复观测点，且权值相等，用求和取平均的方法进行合并.然后对重力异常数据进行网格化，网格大小宽度为1km×1km.经整理、网格化处理后的试验区高程、空间异常、布格异常如图 1所示.
|
图 1 试验区高程(m)(a)、空间异常(10-5m. s-2)(b)和布格异常(10-5m. s-2)(c)(横轴表示经度方向, 纵轴表示纬度方向) Fig. 1 Height of test area (m)(a), free air anomalies (10-5m. s-2)(b) and Bouguer anomalies (10 -5m. s-2)(c) |
从图 1c可以明显看出，在试验区内存在有数据粗差点，在进行统计分析前先剔除了12090~12093号、11309和11312号六个粗差点，整个试验区划分成50km×50km面积范围大小的4个子区，分别进行重力异常的统计计算.采用地形改正快速计算[15, 16]软件进行了相应区域的地形改正插值计算，将试验区内12563个点的空间异常换算成完全空间异常，并进行统计.为保证统计工作的可靠性，只对25km以内距离上的重力异常协方差进行统计，不足整数时采取四舍五入的原则.表 1是根据统计计算给出的试验区4个子区不同长度时协方差的函数值.
|
表 1 试验区不同距离(l)上的完全空间异常协方差统计计算值 Table 1 Statistics of Faye anomaly covariance in different distance of test area |
图 2圈点线为根据统计计算值刻画出试验区内4个子区域的协方差特征参数模型，叉虚线是根据计算值按照最小二乘平差的原则，拟合出这4块区域的协方差模型.表 2给出了试验区4个子区的协方差统计参数值.
|
图 2 试验区4个子区的完全空间异常的统计协方差与拟合协方差模型 Fig. 2 Fitting curves and statistical curves of Faye anomaly covariance model in the four sub-areas |
|
|
表 2 试验区协方差模型统计参数 Table 2 Statistical parameters of covariance in test areas |
从图 2可以看出，协方差函数所表现出的数据点间的相关性随距离的增大而迅速减小，减小的速度在各个地区也不相同.在不同地区，实际统计的协方差值与统计计算出的协方差函数值有一定的差异，在一些地区统计值与函数值符合得较好，而另一些地区符合较差.
3.1.3 各类地区重力异常协方差函数的统计在试验试算的基础上，把计算实施方案应用到全国各类典型地区空间异常、布格异常、完全空间异常、完全布格异常的协方差函数的统计计算中，统计区域都是1°×1°的范围，异常值是利用实际测量的重力点数据计算的1′×1′平均重力异常值.由于重力异常协方差模型函数应该是正定函数，统计中也是采用空间调和正定函数.但在实际统计计算中，统计到一定的点距上所选择统计区域的协方差函数值可能出现负值，在这种情况下，本文不统计负值部分，仅统计正值部分.在实际统计中，共统计了不同典型地形区域17个1°×1°的空间异常和布格异常协方差，统计情况见表 3、4；22个不同典型地形1°×1°的完全空间异常和完全布格异常协方差，统计情况见表 5、6；表 7是统计地区的地形类别情况.
|
|
表 3 空间异常协方差模型特征参数统计表 Table 3 Statistics of characteristic parameters of free air anomaly covariance model |
|
|
表 4 布格异常协方差模型特征参数统计表 Table 4 Statistics of characteristic parameters of Bouguer anomaly covariance model |
|
|
表 5 完全空间异常协方差模型特征参数统计表 Table 5 Statistics of characteristic parameters of Faye anomaly covariance model |
|
|
表 6 完全布格异常协方差模型特征参数统计表 Table 6 Statistics of characteristic parameters of complete Bouguer anomaly covariance model |
|
|
表 7 地形类别情况 Table 7 Terrain sorts |
上述统计表 3~7中，最小值为统计区域内重力异常的最小值，最大值为统计区域内重力异常的最大值，步长为协方差函数统计时出现负的最大距离，C0为方差，S为相关长度，B为协方差函数中的模型系数，m为利用统计计算拟合的协方差函数计算的协方差值与实际统计的协方差值之间的均方根误差.
对全国不同地形类别的完全空间异常和完全布格异常的方差C0，相关长度S，协方差函数中的系数B，利用统计计算拟合的协方差函数计算的协方差值与实际统计的协方差值之间的均方根误差m进行了统计.给出了它们的最大值（max）、最小值（min）及平均值（average），认为1°×1°内的同一地形类别的5′块大于100块，该1°×1°的地形类别为这一地形类别，统计情况见表 8、表 9，表中C0、S、B的单位与表 3~6中一致.
|
|
表 8 完全空间异常的C0、S、系数B分布 Table 8 Distribution of variance C0, correlation length S and model parameter B of Faye anomaly |
|
|
表 9 完全布格异常的C0、S、系数B分布 Table 9 Distribution of variance C0, correlation length S and model parameter B of Bouguer anomaly |
从统计表 8、9中可以看出，不同地区的相关长度S、方差C0、模型参数B也表现不同.由于它们均考虑了地形因素的影响，因此有利于我们进行拟合推估未知点重力异常，提高推值精度.以特大山区为例，利用表 9的参数统计值，其局部完全布格异常的协方差函数模型可写为：
|
(15) |
调和到外部空间：
|
(16) |
同理可得平原、丘陵、小山区等典型区域的完全空间异常与完全布格异常的协方差函数模型.但值得注意的是由于在中山区与大山区统计数据块少，这两类地形区域协方差还需进一步统计计算.
3.2 代表误差和代表误差系数根据代表误差的计算公式，本文的具体做法是：仅对具有4个以上重力点的5′×5′网格按公式（7）计算布格重力异常代表误差和空间重力异常代表误差；按地形类别计算5′×5′代表误差平均值；根据亨特经验公式计算代表误差系数，即
|
(17) |
式中E为计算的5′×5′代表误差平均值.依此即可分别得到5′×5′空间异常代表误差系数c和5′×5′布格异常代表误差系数cB，表 10是对它们计算结果的统计分析.
|
|
表 10 5'×5'空间重力异常和布格重力异常代表误差系数统计结果 Table 10 Statistics of commission error coefficients of 5'×5' free air anomaly and Bouguer anomaly |
从表 10中的数据可以看出，现行规定[11]中的重力异常代表误差系数与计算结果在平原、丘陵和小山区中的吻合情况较好，这既说明这些地区的重力点分布较好，也说明在该地区规定的代表误差系数是可靠的.而在中等、特大山区有一定的偏差，这说明随着地形起伏的增大，重力点在网格内的水平分布和高程分布都越来越不均匀.实际情况正好也是山区的重力点均匀性差，越是地形复杂的地区，山顶和山腰处测定的重力点越少.另外，表 10的数据还显示，新算的代表误差系数在不同地形类别地区的差别不大；规定值与本次计算值之差除有随高差增大而增大的情况外，差值的量级也普遍大于表 10.出现这种情况的原因除上述重力点分布不均匀外，还有布格重力异常在山区的变化比较平缓，与高程的相关性小，与地质密度的变化有较大关系；相同地形类别的不同地区，往往布格重力异常代表误差系数有较大的差异；此外，亨特经验公式用于小网格的计算往往使结果偏小也是原因之一.
采用前面空间异常与布格异常代表误差系数统计的计算方案，利用全国1′×1′平均高数据，统计了我国不同地形区域的完全空间异常、完全布格异常的代表误差系数，统计结果见表 11.表 12统计给出了不同分辨率地形类别的高程代表误差.此次统计工作基于我国不同典型地形区域大量重力与地形数据，具有一定的可靠性，统计结果将对我国重力测量合理布点、补点[17]提供指导和参考.
|
|
表 11 5'×5'完全空间重力异常和完全布格重力异常代表误差系数统计结果 Table 11 Statistical results of commission error coefficients of 5'×5' Faye anomaly and complete Bouguer anomaly |
|
|
表 12 1'×1'和5'×5'不同地形类别高程代表误差统计结果 Table 12 Statistical results of height commission error coefficients of 1'×1' and 5'×5' different terrain sort |
本文针对我国重力场以往统计分析工作的不足，构造了完全空间异常、完全布格异常的协方差和代表误差模型参数等新统计量，统计给出了平原、丘陵、小山区、中山区、大山区、特大山区6种不同地形类别区域的完全空间异常、完全布格异常的协方差和代表误差模型参数.试验表明顾及数据分布情况并优化软件实现方法后，依据本文的统计模型、算法与思路，在实际测量数据的支撑下，可以给出全国范围的统计参数的网格值、等值线图等.当全国高精度、高分辨率地形高程模型建立后[18]，结合地形信息，利用完全空间异常和完全布格异常进行拟合、推估空白区域的平均空间异常，以及其他扰动场元的最小二乘配置计算，将有利于精度的提高.
| [1] | Moritz H.Advanced Physical Geodesy.Abacus Press, England, 1980 |
| [2] | Heiskanen W A, Moritz H. Physical Geodesy. San Francisco: W.H.Freeman and Co., 1967 . |
| [3] | 孙风华. 陆地平均空间重力异常推值方法的研究. 解放军测绘学院学报 , 1997, 14(2): 97–102. Sun F H. Research of interpolation and extrapolation of land mean free air anomalies. Journal of Institute of PLA Surveying and Mapping (in Chinese) , 1997, 14(2): 97-102. |
| [4] | 方俊. 重力测量与地球形状学. 北京: 科学出版社, 1975 . Fang J. Gravimetry and the Earth's Shape (in Chinese). Beijing: Science Press, 1975 . |
| [5] | 管泽霖, 宁津生. 地球形状及外部重力场. 北京: 测绘出版社, 1981 . Guan Z L, Ning J S. The Earth's Shape and Outside Gravity Field (in Chinese). Beijing: Surveying and Mapping Press, 1981 . |
| [6] | Forsberg R A.Study of terrain reductions, density anomalies, and geophysical inversion methods in gravity field modeling.Reports of the Department of Geodetic Science and Surveying, Ohio State University, Report No.355, 1984 |
| [7] | Tscherining C C, Rapp R H.Closed covariance expressions for gravity anomalies, geoid undulations, and deflection of the vertical implied by anomaly degree variance models.Reports of the Department of Geodetic Science and Surveying, Ohio State University, Report No.208, 1974 |
| [8] | Cruz J Y.Disturbance vector in space from surface gravity anomalies using complementary models.Reports of the Department of Geodetic Science and Surveying, Ohio State University, Report No.366, 1985 |
| [9] | 朱长青. 计算方法及其在测绘中的应用. 北京: 测绘出版社, 1988 . Zhu C Q. Method of Calculation and Application in Surveying and Mapping (in Chinese). Beijing: Surveying and Mapping Press, 1988 . |
| [10] | 吴宗敏. 散乱数据拟合的模型、方法和理论. 北京: 科学出版社, 2008 . Wu Z M. Model, Method and Theory of Scattering Data Fitting (in Chinese). Beijing: Science Press, 2008 . |
| [11] | 陆仲连. 地球重力场的理论与方法. 北京: 解放军出版社, 1996 . Lu Z L. Theory and Method of the Earth's Gravity Field (in Chinese). Beijing: PLA Publishing House, 1996 . |
| [12] | 陆仲连, 吴晓平, 丁行斌, 等. 弹道导弹重力学. 北京: 八一出版社, 1992 . Lu Z L, Wu X P, Ding X B, et al. Theory of Trajectory Missile Gravity (in Chinese). Beijing: Ba-Yi Press, 1992 . |
| [13] | 黄谟涛, 翟国君, 管铮, 等. 海洋重力场测定及其应用. 北京: 测绘出版社, 2005 . Huang M T, Zhai G J, Guan Z, et al. Surveying of Ocean Gravity Field and Its Applications (in Chinese). Beijing: Surveying and Mapping Press, 2005 . |
| [14] | 黄维彬. 近代平差理论及其应用. 北京: 解放军出版社, 1992 . Huang W B. Theory and Application of Latter-day Adjustment (in Chinese). Beijing: PLA Publishing House, 1992 . |
| [15] | Lhaagmans R E, Erik D M, Gelderen M. Fast evaluation of convolution integrals on the sphere using 1D FFT, and a comparison with existing methods. Manuscripta Geodaetica , 1999, 18(5): 227-240. |
| [16] | Sideris M G. A fast Fourier transform method of computing terrain corrections. Manuscript Geodaetica , 1985, 10(1): 4-6. |
| [17] | 孙风华. 我国陆地均匀重力测量补点问题的研究. 武汉大学学报.信息科学版 , 2001, 26(4): 349–353. Sun F H. Research of added points problem of land gravimetric surveying in China. Geomatics and Information Science of Wuhan University (in Chinese) , 2001, 26(4): 349-353. |
| [18] | 胡鹏, 杨传勇, 吴艳兰, 等. 新数字高程模型. 北京: 测绘出版社, 2007 . Hu P, Yang C Y, Wu Y L, et al. New Digital Elevation Model (in Chinese). Beijing: Surveying and Mapping Press, 2007 . |