如何精确确定地球重力场模型是当前大地测量学及其相关学科研究的重点和热点问题[1-2]．高精度的重力场模型对地球动力学、海洋学、地球物理学和地震学等相关学科研究具有重要的作用[2-5]．目前，超高阶重力场模型的建立主要依靠解算某类边值问题，如斯托克斯边值问题、莫洛金斯基边值问题[5-6]．然而这些边值问题所要求的边界面是一个不规则的曲面，如斯托克斯边值问题要求的边界面为大地水准面；莫洛金斯基边值问题要求的边界面为似地球表面．直接处理这样的不规则曲面，从数学上是难以实现的，因此通常采用近似的方法简化这些不规则曲面．早在20世纪五六十年代，受制于数据观测精度低，观测区域有限等因素，处理这些边值问题时通常将边界面近似成平均球面．随着观测技术的进步和数据观测精度的提升，已有的球近似边值问题不再满足要求，于是建立了椭球边值问题．对Dirichlet，Neumann及混合类型的椭球面边值问题已有众多学者进行了研究[6-12]，可见椭球边值问题相较于球边值问题而言，当构建重力场模型时精度可以提高椭球扁率平方量级，即O(ε4)．另外，当构建超高阶重力场模型时，所采用的斯托克斯边值问题或莫洛金斯基边值问题都是线性形式，舍去了非线性项，即O(T2)[10,12]．超高阶地球重力场模型的建立是对全球重力观测数据进行调和分析或最小二乘平差的结果．目前，构建超高阶地球重力场模型的数据主要包括地面重力、海洋重力、航空重力和卫星重力等．对每一种观测数据而言，受观测技术、观测方法等多种因素的影响，每种类型数据的精度并不一致．因此地球重力场模型的精度除了理论模型精度外，还取决于输入模型的观测数据的精度[2,5]．为了获得高精度、高分辨率的超高阶地球重力场模型，国内外机构投入了巨大的人力、物力和财力进行了相关研究．目前重力场模型已经发展到 2 159完全阶次[13]，例如目前应用最广的EGM2008模型、EIGEN系列模型等．虽然超高阶地球重力场模型达到2 159阶次，但依然存在一个问题，即不同数据观测精度对超高阶重力场模型的影响尚不清楚．为此首先须要验证线性化的边值问题是否满足超高阶重力场模型构建的要求；再进一步分析不同观测数据精度对构建超高阶重力场模型阶数的影响．1 线性化的斯托克斯边值问题根据定义，大地水准面上的重力异常可写为[6]Δg=gp-γq，(1)式中：gp为大地水准面上的重力；γq为参考椭球上的正常重力；p和q分别为其对应的点；γ为大地水准面上的正常重力．那么，在p点的正常重力沿着q点进行泰勒展开，即γp=γq+∂γ∂hN+12!∂2γ∂h2N2+⋯，(2)式中：N为大地水准面高；h为大地水准面法线方向．由布隆斯公式可知Tγ=N+12!1γ∂γ∂hN2+⋯，(3)式中T为扰动位．由于大地水准面高在±100  m的范围内，那么N/R的量级小于椭球扁率的量级[14]，因此式(3)只须考虑线性项，即T/γ=N．(4)将式(2)和(4)代入式(1)中有Δg=gp-γp+1γ∂γ∂hT+12!1γ2∂2γ∂h2T2+⋯，(5)根据gp-γp=-∂T/∂h，式(5)可写成Δg=-∂T∂h+1γ∂γ∂hT+12!1γ2∂2γ∂h2T2+⋯．(6)忽略掉非线性项O(T2)，线性化的斯托克斯边值问题可写成ΔT=0    (在Ω外);-∂T∂h+1γ∂γ∂hT|Ω=Δg    (在Ω上);T=O(r-1)    (在无穷远处), (7)式中Ω为大地水准面．事实上对于ΔT=0，斯托克斯边值问题(7)的解可写成扰动位T的球谐级数的形式，即T=∑n=2∞∑m=-nn(R/r)n+1c¯nmWnm(θ,λ)，(8)式中：R为地球平均半径；(r,θ,λ)为球坐标系下的分量；r为原点矩；θ为余纬；λ为经度；c¯nm为球谐系数，具体为c¯nm=a¯nm(m≥0),b¯nm(m0);Wnm(θ,λ)=P¯nmcosθcos(mλ)(m≥0),P¯nmcosθsin(mλ)(m0),其中P¯nm(cosθ)为归一化的Legendre函数．为了叙述方便，对于边值问题(7)中的边界面为球面近似的情况则称为斯托克斯球边值问题，对于边界面为椭球面近似的情况则称为斯托克斯椭球边值问题，也可称为椭球改正．2 斯托克斯边值问题2.1　斯托克斯球边值问题解对线性化的斯托克斯边值问题边界面为球近似的情况下，其形式可写为ΔT=0    (在S外);-∂T/∂r-2R-1T|S=Δg    (在S上);T=O(r-1)    (在无穷远处), (9)式中S表示{r=R}平均球面．根据式(8)，此时重力异常可写为Δg=∑n=2∞(n-1)∑m=-nnR-1c¯nmWnm(θ,λ)，(10)那么c¯nm=R4π(n-1)∬σΔgWnm(θ,λ)dσ．(11)c¯nm为边值问题(11)扰动位T的球谐系数，即知道{r=R}上的重力异常Δg即可求出球谐系数．2.2　斯托克斯椭球边值问题解关于椭球边值问题的解算，有很多学者进行了研究，可以分为两种方法：第一种方法是直接计算第二类勒让德函数，获得椭球谐系数，然后转换到球谐系数[9,15]；第二种方法是将椭球面边值问题在椭球扁率平方的量级下转换成球边值问题，进行迭代求解[14,16]．如果采用第一种方法，须要考虑第二类Legendre函数的计算，关于第二类Legendre函数计算，相较于第一类Legendre函数而言，计算较为复杂，因此这里选择后者．对于参考椭球面，有三轴不一致的椭球和旋转椭球两个表示方式[9,15,17]．然而三轴不一致的椭球须要计算拉梅函数较为复杂，因此选择使用旋转椭球作为边界面，其形式为x2+y2a2+z2b2=1，(12)式中a和b分别为椭球的长、短半轴，并记Σ为旋转椭球面．那么斯托克斯椭球边值问题(7)可以写成ΔT=0    (在Σ外);-∂T∂h+1γ∂γ∂hT|Σ=Δg    (在Σ上);T=O(r-1)    (在无穷远处), (13)式中h为椭球面的法线方向．为了便于描述，令ε2=E2/b2，即椭球的第二偏心率，(u,θ,λ)为椭球坐标分量，分别表示椭球的短半轴、归化余纬和经度．事实上，对斯托克斯椭球边界条件-∂T∂h+1γ∂γ∂hT|Σ=Δg进行解算总是保留ε2量级，而忽略掉ε4以上的量级．在椭球坐标中椭球面Σ可以表示为u=b的形式，因此在保留ε2量级下∂T/∂h可写为[14]：∂T∂hΣ=∂T∂r-ε2sinθcosθ∂TR∂θΣ；(14)∂T∂hΣ=1+12ε2sin2θ∂T∂uΣ．(15)式(14)表示球坐标系的形式，式(15)表示椭球坐标系下的形式．对于1γ∂γ∂hΣ有1γ∂γ∂hΣ=1+12ε2sin2θ1γ∂γ∂uΣ=-1+12ε2sin2θε2sin2θb+db, (16)式中d=2b2/a2+2ω2b3/βGM，其中：βGM为地球引力常数；ω为地球自转角速度．结合式(13)和(15)，在ε2量级下，边界条件-∂T∂h+1γ∂γ∂hT|Σ=Δg可以写成-∂T∂h+1γ∂γ∂hT|Σ=-∂T∂r+ε2sinθcosθ∂TR∂θ-1+12ε2sin2θε2sin2θb+dbT|Σ=Δg. (17)最后须要将椭球面Σ延拓到{r=R}的球面S上．记Q是Σ上任意一点，而Q1是Q沿径向方向在平均球面{r=R}上的对应点，那么有QQ1=rQ-rQ1=b2+E2(1+ε2)sin2θ1+ε2sin2θ1/2-R ．(18)保留ε2量级，忽略掉ε4以上量级，则有QQ1=-13+12sin2θRε2．(19)通过泰勒展开公式，略去O(Tε4)小量之后有：T|Q=T|Q1+∂T∂rQ1-13+12sin2θRε2；(20)∂T∂rQ=∂T∂rQ1+∂2T∂r2Q1-13+12sin2θRε2．(21)结合式(20)和(21)，式(17)可以写成-∂T∂h+1γ∂γ∂hT|Σ=-1+db-13+12sin2θRε2∂T∂r-∂2T∂r2-13+12sin2θRε2+ε2sinθcosθ∂TR∂θ-ε2sin2θb+db+dε2sin2θ2bT|S=Δg. (22)因此斯托克斯椭球边值问题可以写成：ΔT=0    (在S外);    -1+db-13+12sin2θRε2∂T∂r-∂2T∂r2-13+12sin2θRε2+ε2sinθcosθ∂TR∂θ-ε2sin2θb+db+dε2cos2θ2bT|S=Δg    (在S上);T=O(r-1)    (在无穷远处).(23)在边值问题(23)中能够发现ε2量级都被保留，ε4量级都被舍去．比较边值问题(9)和(23)，如果在(23)中忽略掉ε2量级，那么斯托克斯椭球边值问题就变成球边值问题．为了求解椭球面下的位系数，须要将边界条件表示成级数的形式，下面将推导边界条件的级数表示形式．根据式(8)，椭球面上扰动位及其偏导数可以写成如下级数形式：T=∑n=2∞∑m=-nnc¯nmWnm(θ,λ)；(24)∂T∂r=-1R∑n=2∞∑m=-nn(n+1)c¯nmWnm(θ,λ)；(25)∂2T∂r2=1R2∑n=2∞∑m=-nn(n+1)(n+2)c¯nmWnm(θ,λ) ；(26)∂T∂θ=-sinθ∑n=2∞∑m=-nnc¯nmdWnm(θ,λ)dcosθ．(27)将式(24)~(27)代入到边界条件(22)中有        -1+db-13+12sin2θRε2∂T∂r-∂2T∂r2-13+12sin2θRε2+ε2sinθcosθ∂TR∂θ-ε2sin2θb+db+dε2sin2θ2bT|S=∑n=2∞∑m=-nnknc¯nmWnm(θ,λ)+∑n=2∞∑m=-nn(ln+hnsin2θ)ε2c¯nmWnm(θ,λ)-ε2sin2θcosθR∑n=2∞∑m=-nnc¯nmdWnm(θ,λ)dcosθ, (28)式中：kn=(n+1)/R-d/b；ln=(n+1)-d3b+n+23R；hn=d(n+1)2b-(n+1)(n+2)2R-2+d2b．引入sin2θP¯nmcosθ=-αnmP¯n+2mcosθ+(1-βnm)P¯nmcosθ-γnmP¯n-2mcosθ; (29)sin2θcosθdP¯nmcosθdcosθ=ξnmP¯n+2mcosθ+(φnm+τnm)P¯nmcosθ+ηnmP¯n-2mcosθ, (30)式中：αnm=12n+3∙(n-m+1)(n-m+2)(n+m+1)(n+m+2)(2n+1)(2n+5)1/2βnm=2n2-2m2+2n-1(2n+3)(2n-1)；γnm=12n-1∙(n+m)(n+m-1)(n-m-1)(n-m)(2n-3)(2n+1)1/2；ξnm=-n2n+3∙(n-m+2)(n+m+1)(n+m+2)(n-m+1)(2n+1)(2n+5)1/2;φnm=n(m-n-1)(n+m+1)(2n+1)(2n+3)；τnm=(n+1)(n+m)(n-m)(2n-1)(2n+1)；ηnm=n+12n-1∙(n-m-1)(n-m)(n+m)(n+m-1)(2n+1)(2n-3)1/2．将式(29)和(30)代入到式(28)中，那么有-1-db13-12sin2θRε2∂T∂r+∂2T∂r213-12sin2θRε2+ε2sinθcosθ∂TR∂θ-ε2sin2θb+db+dε2sin2θ2bT|S=∑n=2∞∑m=-nnknc¯nmWnm(θ,λ)+∑n=2∞∑m=-nn(1-βnm)ln+hn-φnm+τnmRε2c¯nmWnm(θ,λ)+∑n=4∞∑m=-nn-αn-2mln-2-ξn-2mRε2c¯n-2mWnm(θ,λ)+∑n=0∞∑m=-nn-γn+2mln+2-ηn+2mRε2c¯n+2mWnm(θ,λ)=Δg.(31)对椭球边值问题(23)进行求解，在式(31)中首先忽略掉ε2量级，有c¯(1)nm=14πkn∬σΔgWnm(θ,λ)dσ．(32)顾及ε2量级下，有c¯nm(2)=c¯nm(1)-kn-1(1-βnm)ln+hn-(φnm+τnm)/Rε2c¯nm(1)+kn-1αn-2mln-2+ξn-2m/R∙ε2c¯n-2m(1)+kn-1γn+2mln+2+ηn+2m/Rε2c¯n+2m(1). (33)总的来说在O(Tε4)量级下对边值问题(23)的球谐系数解只须按照式(32)和(33)计算即可．3 模拟计算选取WGS84为参考椭球，其长半轴为a=6.378 137×106 m，扁率为f=1/298.257 223 563，参考椭球的引力常数为βGM=3.986 600 441 8×1014 m3/s2，与地球总质量相等，自转角速度为ω=7.292 115×10-5 rad/s，与地球自转角速度相等．选取EGM2008模型的前2 159阶次模型为真实重力场．根据参考椭球的正常重力位计算相应的大地水准面，并计算大地水准面上的重力异常．设置模拟计算的大地水准面上的重力异常为5'×5'格网．椭球面上的正常重力位以及正常重力按照Somigliana公式给出，其形式为[18]U=GMEarctanEu+12ω2a2qq0cos2θ-13+12(u2+E2)sin2θ,式中：q=121+3u2E2arctanEu-3uE；q0=q|u=b．对于式(31)也可以写成级数形式，即参考椭球外部任意点正常重力位为U=GMr1+∑n=1∞J¯2nRr2nP¯2ncosθ+12ω2sin2θ，式中J¯2n为椭球谐系数．令u=b，那么椭球面上的正常重力位为U0=GMEarctanEb+13ω2a2．椭球面上的正常重力为γ=GMaa2cos2θ+b2sin2θ∙1+ω2a2EGMq0'q012cos2θ-16-ω2a2bGMsin2θ．大地水准面上任意点的重力位为V=GMr1+∑n=2∞Rrn∑m=0nC¯nm08Wnmθ,λ+12ω2r2sin2θ,式中：C¯nm08和S¯nm08为EGM2008模型的球谐系数；r为大地水准面到地心的距离；R为平均球半径，根据V=U0即可以得出r．大地水准面上的重力g为g=1rsin∂V∂λ2+1r∂V∂θ2+∂V∂r21/2．大地水准面上的重力异常为Δg=g-γ，大地水准面上的扰动位为T=GM∑n=2∞∑m=-nnRrn+1ΔC¯nmWnm(θ,λ)，式中：ΔC¯nm=C¯nm-J¯2n (m≥0)；ΔS¯nm=Snm (m0)，其中C¯nm，S¯nm为待求的模型系数．对ΔC¯nm，ΔS¯nm将根据斯托克斯球或椭球边值问题，分别依据式(11)、(32)和(33)进行计算，然后求出C¯nm和S¯nm模型系数．采用阶方差表示恢复出重力场位系数精度，即σn=12n+1∑m=0n[(Cnm-Cnmegm)2+(Snm-Snmegm)2]1/2，式中：{Cnm,Snm}表示采用斯托克斯边值问题恢复的球谐系数；{Cnmegm,Snmegm}表示参考重力场模型的球谐系数．首先在观测数据无误差的情况下，直接采用大地水准面上的重力异常，分别按照球和椭球斯托克斯边值问题反演重力场位系数，结果如图1所示．可以看出用斯托克斯球边值问题和椭球边值问题均能用于完整的恢复2 159地球重力场位系数，并且斯托克斯椭球边值问题恢复重力场模型精度要更高，特别在低阶部分精度可以提高到一个量级以上．此外，目前对于超高阶地球重力场模型构建(如2 159阶)线性化理论已经满足要求，无须考虑非线性影响．如果进一步考虑非线性项的影响，理论上会提高重力场模型精度，鉴于非线性项的复杂性暂不考虑非线性项的计算．10.13245/j.hust.230303.F001图1大地水准面上5'×5'格网上的重力对格网上的重力异常Δg分别添加均值0，标准差分别为1和3 mGal的随机误差，表示观测数据的精度．根据斯托克斯球边值问题和椭球边值问题恢复重力场模型位系数的阶方差如图2和3所示．10.13245/j.hust.230303.F002图2大地水准面上5'×5'格网上的重力异常加1 mGal的随机误差10.13245/j.hust.230303.F003图3大地水准面上5'×5'格网上的重力异常加3 mGal的随机误差可以看出：在不同随机误差下经过椭球改正所恢复的重力场模型位系数的精度在中低阶范围内要比球近似解精度要高．随着误差的增大斯托克斯球边值问题所恢复的最大阶数分别为1 600和1 400阶，椭球边值问题恢复重力模型最大阶数分别为 2 000和1 450阶，即随着误差的增大斯托克斯球边值问题和椭球边值问题所能恢复的重力场模型最大阶数都在减少，并且随着观测噪声的增大椭球边值问题所恢复的重力场模型的最大阶数减少最快，球边值问题所恢复出的重力场模型的最大阶数减少相对较慢．这主要由式(33)导致的，即在解算斯托克斯椭球边值问题中，在式(33)中要得出c¯nm(2)需要知道c¯nm(1)，c¯n-2m(1)和c¯n+2m(1)的值，而通过式(32)得出c¯nm(1)，c¯n-2m(1)和c¯n+2m(1)的值是包含误差．在式(33)中存在nε2项，随着n的增大，当nε21时就抵消了由椭球扁率带来ε2量级的精度，所以椭球变值问题所恢复的重力场模型最大阶数下降比较快．为了更加详细阐明该观点，下面将给出理论上的说明．假设大地水准面上第i个格网未加入随机误差时平均重力异常为Δgi，每一个格网加入的随机误差为εi，由式(27)可知c¯nm(1)*=R24π(n-1)∬σ(Δg+εi)Wnm(θ,λ)dσ=R24π(n-1)∬σΔgWnm(θ,λ)dσ+R24π(n-1)∬σεiWnm(θ,λ)dσ,式中c¯nm(1)*为带有随机误差εi所恢复的重力场模型系数．为了简化表示形式，令ϕnm*=R24π(n-1)∬σεiWnm(θ,λ)dσ，并且在n很大时可以认为ϕnm*，ϕn-2m*和ϕn+2m*近似相等．那么c¯nm(1)*可以写成c¯nm(1)*=c¯nm(1)+ϕnm*．由式(33)可知：c¯nm(2)*=c¯nm(2)+ϕnm*-kn-1[(1-βnm)ln+hn-(φnm+τnm)/R]ϕnm*+(αn-2mln-2+R-1ξn-2m)ε2ϕn-2m*+kn-1(γn+2mln+2+ηn+2m/R)ε2ϕn+2m*.在这里令Wnm*=ϕnm*-kn-1[(1-βnm)ln+hn-(φnm+τnm)/R]ϕnm*+kn-1(αn-2mln-2+R-1ξn-2m)ε2ϕn-2m*+kn-1+(γn+2mln+2+ηn+2m/R)ε2ϕn+2m*.当n很大时，Wnm*可以简写成Wnm*≈ϕnm*-nε2(ϕnm*/3+ϕn-2m*/2+ϕn+2m*/2)．可见随着n增大，nε2抵消了椭球边值问题所带来的优点，即Wnm*ϕnm*，即随着误差的增大椭球边值问题所恢复的重力场模型在高阶部分的精度将低于球边值问题．4 结论从模拟计算结果可见：目前构建超高阶地球重力场模型(2 159阶次)线性化理论已经满足要求，无须考虑非线性项的影响，经过椭球改正所恢复出的重力场模型要优于球边值问题，特别在中低阶部分精度大约要高出一个量级．对于数据观测精度与恢复重力场模型阶数的关系，给模拟的大地水准面上的观测值分别添加1和3 mGal的随机误差，可见随着数据精度的降低，经过球和椭球边值问题所恢复的重力场模型阶数降低，在中低阶部分经过椭球改正的精度要优于球边值问题，但中高阶部分椭球边值问题抗误差能力要比球边值问题差，即经椭球改正反演的重力场位系数具有更高的噪声，并对这一点进行了理论上的说明．总体而言，在构建2 159阶次超高阶重力场模型过程中，舍去的非线性项不影响重力场模型精度，经过椭球改正的精度要优于球边值问题．数据观测精度影响构建重力模型的阶数，并且椭球边值问题的抗误差干扰能力要弱于球边值问题．未来构建超高阶边值问题时，高阶部分可使用球边值问题进行替代．