《1 引言》

1 引言

地球重力场是地球的一种物理场, 它反映了地球物质在时间和空间的分布状况。人类栖息在地球重力场之中, 一切活动均受重力场的影响。人类对地球重力场的了解不断深入, 重力场信息将作为“数字地球”的重要物理信息在大地测量、地球物理、空间技术、海洋学等诸多领域发挥重要作用。

地球重力场模型是描述和表示地球重力场的数学式及其相应参数的集合, 是对地球重力场和大地水准面形状的拟合逼近。用球谐函数的级数展开式表示的地球重力场模型由一组位系数予以体现, 是广泛采用的一种表示方法, 也称地球位系数模型。这种模型的优点之一是能够用位系数和相应的球谐展开式计算地球扰动引力场的其它场元, 如重力异常、垂线偏差、大地水准面高和空中扰动引力等。另一个优点是截断位系数模型到不同阶次能表示不同的空间分辨率, 便于在不同尺度上研究地球重力场和地球形状。

地球重力场模型的构制是一个不断发展的过程。20世纪初, 最简单的模型只有几个参数即所谓的椭球重力场模型, 在跨入21世纪时, 重力场模型参数已达数十万个, 空间分辨率大大提高。空间技术的发展, 需要用重力场模型精确计算对卫星轨道的重力摄动, 对卫星轨道的精确跟踪测量又为精化重力场模型提供了先进手段, 其中卫星海洋测高对高阶重力场的贡献尤为突出。在卫星上天之前, 是用地面重力数据解算位系数, 比较著名的是前苏联学者让戈洛维奇解算的8阶模型, 所用南半球的重力数据仅为北半球数据的百分之三, 用该模型计算的全球大地水准面最低处为-160 m, 最高处达到180 m, 与近期的数值相差近100 m。目前的重力场模型计算大地水准面的精度已达到1 m左右 [1]

重力场模型表示重力场和大地水准面的精度是很不均匀的。卫星跟踪站布设不足, 地面重力测量在两极和一些陆地覆盖不好, 都影响精度的均匀性。即使当前常用的360阶重力模型, 也只能在重力覆盖好的地区有较好的精度, 而在欠缺重力数据的地区, 用模型计算大地水准面的误差甚至达到9m之多 [2]

由于众所周知的原因, 最权威的国际机构也不可能掌握世界各地最详尽的重力数据。从20世纪80年代开始, 利用局部重力数据改进已知重力模型的实践和理论应运而生 [3,4,5]。我们在1984年构制的地球重力场模型DQM84首先采用了局部积分改进的方法, 当时是固定低阶位系数, 只改进高阶。1994年的DQM94模型构制时用了局部积分改进的谱权综合方法, 得到了更好的解算结果 [6]。近年来, 这种方法被称之为“剪接法”, 英文为“tailored approach”。

剪接法是相对于全球积分而言的, 因为最初人们怀疑, 位系数应采用全球重力异常积分解算, 局部积分改进会损害已知模型的精度。而我们改进的目的是, 使被改进的模型在描述局部重力场时, 比没有改进的要好, 而在其它地区, 改进后的模型应该与初始模型相当。局部改进总是必然的。正像全球重力异常, 都是先由局部改进, 然后才有全球的改进。

地球位系数又称地球质量分布系数, 低阶项的物理意义可以证明, 如零阶项与地球总质量有关, 一阶项与地球质量中心位置有关, 二阶一次项与地球旋转轴方向有关等。但位系数的阶数与异常质量的深度相关 [7], 所以位系数与重力异常一样, 都可以局部改进, 只是由于频域分析中要抑制局部波的传播。

基于以上思路, 我们在构制DQM99时突破了在同样阶次基础上进行局部改进的方法, 而是以360阶模型为初始模型, 利用15′×15′局部重力异常改进到720阶, 并利用5′×5′进行了局部积分改进到2 160阶的试验。对于初始模型为360阶的情况, 这种改进可称之为超阶改进。

《2 基本数学模型》

2 基本数学模型

地球引力位V的球谐级数表示式为

\(\begin{array}{l} V(r, \theta, \lambda)=\frac{G M}{r}\left[1+\sum_{n=2}^{\infty} \sum_{m=0}^{n}\left(\frac{a}{r}\right)^{n}\right. \\ \left.\left(C_{n m} \cos m \lambda+S_{n m} \sin m \lambda\right) P_{n m}(\cos \theta)\right] \end{array}\), (1)

 

