地球重力场的空间和时间变化以不同的方式影响卫星的轨道。以GRACE卫星为例,当GRACE卫星围绕地球旋转时,这种重力变化以GRACE-A和GRACE-B两卫星间距离变化的形式表现出来,这些距离变化与其它的附加观测数据可用来分析和提取地球重力场时空变化模型参数。当坐标原点与地球质心重合时,地球重力位可通过球谐系数表示为

V(r,φ,λ)=(GM)/r{1+∑lmaxl=2∑lm=0(a/r)lP^-lm(sinφ)

[C^-_lmcosmλ+S^-_lmsinmλ]}.(1)

式中,GM为地球引力常数(单位m³·s^-2); r是外部一点到地心的距离; φ、λ分别是地心纬度和地心经度; a是地球平均赤道半径; l_max是位系数能够恢复的最大阶数; P^-_lm(sinφ)为l阶m次的规格化勒让德函数; C^-_lm、S^-_lm为规格化球谐系数。

地球对外部空间的引力g表示为

g(r,φ,λ)=-(GM)/(r2){1+(l+1)∑lmaxl=2∑lm=0(a/r)l

P^-_lm(sinφ)[C^-_lmcosmλ+S^-_lmsinmλ]}.(2)

将GRACE卫星获得的重力场与背景重力场求差,得到重力场的变化值,即

Δg=g-g_背景.(3)

相应的重力变化值为

Δg=-(GM)/(r2)∑lmaxl=2(l+1)(R/r)l∑lm=0WlP^-lm(sinφ)[ΔC^-lmcosmλ+ΔS^-lmsinmλ].(4)

其中,ΔC^-、ΔS^-分别为重力场与背景重力场模型的球谐系数之差; W₁为高斯平均函数。理论上地球重力场模型只有在其包含从0到无穷的所有阶次的球谐系数时,才能反映全波段的重力场信息。而实际上只能得到某一最大阶数,GRACE卫星恢复的重力场只能恢复某一固定阶次,因此只能对部分频谱进行分析,而且GRACE恢复的重力场模型的阶数越大,其误差也越大,加入算子后,(4)式高阶系数的权重减小,高阶系数误差的影响减弱,截断的误差影响减小,计算结果更可靠(段虎荣等,2009)。

在实际研究地质体所在的空间时,一般将区域分成若干个大小形状相同的规则直立长方体,直立长方体各边分别平行于三个坐标轴,可以认为每个直立长方体由许多个点源组成,在单一的长方体中密度是均匀的,用积分的方法可以得到直立长方体产生的重力。

现介绍直立长方体任意一点重力的计算。取1个直立长方体并置于图1所建立坐标系中,假设Q(ε,η,τ)为直立长方体内任一点坐标,其直立长方体的空间取值范围为:ε∈[a₁,a₂], η∈[b₁,b₂],

图1 直立长方体外任意一点的重力
Fig.1 Gravity of any point outside the erect cube

V=(Gm)/r=Gρk∫ a2a1∫ b2b1∫ h2h1

1/(((ε-x)²+(η-y)²+(τ-z)²)^1/2)dεdηdτ.(5)

其中,r=((ε-x)²+(η-y)²+(τ-z)²)^1/2; G为万有引力常数; m为直立长方体的质量。

通过对上式求导g(x,y,z)=(V)/(Z),进行积分后的的点P(x,y,z)的重力为

g(x,y,z)=-Gρ〖JB<2*|〗〖JB<2|〗〖JB<1*|〗(ε-x)·ln[(η-y)+r]+(η-y)·ln[(ε-x)+r]+(τ-z)arctan[((τ-z)·r)/((η-y)(ε-x))]〖JB>1*|〗^a2_a1〖JB>2|〗^b2_b1〖JB>2*|〗^h2_h1.(6)

由(6)式假定ρ为常数,那么点P(x,y,z)的重力就是直立长方体大小的函数,可以表示为

g(x,y,z)=f(a₁,a₂,b₁,b₂,h₁,h₂).(7)

