卫星重力梯度数据用于精化地球重力场的研究

摘要

确定厘米级大地水准面和发展超高阶地球重力场模型是现代物理大地测量的主要科学目标之一。卫星重力梯度测量的实现将为这一目标做出重大贡献。文章着重评述这一领域的研究进展，并讨论利用卫星重力梯度数据精化地球重力场的若干理论和方法问题。

正文

《1 引言》

1 引言

地球重力场反映地球物质的空间分布、运动和变化, 制约着地球上及其邻近空间所发生的一切物理事件。确定地球重力场的精细结构及其时间相依变化不仅是物理大地测量的主要科学目标之一, 而且也将为现代地球科学解决人类面临的资源、环境和灾害等紧迫课题, 以及为有关大型工程 (如水利、桥梁工程等) 提供重要的基础地球空间信息。

现代大地测量、地球物理、地球动力学和海洋学等相关地学学科的发展迫切需要更加精细的地球重力场支持, 表1列出主要相关地学学科对重力场的要求。目前世界公认的最好的全球重力场模型EGM96与表1的要求相距甚远, 因此物理大地测量无论从重力探测技术还是地球重力场逼近理论和方法上均面临着新的挑战, 其中主要是确定厘米级大地水准面和发展超高阶全球重力场模型。

表1 相关地学学科对重力场的要求

Table 1 Requirements on the earth's gravity field from related earth science

《表1》

*1 mgal=10^-5 m/s²

地面重力测量、航空重力测量或航空重力梯度测量是在局部或区域范围内获取高精度重力场信息的有效技术手段, 适用于恢复短波重力场。卫星测高目前是精化海洋大地水准面最有效的技术手段, 但由于受海潮残差、物理环境影响的残差和定轨误差等多种难于做精确模拟的误差影响, 利用该技术恢复海洋大地水准面优于±10 cm量级的准确度十分困难。由地面跟踪观测卫星轨道摄动适于确定长波重力场 (2～36阶) , 但轨道摄动不仅受重力异常场的影响, 而且还受日月引力、地球潮汐、大气和太阳光压等物理环境的影响, 其残余误差目前仍大于厘米级;其次, 由于卫星轨道的全球覆盖不全, 若干高阶位系数之间强相关, 会产生混叠和共振影响等, 只有精心设计多种不同卫星轨道倾角和运行周期才能解决, 显然要发射大量有特定轨道根数的卫星是不现实的, 这种影响也大于厘米级水平。因此, 现有重力探测技术获取全球均匀分布的高精度重力场信息的能力受到了限制, 迫切需要新的技术突破。在这一背景下出现了发展卫星重力梯度测量 (SGG) 和卫星跟踪卫星 (SST) 技术的设想, 其中SGG是利用卫星携带的重力梯度仪直接测定引力位的二阶导数张量来确定地球重力场。下面着重讨论卫星重力梯度测量技术的进展及其在精化地球重力场中的若干理论和方法问题。

《2 卫星重力梯度测量技术的进展》

2 卫星重力梯度测量技术的进展

自20世纪初匈牙利物理学家R.E⌀tv⌀s设计出第一台重力梯度仪 (E⌀tv⌀s扭秤) 以来, 重力梯度仪大致经历了从单轴旋转到三轴定向, 从室温到低温 (低于4.2 K) , 从扭力、静电悬浮到超导的发展过程, 仪器精度日益提高。特别是随着现代电子技术、计算机技术、超导量子干涉技术、低温微波空腔谐振技术、“超导负弹簧”技术的发展及应用, 使得重力梯度仪在灵敏度和稳定性方面取得了突破性进展。重力梯度仪的精度从早期的1E已发展到10^-4E, 未来超导重力梯度仪的精度可望达到10^-6E ^[1]。重力梯度测量原理从早期的扭力测量发展到目前的差分加速度测量, 前者测定作用于检测质量的力矩来间接获取重力梯度数据, 后者通过观测两检测质量 (或加速度计) 之间的加速度差来获取重力梯度观测值。目前卫星重力梯度测量主要采用差分加速度测量技术。表2列出了目前重力梯度仪的主要研制机构、类型及精度。正是重力梯度测量技术的突破, 推动了卫星重力梯度测量的迅速发展。

表2 重力梯度仪的类型及精度

Table 2 Types and precision of the gravity gradiometer

《表2》