式中, G为引力常数, M为地球总质量, r、θ、λ分别为计算点的地心向径、地心余纬和经度, a为参考椭球的长半径, Cnm、Snm是n阶m次完全正常化地球位系数, Pnm (cosθ) 是n阶m次完全正常化的缔合勒让德多项式。

参考椭球的正常位U的球谐展开式为

\(U(r, \theta)=\frac{G M}{r} \sum_{n=2}^{k}\left(\frac{a}{r}\right)^{n} C_{n 0}^{0} P_{n 0}(\cos \theta)\),  (2)

 

C0n0为参考椭球的n阶带谐位系数。

(1) 、 (2) 式中n从2开始, 是因为地球质量中心作为坐标系中心, 一阶项三个位系数为零。选择参考椭球质量与地球总质量相同, 旋转角速度相同, 则地球重力位与参考椭球正常重力位之差即为两者引力位之差, 称为扰动位, 以T表示, T=V-U, 球谐级数式为

\(\begin{array}{c} T(r, \theta, \lambda)=\frac{G M}{r} \sum_{n=2}^{\infty}\left(\frac{a}{r}\right)^{n} \sum_{m=0}^{n}\left(C_{n m}^{*} \circ\right. \\ \left.\quad \cos m \lambda+S_{n m} \sin m \lambda\right) P_{n m}(\cos \theta), \end{array}\),  (3)

 

C*nm表示相应位系数之差:

\(C_{n m}^{*}=C_{n m}-C_{n 0}^{0}\),  (4)

 

重力异常Δg、大地水准面高N、垂线偏差分量ζ、η等都是扰动位T的函数:

\(\begin{aligned} \Delta g &=-\left(\frac{\partial T}{\partial \cdot}+\frac{2}{r} T\right), \\ N &=\frac{T}{\gamma}, \\ \zeta &=\frac{1}{\gamma^{\circ} r} \frac{\partial T}{\partial \theta} \\ \eta &=\frac{1}{\gamma^{\circ} r \sin \theta} \frac{\partial T}{\partial \theta} \end{aligned}\) ,   (5)

 

上式中γ为地球平均重力。由位系数计算重力异常和大地水准面高的公式为

\(\begin{array}{c} \Delta g=\frac{G M}{r^{2}} \sum_{n=2}^{\infty}(n-1)\left(\frac{a}{r}\right)^{n+2} \sum_{m=0}^{n}\left(C_{n m}^{*}\right. \\ \left.\cos m \lambda+S_{n m} \sin m \lambda\right) P_{n m}(\cos \theta), \end{array}\)   (6)

\(\begin{aligned} N=& \frac{a^{2}}{r} \sum_{n=2}^{\infty}\left(\frac{a}{r}\right)^{n} \sum_{m=0}^{n}\left(C_{n m}^{*} \cos m \lambda+\right.\\ &\left.S_{n m} \sin m \lambda\right) P_{n m}(\cos \theta) \end{aligned}\),  (7)

 

由球谐函数的正交性, 已知全球的Δg或N, 都可以用球谐分析计算C*nm和Snm, 比较常用的是由重力异常计算高阶位系数, 而低阶位系数用卫星数据推求, 或者是卫星数据与地面重力数据的联合平差。由于高阶位系数个数较多, 一般用球谐分析数值积分计算。

《3 局部积分改进》

3 局部积分改进

参考文献[8]给出了局部积分谱权综合改进位系数的理论和方法, 这里仅给出构制DQM99的主要公式。

基于已知初始重力场模型的局部积分改进可以看成是位系数和局部重力异常的综合应用, 由于地球重力场扰动场元中重力异常Δg和大地水准面高N的频谱特性不同, 前者主要反映重力场的中、高频信息, 后者则反映中、低频信息。由位系数阶方差模型可以计算两者的谱信息量, 如表1所示。

  

《表1》

表1 不同频段Δg和N的能量分布

Table 1 Spectral power ofΔg and N for different

表1 不同频段Δg和N的能量分布Table 1 Spectral power ofΔg and N for different

设初始位系数为Cnm和Snm, 用局部重力数据计算的位系数改正数为ΔCnm和ΔSnm, 则改进后的位系数

\( \hat{\bar{C}}_{m m} 、 \hat{\bar{S}}_{n m}\)