当直立长方体发生运动(dx,dy,dz)时,那么重力变化可表示为

Δg(x,y,z)=f(a₁+dx,a₂+dx,b₁+dy,b₂+dy,h₁+dz,h₂+dz)-f(a₁,a₂,b₁,b₂,h₁,h₂).(8)

当研究区域被划分成n个直立长方体,因为重力的可加性,由所有直立长方体运动在p(x,y,z)点产生的重力变化为

Δg^→(x,y,z)=∑ⁿ_k=1Δg^→_k=F(dx,dy,dz).(9)

美国普渡大学的Kennedy和Eberhart(1995)提出了PSO算法,该算法具有很强的全局优化能力。为了证明计算地壳水平运动引起重力变化的正确性以及PSO算法的有效性,笔者在图4中将模拟的研究区域划分成5×5个直立长方体(100 km×100 km),假设其密度均匀,每个直立长方体Z₁=0.001 m、Z₂=30 000 m,对每个直立长方体赋予水平运动速率,如图4中黑色箭头所示,箭头方向为直立长方体的水平运动方向,大小为运动速率。则由于水平运动引起的重力场变化(计算点为每个直立长方体在Z=0 的平面上的投影正方形的中心坐标(x,y,0)。由图4可以看出重力变化的数值大小与水平运动速率相关,在物质水平运动方向上,物质流入使得空间位置质量增大,导致观测点重力增加,相反,物质流出使得空间位置质量减少,导致重力减小。整体来看重力变化的梯度与水平运动方向相一致。

将重力变化、直立长方体的大小、位置、密度等作为已知量,根据公式(8)采用PSO算法反演直立长方体的水平运动量参数,在算法中取粒子个数为100,粒子维数为100,最大速度V_max=0.2。图4中黑色箭头为水平运动速率的理论值,蓝色为反演结果,总体结果与理论吻合,个别点的大小有偏离,这与目标适应度等因素有关。

图4 直立长方体的运动速率及其产生的重力变化
Fig.4 Horizontal movement velocity of erect cube and it cause gravity changes

笔者将中国西部地区(21°～45°N,75°～105°E)划分成7行10列(共70个直立长方体),计算了2008年12个月的重力场模型平均值相对于2003年平均值的变化,模型的最大阶为120,高斯平滑半径为390 km。设目标函数为

minΦ(P)=∑^N_i=1((Δg^a_i-Δg^b_i)/(Δg^a_i))².(11)

式中,Δg^a为卫星重力模型计算的重力变化值; Δg^b为地壳水平运动引起的重力变化; N为划分的直立长方体个数。然后根据目标函数评价适应度函数

Q_t=(1.0)/(Φ(P)_t+1.0).(12)

本文采用PSO算法反演地壳水平运动速率时,取粒子的维数为140,取粒子的最大速度取V_max=0.2,加速度因子取c₁=c₂=0.5。当适应度为0.95时结束程序运算,反演结果见图5。

图5 中国西部地壳水平运动速度场的反演结果
Fig.5 Inversion of crustal horizontal movement velocity field in West China

表1 卫星重力反演的地壳速率场与GPS观测结果对比
Tab.1 Comparision of crustal horizontal movement velocity by satellite gravity and GPS observation

从图5可看出:中国西部地壳运动速率的反演结果与GPS观测结果大体相符,整体具有运动速率不均匀的特点。地壳运动由南向北逐渐递减,且存在向东运动的趋势,在昆明区域,地壳由北向南运动,青藏高原南北方向运动速率为10～35 mm/a,东西方向运动速率1～25 mm/a; 塔里木盆地南北方向运动速率为5～15 mm/a,东西方向运动速率为4～7 mm/a。为了进一步验证由卫星重力计算的地壳运动速率,我们对部分地区的地壳运动速率与GPS观测结果(王琪等,2002)进行对比(表1),除了1号直立方体与GPS观测结果差值较大外,其余的直立方体的结果与GPS观测结果的趋势一致,只是计算结果总体偏大。