国外自20世纪70年代末就开展SGG技术的论证, 对其测量原理、技术模式、误差源、数据处理的理论和方法等进行了大量的数值模拟实验和分析。80年代开始研究制定国际卫星重力梯度计划, 其中有欧洲空间局 (ESA) 的ARISTOTELES计划 (application and research involving space techniques observing the earth field from low earth orbit satellite) 、美国宇航局 (NASA) 的超导重力梯度测量任务SGGM (superconducting gravity gradiometer mission) 和美意两国的系留卫星系统TSS (tethered satellite system) 的重力梯度测量等。经过近30多年的潜心研究, 卫星重力梯度测量技术已趋向成熟。欧洲空间局和美国宇航局今后几年将陆续发射具有测定地球重力场能力的卫星, 如CHAMP (challenging minisatellite payload for geoscience and application) 、GRACE (gravity recovery and climate experiment) 和GOCE (gravity field and steady-state ocean circulation explorer mission) , 其中CHAMP是用于地球物理研究的小卫星, 采用高-低卫星跟踪卫星技术模式, 其观测量可认为是低卫星处的重力位一阶导数。该卫星已于2000年7月15日 (UTC时间) 在俄罗斯的Plesetsk卫星发射基地成功发射, 目前运行状况良好。GRACE是“探测重力场和气象实验”的卫星探测计划, 计划于2001年11月实施, 同时采用高-低和低-低SST技术, 其中低-低SST实际上相当于一个重力梯度仪, 其观测量为两卫星视线方向上的重力梯度分量, 可用于恢复150阶地球重力场, 大地水准面的精度可达到20 cm。GOCE是所谓“重力场和静态洋流探索"的卫星探测计划, 同时实施高-低SST技术和卫星重力梯度测量技术, 预计于2005年发射, 期望恢复250阶地球重力场和1 cm精度的大地水准面。SST和SGG被认为是目前最有价值和最有应用前景的两项互为补充的高效重力探测技术, 其主要科学目标不仅测定地球重力场的精细结构及长波重力场随时间的变化, 而且还将以全球尺度精密测定磁场及探测全球大气层和电离层。这一系列卫星重力探测计划的实施将为重力场的研究树立新的里程碑, 同时对现代地球科学研究岩石圈、水圈和大气圈及其相互作用具有重大科学和现实意义。

《3 卫星重力梯度边值问题》

3 卫星重力梯度边值问题

卫星重力梯度测量有不同的技术模式, 它们可利用卫星上的重力梯度仪直接测定引力梯度张量的所有分量或部分分量或某些分量的组合。在某一给定的坐标系下, 引力梯度张量可表示为

$Γ = [\begin{matrix} V_{x x} & V_{x y} & V_{x z} \\ V_{y x} & V_{y y} & V_{y z} \\ V_{z x} & V_{z y} & V_{z z} \end{matrix}] 。 (1)$ $Γ = [\begin{matrix} V_{x x} & V_{x y} & V_{x z} \\ V_{y x} & V_{y y} & V_{y z} \\ V_{z x} & V_{z y} & V_{z z} \end{matrix}] 。 (1)$

根据重力场的性质, 该张量是一个对称无迹的偏差张量, 只有五个独立分量。利用这一新型重力场数据确定地球重力场, 理论上应归结为求解卫星重力梯度边值问题。

按照物理大地测量的传统作法, 选择一合适的参考重力场, 则重力场的研究归结为对扰动场的研究。因卫星的轨道面非常接近于一个球面, 因此选择以卫星标称轨道R_S为半径、地心为球心的球面S作为边界面, 同时假设卫星轨道面上的重力梯度数据已归算到该球面上并覆盖整个边界面, 则卫星重力梯度边值问题可表述如下:

${\begin{cases} Δ Τ = 0, ? 在 S 外； \\ B_{i} Τ = f_{i}, 在 S 上 ‚ i = 1, 2, \dots, 6 ； \\ Τ = 0 (r^{- 3}), 在无穷远处。 \end{cases} (2)$ ${\begin{cases} Δ Τ = 0, ? 在 S 外； \\ B_{i} Τ = f_{i}, 在 S 上 ‚ i = 1, 2, \dots, 6 ； \\ Τ = 0 (r^{- 3}), 在无穷远处。 \end{cases} (2)$

这里T为扰动位, B_i为二阶边界微分算子, f_i为扰动重力梯度数据 (边值) 。式 (2) 也可简记为

${Δ; B} Τ = {0; f}, Τ = 0 (r^{- 3}) ‚ (3)$ ${Δ; B} Τ = {0; f}, Τ = 0 (r^{- 3}) ‚ (3)$

其中B为边界微分算子组, f为边值组。

在式 (2) 或式 (3) 中, 如果只考虑某一类重力梯度数据 (如径向二阶重力梯度) , 则该边值问题可称为卫星重力梯度单定边值问题 ^[2,3,4];如果同时考虑多类重力梯度数据, 则应归结为求解卫星重力梯度超定边值问题 ^[5,6,7,8,9]。下面着重讨论卫星重力梯度超定边值问题的准解模型。

考虑Sobolev空间H² (S) 的迹定理B_iT∈H^p_i (S) , 其中 $p_{i} = \frac{1}{2}, i = 1, 2, \dots ‚ 6$ $p_{i} = \frac{1}{2}, i = 1, 2, \dots ‚ 6$ , 定义乘积空间为

$\begin{array}{l} J \equiv Η^{p_{1}} (S) \times Η^{p_{2}} (S) \times Η^{p_{3}} (S) \times \\ Η^{p_{4}} (S) \times Η^{p_{5}} (S) \times Η^{p_{6}} (S) \\ = {f = (f_{1}, f_{2}, f_{3}, f_{4}, f_{5}, f_{6}); \\ f_{i} \in Η^{p_{i}} (S), i = 1, 2, \dots ‚ 6} 。 (4) \end{array}$ $\begin{array}{l} J \equiv Η^{p_{1}} (S) \times Η^{p_{2}} (S) \times Η^{p_{3}} (S) \times \\ Η^{p_{4}} (S) \times Η^{p_{5}} (S) \times Η^{p_{6}} (S) \\ = {f = (f_{1}, f_{2}, f_{3}, f_{4}, f_{5}, f_{6}); \\ f_{i} \in Η^{p_{i}} (S), i = 1, 2, \dots ‚ 6} 。 (4) \end{array}$

再定义J的闭子空间为

$\tilde{J} \equiv {B Τ; Τ \in Η^{2} (S) ‚ Δ Τ = 0, Τ = 0 (r^{- 3})} 。 (5)$ $\tilde{J} \equiv {B Τ; Τ \in Η^{2} (S) ‚ Δ Τ = 0, Τ = 0 (r^{- 3})} 。 (5)$

所谓卫星重力梯度边值问题的准解, 实际上是求f∈J在闭子空间 $\tilde{J}$ $\tilde{J}$ 中的投影。省略复杂的推导过程, 直接给出准解的简化模型为

$A X = L, (6)$ $A X = L, (6)$

其中A是阶数为[ (N_max-n) /2+1]、带宽为5的带状方阵, 其元素是关于重力场模型的阶次n和m的代数表达式, N_max是重力场模型的最大阶数;L是以边值f_i的球谐展开系数 (由球谐分析得到) 的代数表达式, 为元素的列向量;X是由扰动位的归格化位系数 $\bar{a}$ $\bar{a}$ _nm或 $\bar{b}$ $\bar{b}$ _nm组成的待求未知列向量。A、X和L的具体表达式详见文献[4]。因此, 卫星重力梯度边值问题的准解归结为求解式 (6) 的线性带状方程组。模拟试算表明, 该模型是有效的, 并便于计算。需要指出, 上述结果是在球近似下导出的, 没有考虑非线性影响和观测数据的噪声问题。模型本身有待进一步完善, 使之更加符合卫星重力梯度测量的实际情况。

对于利用卫星重力梯度数据恢复地球重力场位系数模型, 还可采用球谐分析 ^[10,11,12]或广义球谐分析 ^[13]方法计算, 但这些方法只能处理某一种重力梯度分量 (如V_zz和V_yz) 或某些重力梯度分量的组合形式, 不能充分利用重力梯度信息。此外, 还有一种球近似下的最小二乘迭代解法 ^[14], 采用重力梯度分量的某些线性组合较好地分离了卫星轨道摄动。模拟研究表明, 该方法通过迭代计算能获得较准确的卫星轨道和地球重力场模型。由于这种方法将重力梯度观测值看成是卫星轨道位置的函数, 因此习惯上称之为空域最小二乘法。沿卫星轨道采集重力梯度数据实际上是一个离散的时间采样过程, 重力梯度观测值可看成是时间的函数, 因此基于卫星轨道摄动分析理论, 空域最小二乘法又发展到时域 (称为时域最小二乘法) ^[15,16]。上述这些方法在理论上各有优缺点, (广义) 球谐分析法计算简单实用, 但由于重力梯度数据实际上不可能连续分布, 因而存在积分公式离散化的误差, 不能估计重力梯度数据误差对求解位系数的影响和同时处理所有重力梯度分量。空域最小二乘法不要求数据连续分布, 因而容易克服两极无数据的缺陷和解决卫星重力梯度测量过程中的数据间断问题, 其缺点是难以反映卫星的轨道特征和顾及有色噪声的影响。时域最小二乘法因基于卫星轨道摄动分析理论, 因而能较准确地反映卫星的轨道特征并顾及卫星的轨道误差、姿态误差及某些有色噪声的影响, 适合于卫星重力梯度测量的误差分析和方案设计, 但该方法难以处理两极地区无重力梯度数据和因卫星操作等原因造成的数据间断问题。

利用覆盖全球的卫星重力梯度数据解上述边值问题可以确定全球重力场, 例如完全阶次至250阶的全球重力场模型, 并能大大改善其中短波分量。

《4 卫星重力梯度的局部解》

4 卫星重力梯度的局部解

所谓卫星重力梯度的局部解, 实际上是将卫星轨道上的重力梯度观测值向下延拓到地球表面或大地水准面上, 以确定扰动位泛函 (如重力异常, 大地水准面差距等) , 从而达到精化局部或区域重力场的目的。目前在解算方法上可归纳为最小二乘法、最小二乘配置法和组合型方法 ^{[17,18,19,20]}。其中最小二乘法的关键问题是不能顾及信号场相关统计信息的影响, 而最小二乘配置法却难于构造结构良好的协方差阵;最佳积分核谱组合解在理论上具有稳定性和最小方差的特征, 这类似于最小二乘配置法的优点, 但实用上存在如何合理确定谱权和计算量过大的问题。针对以上几种方法存在的问题, 笔者曾提出一种卫星重力梯度向下延拓的谱组合法 ^[21]。

在平面近似时, 利用卫星重力梯度分量确定大地水准面或地面上扰动位T (x, y, 0) 的谱表达式为

$\begin{array}{l} F [Τ (x, y, 0)] = F [Τ_{i j} (x, y, h)] \cdot \\ G (u, v)_{i j} ‚ i, j = x, y, z (7) \end{array}$ $\begin{array}{l} F [Τ (x, y, 0)] = F [Τ_{i j} (x, y, h)] \cdot \\ G (u, v)_{i j} ‚ i, j = x, y, z (7) \end{array}$

计算扰动重力分量T_x (x, y, 0) 、T_y (x, y, 0) 和T_z (x, y, 0) 的谱表达式为

$\begin{array}{l} F [Τ_{i} (x, y, 0)] = F [Τ_{i j} (x, y, h)] \cdot \\ G_{i} (u, v)_{i j}, i, j = x, y, z (8) \end{array}$ $\begin{array}{l} F [Τ_{i} (x, y, 0)] = F [Τ_{i j} (x, y, h)] \cdot \\ G_{i} (u, v)_{i j}, i, j = x, y, z (8) \end{array}$

这里F为二维Fourier变换算子, T_ij为扰动重力梯度分量, u、v为对应空域直角坐标变量x和y的频域变量, G (u, v) _ij和G_i (u, v) _ij的表达式见表3 (表3T式中 $i = \sqrt{- 1} ‚ q = \sqrt{u^{2} + v^{2}})$ $i = \sqrt{- 1} ‚ q = \sqrt{u^{2} + v^{2}})$ 。

表3 谱函数G (u, v) _ij, G_x (u, v) _ij, G_y (u, v) _ij和G_z (u, z) _ij的表达式Table 3 Representative of spectral function G (u, v) _ij, G_x (u, v) _ij, G_y (u, v) _ijand G_z (u, v) _ij

《表3》

假设由单一重力梯度分量按式 (7) 计算扰动位, 其谱用T^j (u, v) 表示, 功率谱密度函数的估值为 $\hat{S}$ $\hat{S}$ _T^j (u, v) (可采用最大熵谱分析或经典谱估计方法求得) , 则由多个重力梯度分量计算F[T (x, y, 0) ]的最佳估值为

$\begin{array}{l} F [\hat{Τ} (x, y, 0)] = [\sum_{j = 1}^{Κ} \hat{S}_{Τ^{j}} (u, v) Τ^{j} (u, v)] / \\ \sum_{j = 1}^{Κ} \hat{S}_{Τ^{j}} (u, v), (9) \end{array}$ $\begin{array}{l} F [\hat{Τ} (x, y, 0)] = [\sum_{j = 1}^{Κ} \hat{S}_{Τ^{j}} (u, v) Τ^{j} (u, v)] / \\ \sum_{j = 1}^{Κ} \hat{S}_{Τ^{j}} (u, v), (9) \end{array}$

其中K为参加计算的重力梯度分量的个数。

对式 (9) 做Fourier逆变换, 即可得到大地水准面或地面上扰动位的最佳估值为

$\hat{Τ} (x, y, 0) = F^{- 1} {F [\hat{Τ} (x, y, 0)]} 。 (10)$ $\hat{Τ} (x, y, 0) = F^{- 1} {F [\hat{Τ} (x, y, 0)]} 。 (10)$

根据卫星重力梯度数据计算扰动重力分量T_x (x, y, 0) 、T_y (x, y, 0) 和T_z (x, y, 0) 的方法与计算扰动位T (x, y, 0) 类似。上述方法的最大优点是利用FFT (或FHT) 技术将平面积分和向下延拓一步完成, 既具有高效快速的特点, 又可分析局部重力场的谱特性。其主要缺陷:因采用球面的平面近似, 重力梯度数据序列的线性卷积与Fourier变换二维循环卷积存在差异, 数据采样不充分和数据长度有限将导致混叠和频谱泄漏现象。此外, 从表3可以看出, 指数运算算子在实际计算中起放大作用, 如果重力梯度数据存在小的误差, 则有可能对计算结果产生较大的歪曲, 这是由卫星重力梯度向下延拓的不适定性造成的, 在实际计算中必须慎重考虑。无论采用上述何种方法, 在重力梯度观测值中消去重力场长波成分和利用地形模型平滑其高频成分, 将有助于提高恢复局部或区域重力场的精度和探测原始数据中的粗差。若要联合卫星重力梯度数据及其他重力场信息确定更高精度的局部重力场模型, 则需要进一步发展高效实用的新方法。

《5 若干建议》

5 若干建议

物理大地测量的理论与实践表明, 不同的重力探测技术只能获取相应频谱的重力场信息, 只有尽可能多地联合利用具有不同频谱特性的重力场信息, 才能达到精化全频段地球重力场的目的。尽管过去已取得大量的研究成果, 但在利用卫星重力梯度数据精化地球重力场方面, 笔者认为还需在以下领域进行深入研究。

1) 卫星重力梯度超定边值问题