\(\begin{array}{l} {\left[\begin{array}{l} \hat{\bar{C}}_{n m} \\ \bar{S}_{n m} \end{array}\right]=\left[\begin{array}{l} \bar{C}_{n m} \\ \bar{S}_{n m} \end{array}\right]+P_{n}\left[\begin{array}{l} \Delta \bar{C}_{n m} \\ \Delta \bar{S}_{n m} \end{array}\right],} \end{array}\)   (8)

 

式中Pn为权

\(\begin{aligned} P_{n} &=\frac{P_{n}(\Delta g)}{1+P_{n}(\Delta g)} \\ \end{aligned}\),   (9)

\(\begin{aligned} P_{n}(\Delta g) &=\left(\frac{r}{R}\right)^{2}(n-1)^{2} \frac{m_{N}^{2}}{m_{\Delta g}^{2}} \end{aligned}\),  (10)

 

(9) 式分母中的1实际上是位系数计算N的权, mN为初始位系数计算全球大地水准面高的中误差, mΔg为局部重力异常平均中误差。用局部重力数据计算位系数改正的公式为:

\(\left|\begin{array}{c} \Delta C_{n m} \\ \Delta S_{n m} \end{array}\right|=\frac{1}{4 \pi \gamma(n-1) S_{n m}(b / E)} \sum_{0} \frac{r}{a} \cdot \frac{1}{q_{n}} . \\ \overline{I P_{n m}(\cos \theta)} \sum_{\lambda}\left|\begin{array}{l} I C \\ I S \end{array}\right|_{m} \cdot\left[\overline{\Delta g}-\overline{\Delta_{g}}{ }^{M}\right],\)   (11)

 

式中

\(\gamma=\frac{G M}{a^{2}}\), Δg为局部区域重力异常, ΔgM为初始模型位系数计算的重力异常, Snm (b/E) 为第二类完全正常化勒让德多项式, qn是n阶平滑算子, IPnm为Pnm (cosθ) 在平均重力异常所在经纬格网中的积分平均值, IC和IS分别为cosmλ和sinmλ在格网中的积分平均值 [9,10]

 

《4 构模数据》

4 构模数据

《4.1 初始位系数模型》

4.1 初始位系数模型

初始位系数模型用了两个, 一个是美国俄亥俄大学的Rapp构制的360阶模型OSU91A, 另一个是美国Lemoine等构制的EGM96模型, 也是360阶。

OSU91A自90年代以来在国际上有广泛的应用, 是综合利用卫星和地面重力数据最好的模型之一。卫星数据是间接利用由31个卫星的1 100弧段推求的卫星解GEM-T2 50阶位系数。地面重力数据有1°×1°重力异常45 932个, 占全球总数的70.8%, 30′×30′重力异常66990个, 占全球总数的25.8%, 积分计算时, 将1°×1°重力异常一分为四, 取相同值为30′×30′重力异常。然后将所有30′×30′重力异常与GEM-T2 50阶位系数联合平差, 用平差后30′×30′重力异常数值积分计算51-360阶位系数, 与GEM-T2前50阶组合为2-360阶重力场模型 [10]

EGM96是1996年由美国国防部NIMA NASA和OSU等共同完成的。卫星解为70阶的JGM-3, 地面数据包括了前苏联、中国地区、格陵兰、南极洲、北极地区等重力数据, 但新增数据精度和覆盖分辨率差异较大, 如同OSU91A一样, 30′×30′重力异常与JGM-3 70阶位系数联合平差, 用平差后重力异常数值积分计算71-360阶位系数 [11]

《4.2 局部重力数据》

4.2 局部重力数据

我国及周边地区构成5′×5′重力异常200 692个, 在此基础上构成15′×15′重力异常22 866个, 30′×30′重力异常6 959个, 1°×1°重力异常1 790个, 重力系统均化为与初始模型相同的系统。

《4.3 误差估计》

4.3 误差估计

EGM96和OSU91A计算全球大地水准面的中误差mN按±0.5m估计, 在我国地区实际误差要大一倍。5′×5′重力异常中误差根据有关单位估计平均为±10×10-5m/s2, 30′×30′平均重力异常和1°×1°平均重力异常中误差分别为±7×10-5m/s2和±4×10-5m/s2, 这个估计可能也比实际精度高。

《5 DQM99解算及精度比较》

5 DQM99解算及精度比较

实际计算采用分频段逐级改进算法, 即先用1°×1°改进到180阶, 再用30′×30′从181改进到360, 361~720用15′×15′重力异常改进, 721~2 160阶用5′×5′重力异常改进。这样做是为了避免频谱的混叠, 便于实际使用时分频段利用位系数。解算结果见表2。

  

