一种精密引潮力的计算及其影响因素分析方法与流程

文档序号:11155748阅读:693来源:国知局
一种精密引潮力的计算及其影响因素分析方法与制造工艺

本发明涉及一种精密引潮力的计算及其影响因素分析方法,属测绘技术领域。



背景技术:

精密引潮力的计算是地球物理、空间科学、天文/测地学研究的基本理论问题。张捍卫(2004)根据弹性力学的基本理论,从地球的微小弹性运动方程出发,给出了引潮力的严格数学定义,抛弃了前人对引潮力定义所给定的所有假设条件,回避了引潮力不同定义容易引起的歧义性。Hartmann&Wenzel(1995a,1995b)、Kudryavtsev(2004)基于引潮位理论公式和天体历表,将月球、太阳、近地行星的引潮位分别展开至6阶、3阶、2阶,计算得到引潮位及引潮力的时间序列。地球扁率对引潮力的计算存在影响,郗钦文(1999)和张捍卫(2010)分别指出,地球扁率对引潮力的影响不超过5×10-11ms-2(10-11ms-2=1nGal),这在高精度的潮汐分析、引潮力计算中是必须要考虑的。其实Bartels早在1957年就曾指出这一问题,随后Wilhelm(1983)和Dahlen(1993)给出了相关的计算公式,Hartmann&Wenzel(1995a,1995b)在此基础上,通过对公式中参数的数量级估计,表明当引潮力计算截断阈值为0.001×10-11ms-2时,需要考虑地球扁率对月球和太阳引潮力的影响,对其它行星的影响则不需考虑。通常情况下,引潮力被认为是引潮位的径向梯度值,而重力仪的测量数据是基于铅垂线方向,因此需将径向元素转换为铅垂线法向的元素,Wenzel(1974)指出此处可忽略垂线偏差的影响,只需将基于球面的径向元素转换为基于椭球面的法向元素即可。

在引潮力的计算过程中,天球参考系转换、勒让德函数计算是两个核心工作,在以往的计算中(Hartmann&Wenzel,1995a,1995b;Roosbeek,1996;Kudryavtsev,2004),对于前者采用经典的基于春分点的岁差章动转换,对于后者采用手工推导出各阶次的具体表达式;这种做法存在两个不足:一是缺乏检核条件,无法快速确定计算过程与结果的准确性;二是勒让德函数表达式推导的难度随着阶次升高将会倍增,从而严重影响了引潮力数据计算准确性和可靠性,同时也导致计算过程相对复杂,不利于引潮力计算方法的推广和掌握,同时,传统的引潮力计算过程中,无法直接多的对引潮力计算结果产生干扰情况进行分析,从而进一步增加了传统方法计算引潮力数据误差率,因此针对这一问题,迫切需要开发一种全新的精密引潮力的计算方法,以满足实际工作的需要。



技术实现要素:

本发明的目的是要提供一种精密引潮力的计算及其影响因素分析方法。

为达到上述目的,本发明是按照以下技术方案实施的:

一种精密引潮力的计算及其影响因素分析方法,其特征在于,所述的精密引潮力的计算及其影响因素分析方法包括如下步骤:

第一步,建立天体对地球上测站点的引潮位计算函数模型,根据待计算天体位置、观测点所在地球位置参数为基础,构建天体对地球上测站点的引潮位计算函数模型为:

其中,GMJ为万有引力常数与天体J的质量之积;

RJ、r分别表示天体、测站点的地心距;

ZJ为天体与测站之间的地心天顶距,

Pn(x)为n阶勒让德函数;

第二步,天球参考系转换,首先基于美国喷气推进实验室发布的DE历表进行天体坐标计算工作,以地球为中心星、各天体为目标星计算得到天体坐标,该天体坐标是天体在GCRS坐标,然后并将天体在GCRS中的坐标转换为ITRS中的坐标,在进行GCRS坐标转换为ITRS中的坐标时,首先分别构建基于春分点的岁差章动转换函数模型和基于CIO的无旋转原点转换函数模型,最后根据IERS规范中的GCRS坐标转换为ITRS坐标的转换参数及计算模型,分别引入到基于春分点的岁差章动转换函数模型和基于CIO的无旋转原点转换函数模型中,并分别按照基于春分点的岁差章动转换函数模型和基于CIO的无旋转原点转换函数模型进行数据运算,并将基于春分点的岁差章动转换函数模型和基于CIO的无旋转原点转换函数模型计算的结果进行比对,且当基于春分点的岁差章动转换函数模型和基于CIO的无旋转原点转换函数模型计算结果一致,则证明数据运算正确,从而进行下一步操作,当基于春分点的岁差章动转换函数模型和基于CIO的无旋转原点转换函数模型计算结果不一致时,则证明数据运算不正确;

