本发明公开了一种极区扰动重力无奇异性详细计算模型及其建模方法，属于大地测量领域
背景技术：
在地球物理问题的求解中，为了使数学模型更接近于研究对象的物理特性，坐标形式的选择往往取决于所研究对象边界面的几何形状。因此，在扰动重力的表示中，更倾向于采用地心球坐标(r,θ,λ)，这使得扰动重力的计算变得更加简单实用。
扰动重力三分量的地心球坐标表达式为
其中，δρ、δλ分别表示扰动重力的径向、纬向和经向分量；fM表示地球重力常数；ρ表示计算点的地心向径；θ和λ分别表示计算点的余纬和地心经度；R表示地球平均半径；表示完全正常化地球扰动重力位系数；n和m分别表示球谐系数的阶和次；表示完全正常化缔合Legendre函数。
公式(2)中的的递推计算公式为
式中是kroneker符号，当m＝0时，当m≠0时，
将(4)式代入(2)式，就可以得到扰动重力的传统计算模型。
从公式(2)-(4)可以看出，当飞行器飞越两极及其附近地区时，由于分母上的sinθ＝0或者接近于0，会使得扰动重力的2个分量的计算出现无穷大。然而，地球两极及其附近地区作为地球表面上的普通区域，扰动重力的计算值应该与其它地球表面一样是有限值，不同区域，量级不会太大的差距。因此，这就是极区扰动重力计算存在的θ奇异性问题。随着重力测量技术的发展，从理论上彻底解决球坐标下极区扰动重力计算存在的θ奇异性问题就显得越来越重要。
为了解决极区扰动重力计算存在的θ奇异性问题，传统方法主要是直角坐标表示法。
直角坐标表示法的基本原理是指将扰动重力直接表示成地心直角坐标的形式，因而分母上不会含有等于0或接近于0的项，避免了扰动重力计算存在的θ奇异性问题。但也由此引入了其他的问题：公式形式和参数比较复杂，对递推计算的稳定性和精度也有较大影响。故本专利提出一种极区扰动重力无θ奇异性的详细计算模型。
技术实现要素：
本发明提供一种极区扰动重力无θ奇异性的详细计算模型，即从Legendre函数的基本性质出发，研究得出Legendre函数对θ的一阶导数和一类多项式的无奇异性计算公式，再将其代入扰动重力的传统计算模型中，同时充分考虑了m等于0、1以及其它量时的情况，最终建立了极区扰动重力无θ奇异性的详细计算模型。该发明使得扰动重力计算模型更接近于研究对象的物理特性，并且公式形式简单，稳定性和实用性较好，彻底解决了两极地区扰动重力计算公式存在的θ奇异性难题，为飞行器飞越两极及其附近地区时的扰动重力计算奠定了理论基础，也可利用该发明提供的模型计算全球任何区域的扰动重力结果。
参阅图7，本发明提出的一种极区扰动重力无奇异性详细计算模型的建立过程主要步骤如下：
第一步：计算得到Legendre函数对θ的一阶导数的无奇异性结果。
用表示伴随Legendre函数，其中Pn(x)是n阶Legendre多项式，将x用cosθ来代替，可得
(5)式两端对θ求导，得
因为，Legendre多项式Pn(x)满足下列关系
结合(6)式和(7)式，可得
对(8)式中的Legendre函数进行规格化处理，得到Legendre函数对θ的一阶导数的无奇异性计算公式为：
第二步：计算得到一类多项式的无奇异性结果。
Legendre函数满足下列伴随Legendre方程
由此可得
在(8)式两端继续对θ求导，得
将(8)式代入(13)式，可以得到Legendre函数二阶导数的计算公式
将(8)式和(14)式代入(12)式，可得
根据Legendre函数的标准前向行推法，有
因此，
将(17)式和(18)式代入(15)式，可得
对(20)式中的Legendre函数进行规格化处理，得到多项式无奇异性计算公式为：
第三步：将上述公式代入扰动重力各分量的传统计算模型中，并充分考虑m等于0、1以及其它量时的情况，建立极区扰动重力无θ奇异性详细计算模型。
将(9)、(10)、(20)式代入(1)-(3)式，可得到扰动重力各分量的无θ奇异性计算公式为：
公式(21)-(23)即为扰动重力在极区的无θ奇异性详细计算模型。
在对本发明新建模型有效性进行测试基础上，本发明与现有技术相比体现出如下优点：
a.彻底解决了极区扰动重力计算公式存在的θ奇异性难题；
b.新建立的扰动重力各分量无θ奇异性详细计算模型，公式形式简单，稳定性和实用性较好；
c.新建立的扰动重力无θ奇异性详细计算模型，为飞行器飞越两极及其附近地区时的扰动重力计算奠定了理论基础，也可利用该发明提供的模型计算全球任何区域的扰动重力结果。
附图说明：
图1是实施例中南极地区扰动重力分量的等值线图；
图2是实施例中南极地区扰动重力δλ分量的等值线图；
图3是实施例中南极地区扰动重力δρ分量的等值线图；
图4是实施例中北极地区扰动重力分量的等值线图；
图5是实施例中北极地区扰动重力δλ分量的等值线图；
图6是实施例中北极地区扰动重力δρ分量的等值线图；
图7是本发明提出的极区扰动重力无θ奇异性详细计算模型建立基本流程。
具体实施方式：
假设飞行器飞行在距离地面300km的被动段轨道上，利用360阶EGM2008地球重力场模型，计算了南、北极圈(66.5°)以内30′×30′格网大小的扰动重力格网点值。两极地区扰动重力各分量的统计结果如表1、2所示：
表1南极地区扰动重力三分量统计结果(mGal)
表2北极地区扰动重力三分量统计结果(mGal)
从表1、2可以看出，扰动重力分量计算结果的标准差最小，因此其结果的变化幅度较小，δλ分量次之，δρ分量变化幅度最大，这是由于地球重力场对扰动重力δρ分量影响最大；从各分量的平均值来看，δλ分量沿格林尼治子午线的对称性最好。
为了更加直观地分析扰动重力各分量的变化趋势，本文绘出了南、北极圈以内扰动重力三分量的等值线图。
图1～6反映了扰动重力各分量的变化趋势，可以看出，扰动重力各分量基本上沿格林尼治子午线对称分布，其中δλ分量的对称性最好，这与表1、2的分析结果一致；北极地区比南极地区的扰动重力各分量的对称性要好，这是因为北极地区陆地面积很少，地形变化比较平缓的原因造成的；δρ分量比δλ分量变化更加剧烈，这与前面表1和表2的分析结果一致；各分量在陆地部分比海洋部分的变化更加剧烈，这是因为该地区陆地上地势起伏较大，地形情况复杂，而该地区海洋地形变化较为平缓的原因造成的。
图1～3反映出扰动重力各分量的变化比较均匀，这是因为南极地区大部分范围都是陆地，而图4～6反映出扰动重力各分量的变化比较复杂，从北极圈到北极点变化程度由剧烈到平缓，这是因为北极地区经历了从陆地到海洋的变化过程。