《表2》

表2 DQM99系列重力场模型

Table 2 DQM99 series gravity models

表2 DQM99系列重力场模型Table 2 DQM99 series gravity models

《5.1 模型重力异常与地面实际值比较》

5.1 模型重力异常与地面实际值比较

由于改进模型是对实际重力异常的拟合逼近, 所以比较的精度是内部符合精度。与初始模型共同比较, 可以了解改进后的模型是否更适合局部重力场, 在其它地区与地面重力异常比较, 则主要证明改进后的模型在表示其它地区重力场时应该与初始模型相等。表3列出各模型在不同地区与实际重力异常比较的数值结果。

  

《表3》

表3 模型重力异常拟合精度1×10-5m/s2

Table 3 Models gravity anomalies fitting accuracy

表3 模型重力异常拟合精度1×10-5m/s2Table 3 Models gravity anomalies fitting accuracy

由表3可以看出, 在局部地区, 改进后的重力模型较初始重力模型在重力异常拟合精度上提高一倍以上, 而在其它地区, 精度与初始模型相同。需要说明的是, 所用国外地区重力异常都是构制OSU91A时所用重力数据, 故OSU91A拟合精度优于EGM96, 这不表明OSU91A比EGM96精度高, 因为重力异常拟合是内部符合精度, 而模型大地水准面与GPS/水准比较是外部符合精度。

《5.2 模型大地水准面与GPS/水准比较》

5.2 模型大地水准面与GPS/水准比较

地球重力场模型应用最广泛的是计算大地水准面, 大地水准面表示了地球的基本形状。GPS技术可精确测定大地高, 当用几何水准联测出海拔高后, 即可得到该点大地水准面高。在以前, 大地水准面高必须由重力数据或位系数计算, 现在可以用GPS和水准测量出来, 而且非常精确, 在我国, GPS/水准确定的大地水准面的精度可达0.1m左右。故以GPS/水准点评价模型大地水准面精度是外部符合精度。我们比较的结果见表4、表5和表6。

从表4可以看出, 在外部符合精度比较中, OSU91A较EGM96有较大误差, 但是经过局部数据改进之后的DQM99A和DQM99B, 精度却是相当的, 这说明局部改进是有意义的。

  

《表4》

表4 与全国GPS一级网37点比较

Table 4 Model's geoid comparison with 37 GPS/m leveling points of the national wide GPS 1st order

表4 与全国GPS一级网37点比较Table 4 Models geoid comparison with 37 GPS/m leveling points of the national wide GPS 1st order

由表4还可以看出, 2 160阶模型DQM99C和DQM99D计算大地水准面的精度比720阶模型稍差, 其原因主要是因为位系数计算中对Legendre函数的递推计算在高阶函数中有一定的积累误差, 从而降低了高阶位系数的计算精度。

除了在全国范围的GPS一级网点上比较外, 还在华北地区GPS二级网的42点上进行了精度比较, 数值结果见表5。

  

《表5》

表5 与华北地区GPS二级网42点比较          m

Table 5 Model's geoid comparison with 42GPS/leveling points of the GPS 2nd order

表5 与华北地区GPS二级网42点比较m Table 5 Models geoid comparison with 42GPS/leveling points of the GPS 2nd order

为了证明改进后的模型在计算其它地区大地水准面时精度不受损失, 我们利用境外的高精度GPS/水准点作为标准对各模型大地水准面进行了精度比较, 结果见表6。

  

《表6》

表6 与其它地区模型大地水准面精度比较       m

Table 6 Models geoid comparison in other areas

表6 与其它地区模型大地水准面精度比较m Table 6 Models geoid comparison in other areas

《6 结论》

6 结论

全球重力数据的改进是由局部的改进到全球改进的漫长过程。地球重力场模型是对全球重力场的拟合逼近, 进行局部改进也是可能的。相同阶次的局部积分改进已有成功的实践, 而超出初始模型阶次的改进尚属首次。大量数据表明, 改进后的高分辨率地球重力场模型DQM99在表示局部重力场和大地水准面时精度有明显提高, 而对其它地区重力场和大地水准面与初始模型精度相同, 完全实现了进行局部改进的初衷。由于勒让德多项式是递推计算, 随着阶数增高因有效数位数有限, 产生误差积累, 所以2 160阶模型DQM99C和D的精度不如人意, 需要解决递推计算的精度。