第三步,勒让德函数及其导数运算,第一步构建的建天体对地球上测站点的引潮位计算函数模型中勒让德函数Pn(x),结合地心纬度参数,建立勒让德函数Pn(sinθ)的递推运算函数:

其中:θ为地心纬度,θ∈[-90°,90°](即对于天体θ=δ、对于测站

n≥2;

同时基于勒让德函数Pn(sinθ),构建勒让德函数Pn(sinθ)对sinθ的一阶导数以及对θ的一阶导数

的递推公式为:

其中和在数据运算时均采用标准向前列递推算法,其运算公式为:

其中n≥2且n>m;

m≥2;

n>m;

m≥0;

并由此得到:

则递推初值为:

第四步,天体对地球上测站点的引潮力计算,首先基于测站点在ITRS中的地心经度、地心纬度得勒让德函数Pn(x)的等价变形式,由球面三角学中边的余弦公式得:

则存在以下关系式

其中为n阶m次完全规格化缔合勒让德函数;

然后将勒让德函数Pn(x)的变形式引入到天体对地球上测站点的引潮位计算函数模型中,并得到天体对地球上测站点的引潮位计算函数模型等价变形式:

其中a为地球参考椭球长半径;

然后根据地球扁率对引潮位的影响计算函数:

其中J2⊕为地球重力场二阶带球谐系数;

最后将地球上测站点的引潮位计算函数模型等价变形式得到的V<1>值和然后根据地球扁率对引潮位的影响计算函数得到的V<2>值求和,以得到天体对测站点总的引力值,其具体表达式为:

则引潮位的梯度值为:

其中:gr,S、gλ,S、分别表示测站点在球坐标系下的径向引潮位径向重力固体潮、径向经度固体潮、径向纬度固体潮;

第五步,天体对测站点总的引力值验算,天体对测站点总的引力值验算检核计算公式:

此时记为测站点处的大地纬度B与球心纬度之差,则测站点在椭球坐标系下的引潮力gr,E、法向经度固体潮gL,E、法向纬度固体潮gB,E验算检核计算公式为:

第六步,数据计算与分析,在进行精密引潮力的计算时,首先需要进行各天体展开阶次及天体取舍,确定运算参数,然后在计算过程中分别对天体坐标转换对计算结果影响、地球扁率对引潮位对计算结果影响、岁差章动模型更新及春分点的岁差章动转换函数模型和基于CIO的无旋转原点转换函数模型运算对计算结果影响、径向法向转换对计算结果影响、参考框架偏差对计算结果影响、计算过程中地球时转换对计算结果影响、极移对计算结果影响及历表更新对计算结果影响进行验证即可。

进一步的,所选择的第二步中进行基于春分点的岁差章动转换函数模型和基于CIO的无旋转原点转换函数模型进行数据计算时,均涉及到了ERA的计算,ERA的计算公式为:ERA(Tu)=2π(0.7790572732640+1.00273781191135448Tu);

其中Tu是从J2000.0起算的儒略日数,时间系统为UT1,这就涉及到TT与UT1之间的转换,记二者的差值为ΔT=TT-UT1,ΔT是基于最小二乘法拟合得到一组六次多项式,其具体表达式为:

进一步的,所述基于春分点的岁差章动转换函数模型为:[ITRS]=W(t)·Reqx(t)·Qeqx(t)·[GCRS];

t为从J2000.0起算的TT儒略世纪数(以下相同),[ITRS]、[GCRS]分别表示相应参考系中的三维直角坐标,Qeqx(t)、Reqx(t)、W(t)分别称为框架偏差-岁差-章动矩阵、地球自转矩阵、极移矩阵;

其中框架偏差-岁差-章动矩阵、地球自转矩阵、极移矩阵表达式为:

B为参考框架偏差矩阵(简称框架偏差),P为岁差矩阵,N为章动矩阵,它们的计算公式分别为:

ξ0、η0表示CIP在GCRS中J2000.0时刻的天极偏差,dα0为瞬时平春分点的赤经偏差值;ψA、ωA为基本赤道岁差参数,χA为黄道岁差导出量,εA为瞬时黄赤交角;Δψ为黄经章动,Δε为交角章动。

进一步的,所述基于CIO的无旋转原点转换函数模型为:[ITRS]=W(t)·RCIO(t)·QCIO(t)·[GCRS];

QCIO(t)、RCIO(t)也分别称为框架偏差-岁差-章动矩阵、地球自转矩阵;

QCIO(t)、RCIO(t)的计算模型分别为:

RCIO(t)=R3(ERA);

其中:s为CIO定位角,表示CIO在CIP赤道上的位置,ERA为地球自转角。

本发明一方面数据计算过程简便。数据运算过程通用性强,便于对数据计算方法的掌握和交流,另一方面有效的克服了传统引潮力计算过程中缺乏检核条件,无法快速确定计算过程与结果的准确性和勒让德函数表达式推导的难度随着阶次升高将会倍增的弊端,于此同时也便于对引潮力计算过程中易出现的干扰因素进行分析排除,从而极大的提高了引潮力计算工作的准确性和可靠性。

附图说明

下面结合附图和具体实施方式来详细说明本发明

图1:为本发明方法流程图;

图2:为天体、测站在地心天球中的位置示意图。

具体实施方式

为使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解,下面结合具体实施方式,进一步阐述本发明。

如图1和2所示,一种精密引潮力的计算及其影响因素分析方法,其特征在于,所述的精密引潮力的计算及其影响因素分析方法包括如下步骤:

第一步,建立天体对地球上测站点的引潮位计算函数模型,根据待计算天体位置、观测点所在地球位置参数为基础,构建天体对地球上测站点的引潮位计算函数模型为:

其中,GMJ为万有引力常数与天体J的质量之积;

RJ、r分别表示天体、测站点的地心距;

ZJ为天体与测站之间的地心天顶距,

Pn(x)为n阶勒让德函数;

第二步,天球参考系转换,首先基于美国喷气推进实验室发布的DE历表进行天体坐标计算工作,以地球为中心星、各天体为目标星计算得到天体坐标,该天体坐标是天体在GCRS坐标,然后并将天体在GCRS中的坐标转换为ITRS中的坐标,在进行GCRS坐标转换为ITRS中的坐标时,首先分别构建基于春分点的岁差章动转换函数模型和基于CIO的无旋转原点转换函数模型,最后根据IERS规范中的GCRS坐标天体坐标转换为ITRS坐标的转换参数及计算模型,分别引入到基于春分点的岁差章动转换函数模型和基于CIO的无旋转原点转换函数模型中,并分别按照基于春分点的岁差章动转换函数模型和基于CIO的无旋转原点转换函数模型进行数据运算,并将基于春分点的岁差章动转换函数模型和基于CIO的无旋转原点转换函数模型计算的结果进行比对,且当基于春分点的岁差章动转换函数模型和基于CIO的无旋转原点转换函数模型计算结果一致,则证明数据运算正确,从而进行下一步操作,当基于春分点的岁差章动转换函数模型和基于CIO的无旋转原点转换函数模型计算结果不一致时,则证明数据运算不正确;

第三步,勒让德函数及其导数运算,第一步构建的建天体对地球上测站点的引潮位计算函数模型中勒让德函数Pn(x),结合地心纬度参数,建立勒让德函数Pn(sinθ)的递推运算函数:

其中:θ为地心纬度,θ∈[-90°,90°](即对于天体θ=δ、对于测站

n≥2;

同时基于勒让德函数Pn(sinθ),构建勒让德函数Pn(sinθ)对sinθ的一阶导数以及对θ的一阶导数

的递推公式为:

其中和在数据运算时均采用标准向前列递推算法,其运算公式为:

其中n≥2且n>m;

m≥2;

n>m;

m≥0;

并由此得到:

则递推初值为:

第四步,天体对地球上测站点的引潮力计算,首先基于测站点在ITRS中的地心经度、地心纬度得勒让德函数Pn(x)的等价变形式,由球面三角学中边的余弦公式得:

则存在以下关系式

其中为n阶m次完全规格化缔合勒让德函数;

然后将勒让德函数Pn(x)的变形式引入到天体对地球上测站点的引潮位计算函数模型中,并得到天体对地球上测站点的引潮位计算函数模型等价变形式:

其中a为地球参考椭球长半径;

然后根据地球扁率对引潮位的影响计算函数:

其中J2⊕为地球重力场二阶带球谐系数;

最后将地球上测站点的引潮位计算函数模型等价变形式得到的V<1>值和然后根据地球扁率对引潮位的影响计算函数得到的V<2>值求和,以得到天体对测站点总的引力值,其具体表达式为:

则引潮位的梯度值为:

其中:gr,S、gλ,S、分别表示测站点在球坐标系下的径向引潮位径向重力固体潮、径向经度固体潮、径向纬度固体潮;

第五步,天体对测站点总的引力值验算,天体对测站点总的引力值验算检核计算公式:

此时记为测站点处的大地纬度B与球心纬度之差,则测站点在椭球坐标系下的引潮力gr,E、法向经度固体潮gL,E、法向纬度固体潮gB,E验算检核计算公式为:

第六步,数据计算与分析,在进行精密引潮力的计算时,首先需要进行各天体展开阶次及天体取舍,确定运算参数,然后在计算过程中分别对天体坐标转换对计算结果影响、地球扁率对引潮位对计算结果影响、岁差章动模型更新及春分点的岁差章动转换函数模型和基于CIO的无旋转原点转换函数模型运算对计算结果影响、径向法向转换对计算结果影响、参考框架偏差对计算结果影响、计算过程中地球时转换对计算结果影响、极移对计算结果影响及历表更新对计算结果影响进行验证即可。

本实施例中,所选择的第二步中进行基于春分点的岁差章动转换函数模型和基于CIO的无旋转原点转换函数模型进行数据计算时,均涉及到了ERA的计算,ERA的计算公式为:ERA(Tu)=2π(0.7790572732640+1.00273781191135448Tu);

其中Tu是从J2000.0起算的儒略日数,时间系统为UT1,这就涉及到TT与UT1之间的转换,记二者的差值为ΔT=TT-UT1,ΔT是基于最小二乘法拟合得到一组六次多项式,其具体表达式为:

本实施例中,所述基于春分点的岁差章动转换函数模型为:[ITRS]=W(t)·Reqx(t)·Qeqx(t)·[GCRS];

t为从J2000.0起算的TT儒略世纪数(以下相同),[ITRS]、[GCRS]分别表示相应参考系中的三维直角坐标,Qeqx(t)、Reqx(t)、W(t)分别称为框架偏差-岁差-章动矩阵、地球自转矩阵、极移矩阵;

其中框架偏差-岁差-章动矩阵、地球自转矩阵、极移矩阵表达式为:

B为参考框架偏差矩阵(简称框架偏差),P为岁差矩阵,N为章动矩阵,它们的计算公式分别为:

ξ0、η0表示CIP在GCRS中J2000.0时刻的天极偏差,dα0为瞬时平春分点的赤经偏差值;ψA、ωA为基本赤道岁差参数,χA为黄道岁差导出量,εA为瞬时黄赤交角;△ψ为黄经章动,△ε为交角章动。

本实施例中,所述基于CIO的无旋转原点转换函数模型为:[ITRS]=W(t)·RCIO(t)·QCIO(t)·[GCRS];

QCIO(t)、RCIO(t)也分别称为框架偏差-岁差-章动矩阵、地球自转矩阵;

QCIO(t)、RCIO(t)的计算模型分别为:

RCIO(t)=R3(ERA);

其中:s为CIO定位角,表示CIO在CIP赤道上的位置,ERA为地球自转角。

本实施例中,第六步中判断天体坐标转换对计算结果影响时,在进行精密引潮力的计算工作之前,首先需要确定参与计算的天体及其展开阶数,基于现阶段重力仪测量精度水平,并兼顾重力仪未来的发展,以|gr,S,J|≥0.001×10-11ms-2(0.001nGal)为计算截断阈值,由于:

取rmax=a,对于x∈[-1,1],|Pn(x)|≤1,对于天体J而言,其第n阶参与计算的条件为:

基于DE431历表,计算得到1000AD~3000AD共计2000年内时间间隔为1天的各天体地心距序列,统计得到最小地心距min(RJ),并根据DE431历表中给出的GMJ数值,即可得到各天体在各阶次引潮力最大值

在计算过程中,理论上引潮力的计算至0.001×10-11ms-2即可,但为了分析各类影响因素对引潮力计算影响的具体量值,实际计算至0.000001×10-11ms-2即可。

本实施例中,第六步中判断天体坐标转换对计算结果影响时,为确保计算过程与结果的准确性,对天体坐标转换,分别采用基于春分点的岁差章动转换、基于CIO的无旋转原点转换两种方法,转换模型中的参数又分别基于IERS 2003、2010规范推荐的四组参数计算模型;对引潮力计算公式,分别采用基于勒让德函数及其导数递推算法、完全规格化缔合勒让德函数及其导数递推算法两种数学模型得到的第四步的天体对地球上测站点的引潮力计算式。

在下列各项影响因素计算过程中,均采用第四步的天体对地球上测站点的引潮力计算式分别独立计算得到V、gr,S、gλ,S、序列,然后由(8)式转换为V、gr,E、gL,E、gB,E序列。计算结果表明:对各个序列而言,两组公式结果间的差值ΔV、Δgr,E、ΔgL,E、ΔgB,E序列均为零,这就从勒让德函数递推计算方面确保了计算过程与计算结果的准确性。

在天体坐标转换方面,计算结果表明:两个规范中转换参数模型间的差异、两种坐标转换方法间的差异对引潮力计算的影响非常微小,可完全忽略不计,这就从天体坐标转换方面确保了计算过程与计算结果的准确性。

本实施例中,第六步中判断地球扁率对引潮位对计算结果影响时,引潮位V<1>是将地球作为均质球体得到的,而真实的地球则是一个两级略扁的旋转椭球体,因此地球扁率就对引潮位的计算带来影响,因此需要计算得到由地球扁率而引起的各天体对测站引潮位与引潮力影响的具体数值,并由此得到由地球扁率引起的月球对测站引潮位与引潮力的影响分别在0.00013m2s-2、1.890×10-11ms-2、1.994×10-11ms-2、1.909×10-11ms-2以内;由地球扁率引起的太阳对测站引潮位与引潮力的影响分别在0.00000011m2s-2、0.002×10-11ms-2、0.002×10-11ms-2、0.002×10-11ms-2以内;对其它天体的影响均为零,因此,在精密引潮力的计算过程中,只需考虑地球扁率对月球、太阳的影响即可。

本实施例中,第六步中判断岁差章动模型更新及春分点的岁差章动转换函数模型和基于CIO的无旋转原点转换函数模型运算对计算结果时,针对两种坐标转换方法,IERS 2003、2010规范分别推荐了两套4组不同的转换参数,为明确岁差章动模型更新、不同坐标转换方法间的差异对引潮力计算的影响,分别基于不同的岁差章动模型和坐标转换方法,计算得到相应序列,并分别设为基于2010A、2010B、2003A、2003B四种坐标转换模型,计算得到V、gr,E、gL,E、gB,E四组序列值,差值序列即反映了岁差章动模型更新、不同坐标转换方法间差异对引潮位与引潮力计算的影响,并得到岁差章动模型更新2010A-2003A、2010B-2003B对引潮位与引潮力的影响分别在2×10-8m2s-2、0.00061×10-11ms-2、0.00047×10-11ms-2、0.00040×10-11ms-2以内;不同坐标转换方法2010A-2010B、2003A-2003B的影响量值分别在2×10-10m2s-2、0.000005×10-11ms-2、0.000004×10-11ms-2、0.000003×10-11ms-2以内;

这说明岁差章动模型更新、两种坐标转换方法间的差异对引潮位与引潮力计算的影响甚微,实际计算时可完全忽略不计。在引潮位与引潮力的计算中,可采用任意岁差章动模型及坐标转换方法。

本实施例中,第六步中判断径向法向转换对计算结果影响时,计算得到的gr,S、gλ,S、是径向梯度结果,需根据将其转换为法向结果。若忽略径向到法向的转换,就对引潮位与引潮力的计算带来影响。

若忽略径向到法向的转换,对V和gL,E的数值没有影响,对gr,E、gB,E影响的最大值分别达到476×10-11ms-2、500×10-11ms-2,该影响量值远远超出允许的精度范围,即计算结果错误。因此,在引潮力计算中,一定不能忽略将径向结果转换为法向结果,否则结果错误。

本实施例中,第六步中判断数据运算参考数据源框架对计算结果影响时,参考框架偏差是由于在标准历元J2000.0时刻岁差章动模型给出的天极、分点与GCRS的天极、经度零点不重合造成的。IERS 2003、2010规范中考虑了参考框架偏差对天体坐标转换的影响,但在IERS 1996规范中并没有考虑该影响。

以2010A模型为例,计算了在基于春分点的岁差章动转换过程中,当忽略参考框架偏差时,对引潮位与引潮力的影响量值分别在1×10-6m2s-2、0.032×10-11ms-2、0.032×10-11ms-2、0.021×10-11ms-2以内。

本实施例中,第六步中判断计算过程中地球时转换对计算结果影响时,在计算ERA时,如果误将时间尺度用作TT,即忽略了地球时与世界时的差值ΔT,就会对引潮位与引潮力计算结果带来影响。

如果忽略地球时与世界时间的转换,对引潮位及引潮力的影响量值分别达到0.03m2s-2、916×10-11ms-2、1144×10-11ms-2、696×10-11ms-2,该影响量值远远超出允许的精度范围,即计算结果错误。

因此,在引潮力计算中,一定不能忽略地球时与世界时转换的影响,否则结果错误。

本实施例中,第六步中判断极移对计算结果影响时,由于极移量xp、yp的具体数值目前无法准确预测,只能根据EOP 08C04模型中已发布的相关数据进行计算。中国天文年历在进行天体坐标转换时,并没有考虑极移的影响,实际上是将天体坐标转换至TIRS。由于在两种坐标转换方法中,极移影响与岁差章动模型无关,因此极移对坐标转换的影响量对两种坐标转换方法、两个规范而言均是相同的。本文计算了当忽略极移时,对引潮位与引潮力的影响量值分别在2×10-5m2s-2、0.586×10-11ms-2、0.543×10-11ms-2、0.341×10-11ms-2以内。

本实施例中,第六步中判断历表更新对计算结果影响时,随着天体测量理论与技术的进步,尤其是精度更高的天体测量数据的积累,JPL不断的对DE历表进行更新,发布了一系列版本的历表,其中比较著名的有DE200、DE403、DE405、DE406、DE413、DE421、DE423等,目前最新的是2014年发布的DE432历表。在IERS 2003、2010规范中,分别将DE405、DE421历表作为ICRS的动力学实现。

在后续的引潮位频谱法展开研究中,采用DE431历表,其覆盖时间范围为30000年(13200BC~17191AD)。为明确历表更新对引潮力的影响,分别基于DE405、DE421、DE431历表,计算了相应的序列,并统计出各序列间的差异。

历表更新对引潮位与引潮力的影响量值分别在4×10-7m2s-2、0.012×10-11ms-2、0.013×10-11ms-2、0.013×10-11ms-2以内。

历表更新对引潮力的影响有两方面原因:一是由于历表自身精度的提高;二是不同历表采用的天文常数系统不同,在DE405、421、431三个历表中,各天体GM值以及定义距离的AU(Astronomical Unit,天文单位)都是不同的。

本发明一方面数据计算过程简便。数据运算过程通用性强,便于对数据计算方法的掌握和交流,另一方面有效的克服了传统引潮力计算过程中缺乏检核条件,无法快速确定计算过程与结果的准确性和勒让德函数表达式推导的难度随着阶次升高将会倍增的弊端,于此同时也便于对引潮力计算过程中易出现的干扰因素进行分析排除,从而极大的提高了引潮力计算工作的准确性和可靠性。

本行业的技术人员应该了解,本发明不受上述实施例的限制,上述实施例和说明书中描述的只是说明本发明的原理,在不脱离本发明精神和范围的前提下,本发明还会有各种变化和改进,这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

当前第1页1 2 3 
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1