如果以米级或分米级精度逼近地球重力场, 则上述关于卫星重力梯度边值问题的理论和解算方法是合理的。但如果以厘米级精度逼近地球重力场, 则必须深入研究卫星重力梯度超定边值问题中边界面的定义、重力梯度数据的归算方法、非线性误差的影响及解的适定性等。由于地壳构造的复杂性和地形表面的粗糙性, 重力位二阶导数在这一边界面上显示强非线性和无序性结构, 经典的数理方程求解所涉及的数学工具已不能严密处理这一边界面, 必须着手深入研究非线性和无序性边值问题的基础理论, 寻求新的适用的数学工具, 并将卫星重力梯度超定边值问题的理论和方法向非线性领域拓展, 使之更贴切于地球重力场实际。

2) 地球重力场的综合确定

卫星重力梯度测量不仅提供重力梯度数据, 而且还可获得SST数据 (如GOCE) , 联合这些数据及其他重力场信息 (如卫星测高数据、地面重力数据和航空重力数据等) 综合确定地球重力场, 除了继续发展上述超定边值问题和非线性边值问题基础理论以外, 在统计法的研究中应发展“非经典配置法”, 即利用现代数学的新理论和新成果 (如模糊论、混沌论和分形论等) , 来研究逼近地球重力场精细结构的新方法。近年来, 小波分析理论、多输入输出法 (multi-input output method) 及输入/输出系统理论算法 (input/output system theory method) 等已在地球重力场精化中初步展示了其应用前景, 特别是小波分析法可望在地球重力场的联合求解中发挥重要作用。目前问题的关键是, 如何根据各类重力场信息的频谱特性建立高效实用的小波分析算法模型。

3) 随机边值问题

在联合多种数据逼近地球重力场精细结构方面, 除了在上述解析法和统计法两个领域里继续开拓新思路以外, 还应在两者的结合上对地球重力场地面边界的无序结构做深入探索。随机边值问题是目前有可能将求解大地测量边值问题的解析法和统计法融为一体的新的研究方向。以厘米级精度水平逼近地球重力场为目标, 研究该问题所涉及的主要内容应包括:如何根据卫星重力梯度数据及其他重力场信息的实际分布定义合适的边界面和合理的数据归算方法;随机边值问题的严密数学表述和具体求解原则;论证该问题解的适定性, 给出适定解的限制性条件和建立实用的解算模型等。随机边值问题的研究及应用, 将为地球重力场的精化提供新的途径。

展示更多