CN106960111B - 一种引潮位的Doodson规格化展开及其精度评定方法 - Google Patents

一种引潮位的Doodson规格化展开及其精度评定方法 Download PDF

Info

Publication number
CN106960111B
CN106960111B CN201710238447.6A CN201710238447A CN106960111B CN 106960111 B CN106960111 B CN 106960111B CN 201710238447 A CN201710238447 A CN 201710238447A CN 106960111 B CN106960111 B CN 106960111B
Authority
CN
China
Prior art keywords
doodson
calculation
tidal
expansion
level
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201710238447.6A
Other languages
English (en)
Other versions
CN106960111A (zh
Inventor
雷伟伟
郑红晓
张煊赫
孙茜
张捍卫
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Henan University of Technology
Original Assignee
Henan University of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Henan University of Technology filed Critical Henan University of Technology
Priority to CN201710238447.6A priority Critical patent/CN106960111B/zh
Publication of CN106960111A publication Critical patent/CN106960111A/zh
Application granted granted Critical
Publication of CN106960111B publication Critical patent/CN106960111B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation

Abstract

本发明公开了一种引潮位的Doodson规格化展开及其精度评定方法，包括建立天体对地球上测站点的引潮位计算函数模型，引潮位计算函数模型由Doodson规格因子展开，构建球面天文学中黄道坐标系、赤道坐标系和时角坐标系间的函数表达式，数据结构与算法设计，“伪波”与“滤波”处理及，计算结果与精度分析等六步。本发明一方面数据计算过程简便，数据运算过程通用性强，数据运算过程规范性和通用性好，便于对数据计算方法的掌握和交流，另一方面有效的克服了传统引潮力计算过程中缺乏检核条件，无法快速确定计算过程与结果的运算精度进行评定，同时在运算过程中也引入了干扰性数据清楚步骤，从而进一步提高了对引潮位数据计算的精度。

Description

一种引潮位的Doodson规格化展开及其精度评定方法
技术领域
本发明涉及一种引潮位的Doodson规格化展开及其精度评定方法，属测绘技术领域。
背景技术
引潮位的展开是地球物理学和大地测量学的基本理论问题之一，高精度的引潮位展开表可为各类地面与空间测量数据的归算、整理和后处理等工作提供重要的参考依据。1921年，Doodson首先基于Brown月球历表和Newcomb太阳历表得到了包含378项展开式的引潮位展开表。在展开过程中，Doodson将展开式表达为“大地系数”和“潮波分量”两部分，为了使“大地系数”在各阶次中的数值保持相对稳定，Doodson定义了一组规格化因子，进行了所谓的“Doodson规格化”处理。随后郗钦文、Cartwright&Tayler、Tamura、Hartmann&Wenzel、Kudryavtsev等人分别基于不同的历表与展开方法，得到了项数各异的引潮位展开表，并被ICET(International Centre for Earth Tide，国际固体潮中心)及IERS(International Earth Rotation and Reference Systems Service，国际地球自转与参考系服务)规范推荐使用，但这些数据分析时，往往数据分析的运算量大，数据运算方式规范性不高，且在计算过程中往往缺乏对干扰性数据缺乏处理，并对运算结构缺乏计算精度评定，从而导致了当前的引潮位计算精度相对较低，因此针对这一问题，迫切需要开发一种全新的引潮位的Doodson规格化展开及其精度评定方法，以满足实际工作的需要。
发明内容
本发明的目的是要提供一种引潮位的Doodson规格化展开及其精度评定方法。
为达到上述目的，本发明是按照以下技术方案实施的：
一种引潮位的Doodson规格化展开及其精度评定方法，包括如下步骤：
第一步，建立天体对地球上测站点的引潮位计算函数模型，根据待计算天体位置、观测点所在地球位置参数为基础，构建天体对地球上测站点的引潮位计算函数模型为：
Figure BDA0001268662910000021
其中，
Figure BDA0001268662910000022
GMJ为万有引力常数与天体J的质量之积；
(αJ,δJ,RJ)、
Figure BDA0001268662910000023
分别表示天体、测站点在国际地球参考系中的地心经度、地心纬度、地心距
ZJ为天体与测站之间的地心天顶距，
Pn(x)为n阶勒让德函数，
HJ为天体地方时角；
第二步，引潮位计算函数模型由Doodson规格因子展开，首先对Doodson规格因子进行定位，然后由Doodson规格因子先展开关于cosZJ的多项式，然后将根据三角函数倍角公式，把与cos HJ幂相关的项转换为cos mHJ角的形式，且将含有cos mHJ的项进行同类项合并，最后得引潮位计算函数模型的Doodson规格化展开到表达式为：
Figure BDA0001268662910000024
其中，
Figure BDA0001268662910000031
分别称为“大地系数”和“潮波分量”；
Figure BDA0001268662910000032
分别为Pn(cos ZJ)的第n阶展开式含cos mHJ的项中与
Figure BDA0001268662910000033
δJ有关的函数项；
Figure BDA0001268662910000034
是在计算过程中产生的常数系数，
m为0和正整数倍；
由于
Figure BDA0001268662910000035
的值域范围各不相同，为了使“大地系数”
Figure BDA0001268662910000036
的数值在不同阶次中保持相对稳定，Doodson规格化因子中引入因子
Figure BDA0001268662910000037
使
Figure BDA0001268662910000038
然后再由Doodson规格化公式形式进行引潮位的展开工作，由此得到：
Figure BDA0001268662910000039
引潮位计算函数模型的Doodson规格化展开到表达式变形为：
Figure BDA00012686629100000310
其中，
Figure BDA00012686629100000311
并基于Doodson规格化的引潮位展开即将
Figure BDA00012686629100000312
表达为如下的“潮波”形式：
Figure BDA00012686629100000313
其中：
Figure BDA00012686629100000314
为无量纲的数值，
ki为整数，可由ki组合得到Doodson编码，
τ、s、h、p、N′、ps为Doodson定义的日月天文辐角参数；
第三步，构建球面天文学中黄道坐标系、赤道坐标系和时角坐标系间的函数表达式，基于ELP/MPP02月球历表和Newcomb太阳历表中，分别将月球、太阳在地心天球中的黄经λ、黄纬β、地心距R的三角函数表达式
Figure BDA0001268662910000041
其中，对应正弦，R对应余弦，D、F、l、l′为Delaunay天文辐角参数，得到Delaunay天文辐角参数与Doodson天文辐角参数之间的关系为：
Figure BDA0001268662910000042
记
Figure BDA0001268662910000043
(对于月球
Figure BDA0001268662910000044
对于太阳
Figure BDA0001268662910000045
)，
Figure BDA0001268662910000046
则有：
Figure BDA0001268662910000047
Figure BDA0001268662910000048
Figure BDA0001268662910000049
则可基于ELP/MPP02和Newcomb历表，经过简单转换后，将xJ、βJ、cJ/RJ三者均表达为形如
Figure BDA00012686629100000410
的三角函数级数的形式，其中振幅均为无量纲的数值，且xJ、βJ对应正弦，cJ/RJ对应余弦；
由此得到球面天文学中黄道坐标系、赤道坐标系和时角坐标系间的关系式为：
Figure BDA00012686629100000411
其中，由于无法基于天文历表直接计算cos mHJ，在计算过程中，需要将cos mHJ按照倍角公式展开，并与cosδJ相乘后得到cosζJ的各次幂，然后再基于
Figure BDA00012686629100000412
进行展开计算；
第四步，数据结构与算法设计，在将
Figure BDA00012686629100000413
展开为潮波的过程，核心是sinδJ、cosζJ、cJ/RJ以及三者各次幂的计算，同时由于：
sinδJ、cosζJ是由sinβJ、sinλJ、cosβJ、cosλJ计算得到；
sinλJ、cosλJ是由sin xJ、cos xJ计算得到；
sinβJ、cosβJ、sin xJ、cos xJ则是由βJ、xJ以及二者的各次幂计算得到；
同时，由于在天文历表中，βJ、xJ、cJ/RJ均可表达由三角函数级数
Figure BDA0001268662910000051
进行表述，因此sinδJ、cosζJ、sinλJ、cosλJ、sinβJ、cosβJ、sin xJ、cos xJ的计算可归结为三角函数乘法和加法的计算工作，因此在进行计算时，首先定义一个数组[A k1 k2 k3 k4 k5 k6 q]，用来代表三角函数A sin/cos(k1τ+k2s+k3h+k4p+k5N′+k6ps)，
其中数组中最后一个元素q为正弦或余弦的标志，当为正弦时q＝1，余弦时q＝-1。则基于该数据结构，两个三角函数相加的算法可定义为上述两个数组的并列，而两个三角函数相乘的算法可根据三角函数积化和差公式表达为两个数组之和，于计算过程中产生大量的数组，为了提高计算效率，在三角函数相加的算法中还需要进行三角函数同类项的合并工作，且需考虑到sin(-θ)＝-sinθ、cos(-θ)＝cosθ这两种特殊情况，其中的乘法算法部分可根据三角函数积化和差公式得到，并由此获得引潮位Doodson规格化展开的计算流程；
第五步，“伪波”与“滤波”处理，在第四步的计算过程中，由于截断的原因，将会产生“伪波”现象，且不同周期的潮波体现在cosmH项中，H表示天体的时角，周期为一日，因此，在计算过程中，需要进行“滤波”处理，即对于包含cosmH项的展开式而言，在展开计算过程中，对于τ的系数不等于m的项需要删除掉，将产生的“伪波”予以消除；
第六步，计算结果与精度分析，欲使引潮位展开表达到10-11ms-2(1ngal)的精度水平，即测站点在球坐标系下的径向引潮力|gr,S|≥10-11ms-2，由于：
Figure BDA0001268662910000052
鉴于ELP/MPP02月球历表和Newcomb太阳历表的截断精度，需将月球、太阳的引潮位分别展开至5阶、3阶，振幅
Figure BDA0001268662910000061
绝对值的截断阈值为10-7；
因此在展开计算时，理论公式中相关天文、测地常数均采用IERS 2010规范推荐值，并在将计算过程中产生的“伪波”进行“滤波”处理后，最终得到一个展开式的引潮位展开表，
然后将Doodson、XI89、RATGP95、HW95以及第四步和第五步中的主要潮波项的振幅大小进行了对比，由于HW95原始展开表理论公式没有采用Doodson规格化，因此需要将其转换为Doodson规格化形式，HW95展开表理论公式为：
Figure BDA0001268662910000062
由于在6阶以内
Figure BDA0001268662910000063
各个阶次极值绝对值均在1～3.606之间，可直接以式
Figure BDA0001268662910000064
为基础进行展开，在HW95展开表中，潮波系数
Figure BDA0001268662910000065
的单位为m2s-2，
由于展开方法、天体历表、展开阶次以及振幅截断阈值的不同，对展开表中主要潮波振幅的比较并不能反映各个展开表的精度水平，对于引潮位展开表的精度评定，目前国际上通用做法是采用Wenzel提出的基于引潮力基准序列的精度评定方法，由于：
Figure BDA0001268662910000066
为在更大时间范围内评价各个引潮位展开表的精度水平，根据JPL最新发布的DE431历表，计算若干时间段内基准序列BFDE431，然后基于Doodson、XI89、RATGP95、HW95展开表以及本文得到的展开表，分别计算各展开表相应的法向引潮力序列值，并求得与BFDE431基准序列间的差值序列，统计结果列表汇总中，并对结果进行汇总比较，并当Doodson的计算值与XI89、RATGP95、HW95中的XI89数值接近时证明计算结果精度达到要求，为否则再次返回到第四部进行数据运算。
进一步的，所述的第六步中，计算若干时间段内基准序列BFDE431时，时间段选择单位为年，且时间段的时间跨度不小于10年，在选定的时间跨度范围内，计算时的数据的时间间隔为1小时。
进一步的，所述的第六步中，计算若干时间段内基准序列BFDE431的具体计算方法为：
1)基于DE431历表计算得到各天体在地心天球参考系中的直角坐标；
2)通过基于春分点的岁差章动转换方法或是基于无旋转原点的转换方法[10]，得到各天体在国际地球参考系中的直角坐标，并转换为球面坐标(αJ,δJ,RJ)；
3)由(αJ,δJ)计算cos ZJ的值，并基于勒让德函数递推算法计算Pn(x)各阶的值[11]；
4)根据式(19)计算得到径向引潮力gr,S的具体数值，并转换为法向引潮力gr，相关天文常数采用DE431历表头文件中提供的数值。
本发明一方面数据计算过程简便，数据运算过程通用性强，数据运算过程规范性和通用性好，便于对数据计算方法的掌握和交流，另一方面有效的克服了传统引潮力计算过程中缺乏检核条件，无法快速确定计算过程与结果的运算精度进行评定，同时在运算过程中也引入了干扰性数据清楚步骤，从而进一步提高了对引潮位数据计算的精度。
附图说明
下面结合附图和具体实施方式来详细说明本发明
图1：为本发明方法流程图；
图2：为天体、测站在地心天球中的位置示意图；
图3：自定义加法算法的流程图；
图4：引潮位计算函数模型Doodson规格化展开流程图；
图5：各个引潮位展开表主要振幅项对比表；
图6：各个引潮位展开表主要振幅项对比表；
图7：引潮力与BFDE431间的差值序列图(Doodson)；
图8：引潮力与BFDE431间的差值序列图(XI89)；
图9：引潮力与BFDE431间的差值序列图(RATGP95)；
图10：引潮力与BFDE431间的差值序列图(HW95)；
图11：引潮力与BFDE431间的差值序列图(由第四步、第五步计算得到)
具体实施方式
为使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体实施方式，进一步阐述本发明。
如图1-2所示，一种引潮位的Doodson规格化展开及其精度评定方法，包括如下步骤：
第一步，建立天体对地球上测站点的引潮位计算函数模型，根据待计算天体位置、观测点所在地球位置参数为基础，构建天体对地球上测站点的引潮位计算函数模型为：
Figure BDA0001268662910000081
其中，
Figure BDA0001268662910000082
GMJ为万有引力常数与天体J的质量之积；
(αJ,δJ,RJ)、
Figure BDA0001268662910000083
分别表示天体、测站点在国际地球参考系中的地心经度、地心纬度、地心距
ZJ为天体与测站之间的地心天顶距，
Pn(x)为n阶勒让德函数，
HJ为天体地方时角；
第二步，引潮位计算函数模型由Doodson规格因子展开，首先对Doodson规格因子进行定位，然后由Doodson规格因子先展开关于cos ZJ的多项式，然后将根据三角函数倍角公式，把与cos HJ幂相关的项转换为cos mHJ角的形式，且将含有cos mHJ的项进行同类项合并，最后得引潮位计算函数模型的Doodson规格化展开到表达式为：
Figure BDA0001268662910000091
其中，
Figure BDA0001268662910000092
分别称为“大地系数”和“潮波分量”；
Figure BDA0001268662910000093
分别为Pn(cos ZJ)的第n阶展开式含cos mHJ的项中与
Figure BDA0001268662910000094
δJ有关的函数项；
Figure BDA0001268662910000095
是在计算过程中产生的常数系数，
m为0和正整数倍；
由于
Figure BDA0001268662910000096
的值域范围各不相同，为了使“大地系数”
Figure BDA0001268662910000097
的数值在不同阶次中保持相对稳定，Doodson规格化因子中引入因子
Figure BDA0001268662910000098
使
Figure BDA0001268662910000099
然后再由Doodson规格化公式形式进行引潮位的展开工作，由此得到：
Figure BDA00012686629100000910
引潮位计算函数模型的Doodson规格化展开到表达式变形为：
Figure BDA00012686629100000911
其中，
Figure BDA00012686629100000912
并基于Doodson规格化的引潮位展开即将
Figure BDA00012686629100000913
表达为如下的“潮波”形式：
Figure BDA0001268662910000101
其中：
Figure BDA0001268662910000102
为无量纲的数值，
ki为整数，可由ki组合得到Doodson编码，
τ、s、h、p、N′、ps为Doodson定义的日月天文辐角参数；
第三步，构建球面天文学中黄道坐标系、赤道坐标系和时角坐标系间的函数表达式，基于ELP/MPP02月球历表和Newcomb太阳历表中，分别将月球、太阳在地心天球中的黄经λ、黄纬β、地心距R的三角函数表达式
Figure BDA0001268662910000103
其中，对应正弦，R对应余弦，D、F、l、l′为Delaunay天文辐角参数，得到Delaunay天文辐角参数与Doodson天文辐角参数之间的关系为：
Figure BDA0001268662910000104
记
Figure BDA0001268662910000105
(对于月球
Figure BDA0001268662910000106
对于太阳
Figure BDA0001268662910000107
)，
Figure BDA0001268662910000108
则有：
Figure BDA0001268662910000109
Figure BDA00012686629100001010
Figure BDA00012686629100001011
则可基于ELP/MPP02和Newcomb历表，经过简单转换后，将xJ、βJ、cJ/RJ三者均表达为形如
Figure BDA00012686629100001012
的三角函数级数的形式，其中振幅均为无量纲的数值，且xJ、βJ对应正弦，cJ/RJ对应余弦；
由此得到球面天文学中黄道坐标系、赤道坐标系和时角坐标系间的关系式为：
Figure BDA0001268662910000111
其中，由于无法基于天文历表直接计算cos mHJ，在计算过程中，需要将cos mHJ按照倍角公式展开，并与cosδJ相乘后得到cosζJ的各次幂，然后再基于
Figure BDA0001268662910000112
进行展开计算；
第四步，数据结构与算法设计，在将
Figure BDA0001268662910000113
展开为潮波的过程，核心是sinδJ、cosζJ、cJ/RJ以及三者各次幂的计算，同时由于：
sinδJ、cosζJ是由sinβJ、sinλJ、cosβJ、cosλJ计算得到；
sinλJ、cosλJ是由sin xJ、cos xJ计算得到；
sinβJ、cosβJ、sin xJ、cos xJ则是由βJ、xJ以及二者的各次幂计算得到；
同时，由于在天文历表中，βJ、xJ、cJ/RJ均可表达由三角函数级数
Figure BDA0001268662910000114
进行表述，因此sinδJ、cosζJ、sinλJ、cosλJ、sinβJ、cosβJ、sin xJ、cos xJ的计算可归结为三角函数乘法和加法的计算工作，因此在进行计算时，首先定义一个数组[A k1 k2 k3 k4 k5 k6 q]，用来代表三角函数A sin/cos(k1τ+k2s+k3h+k4p+k5N′+k6ps)，
其中数组中最后一个元素q为正弦或余弦的标志，当为正弦时q＝1，余弦时q＝-1。则基于该数据结构，两个三角函数相加的算法可定义为上述两个数组的并列，而两个三角函数相乘的算法可根据三角函数积化和差公式表达为两个数组之和，于计算过程中产生大量的数组，为了提高计算效率，在三角函数相加的算法中还需要进行三角函数同类项的合并工作，且需考虑到sin(-θ)＝-sinθ、cos(-θ)＝cosθ这两种特殊情况，其中的乘法算法部分可根据三角函数积化和差公式得到，并由此获得引潮位Doodson规格化展开的计算流程；
第五步，“伪波”与“滤波”处理，在第四步的计算过程中，由于截断的原因，将会产生“伪波”现象，且不同周期的潮波体现在cosmH项中，H表示天体的时角，周期为一日，
以月球为例：
cos H＝cos(τ+s-180°-L)
＝cos(τ+υ)
＝A cosτ+B sinτ (14)
因此，当m＝1时，展开结果为周日波；同理，理论上cos2H＝A′cos2τ+B′sin2τ，即当m＝2时，展开结果理论上为半日波，不可能出现长周期波。而在计算中，cos2H无法直接计算得到，实际采用的是以下公式：
cos2H＝2cos2H-1
＝2(A cosτ+B sinτ)2-1
＝(A2-B2)cos2τ+2AB sin2τ+A2+B2-1 (15)
由于计算过程中的截断，A2+B2-1≠0，即在半日波中产生了长周期波，这就是产生了所谓的“伪波”。同理，在cos3H、cos4H、cos5H项的展开中，均会产生周期各异的“伪波”。
因此，在计算过程中，需要进行“滤波”处理，即对于包含cos mH项的展开式而言，在展开计算过程中，对于τ的系数不等于m的项需要删除掉，将产生的“伪波”予以消除；
第六步，计算结果与精度分析，欲使引潮位展开表达到10-11ms-2(1ngal)的精度水平，即测站点在球坐标系下的径向引潮力|gr,S|≥10-11ms-2，由于：
Figure BDA0001268662910000121
鉴于ELP/MPP02月球历表和Newcomb太阳历表的截断精度，需将月球、太阳的引潮位分别展开至5阶、3阶，振幅
Figure BDA0001268662910000122
绝对值的截断阈值为10-7；
因此在展开计算时，理论公式中相关天文、测地常数均采用IERS 2010规范推荐值，并在将计算过程中产生的“伪波”进行“滤波”处理后，最终得到一个展开式的引潮位展开表，
然后将Doodson、XI89、RATGP95、HW95以及第四步和第五步中的主要潮波项的振幅大小进行了对比，由于HW95原始展开表理论公式没有采用Doodson规格化，因此需要将其转换为Doodson规格化形式，HW95展开表理论公式为：
Figure BDA0001268662910000131
由于在6阶以内
Figure BDA0001268662910000132
各个阶次极值绝对值均在1～3.606之间，可直接以式
Figure BDA0001268662910000133
为基础进行展开，在HW95展开表中，潮波系数
Figure BDA0001268662910000134
的单位为m2s-2，
由于展开方法、天体历表、展开阶次以及振幅截断阈值的不同，对展开表中主要潮波振幅的比较并不能反映各个展开表的精度水平，对于引潮位展开表的精度评定，目前国际上通用做法是采用Wenzel提出的基于引潮力基准序列的精度评定方法，由于：
Figure BDA0001268662910000135
以德国Black Forest Observatory(BFO)观测站(L＝8.3300°E、
Figure BDA0001268662910000136
r＝6366836.969m)为例，首先根据美国JPL发布的DE数值历表，基于式(19)计算得到某段时间范围内时间间隔为1小时的径向引潮力gr,S时间序列，并将径向结果转换为法向结果[16]，得到法向引潮力gr时间序列，命名为引潮力的BFDE基准序列；然后根据引潮位展开表及式(16)，结合各个天文辐角计算公式，计算得到BFO测站相应的法向引潮力序列，将其与BFDE基准序列间差值序列的统计结果作为评定引潮位展开表精度的指标。Wenzel曾根据DE200、DE403历表分别给出了1987年～1994年及2017年～2024年间的BFDE200、BFDE403两个引潮力基准序列，并嵌入固体潮调合分析软件Eterna中。
为在更大时间范围内评价各个引潮位展开表的精度水平，根据JPL最新发布的DE431历表，本文计算得到了1950年～2050年间，时间间隔1小时的gr基准序列BFDE431，步骤如下：
1)基于DE431历表计算得到各天体在地心天球参考系中的直角坐标；
2)通过基于春分点的岁差章动转换方法或是基于无旋转原点的转换方法[10]，得到各天体在国际地球参考系中的直角坐标，并转换为球面坐标(αJ,δJ,RJ)；
3)由(αJ,δJ)计算cos ZJ的值，并基于勒让德函数递推算法计算Pn(x)各阶的值[11]；
4)根据式(19)计算得到径向引潮力gr,S的具体数值，并转换为法向引潮力gr，相关天文常数采用DE431历表头文件中提供的数值。
然后基于Doodson、XI89、RATGP95、HW95展开表以及本文得到的展开表，分别计算各展开表相应的法向引潮力序列值，并求得与BFDE431基准序列间的差值序列，统计结果列表汇总中，从计算结果可知，HW95展开表精度最高，差值序列的数值在±1.3×10-11ms-2以内，均方差达到0.1×10-11ms-2量级，这是由于HW95展开表是将月球、太阳的引潮位分别展开至6阶、3阶，水星、金星、火星、木星、土星展开至2阶，并考虑了地球扁率的影响，基于DE200历表采用频谱分析方法获得的结果；其次是RATGP95展开表，该展开表虽然是将月球、太阳的引潮位分别展开至5阶、3阶，但在天文历表计算中，考虑了行星对月球轨道的摄动效应以及地球扁率的影响，因此也获得了较高的精度；精度最低的是Doodson展开表；而本文计算得到的展开表在差值序列极值方面稍微优于XI89展开表，而均方差基本相等，因此，本文展开表和XI89展开表基本上精度一致。
本发明一方面数据计算过程简便，数据运算过程通用性强，数据运算过程规范性和通用性好，便于对数据计算方法的掌握和交流，另一方面有效的克服了传统引潮力计算过程中缺乏检核条件，无法快速确定计算过程与结果的运算精度进行评定，同时在运算过程中也引入了干扰性数据清楚步骤，从而进一步提高了对引潮位数据计算的精度。
本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (3)

1.一种引潮位的Doodson规格化展开及其精度评定方法，其特征在于，精密引潮力的计算及其影响因素分析方法包括如下步骤：
第一步，建立天体对地球上测站点的引潮位计算函数模型，根据待计算天体位置、观测点所在地球位置参数为基础，构建天体对地球上测站点的引潮位计算函数模型为：
Figure FDA0002473914420000011
n为阶数；
其中，
Figure FDA0002473914420000012
GMJ为万有引力常数与天体J的质量之积；
(αJ,δJ,RJ)、
Figure FDA0002473914420000013
分别表示天体、测站点在国际地球参考系中的地心经度、地心纬度、地心距
ZJ为天体与测站之间的地心天顶距，
Pn(x)为n阶勒让德函数，
HJ为天体地方时角；
第二步，引潮位计算函数模型由Doodson规格因子展开，首先对Doodson规格因子进行定位，然后由Doodson规格因子先展开关于cosZJ的多项式，然后将根据三角函数倍角公式，把与cosHJ幂相关的项转换为cosmHJ角的形式，且将含有cosmHJ的项进行同类项合并，最后得引潮位计算函数模型的Doodson规格化展开到表达式为：
Figure FDA0002473914420000014
其中，
Figure FDA0002473914420000021
分别称为“大地系数”和“潮波分量”，m为次，a为椭球长半径，cJ为天体平均地心距，a/cJ为天体的地平视差，DJ为Doodson常数；
Figure FDA0002473914420000022
分别为Pn(cos ZJ)的第n阶展开式含cos mHJ的项中与
Figure FDA0002473914420000023
δJ有关的函数项，
Figure FDA0002473914420000024
是在计算过程中产生的常数系数；
由于
Figure FDA0002473914420000025
的值域范围各不相同，为了使“大地系数”
Figure FDA0002473914420000026
的数值在不同阶次中保持相对稳定，Doodson规格化因子中引入因子
Figure FDA0002473914420000027
使
Figure FDA0002473914420000028
然后再由Doodson规格化公式形式进行引潮位的展开工作，由此得到：
Figure FDA0002473914420000029
引潮位计算函数模型的Doodson规格化展开到表达式变形为：
Figure FDA00024739144200000210
其中，
Figure FDA00024739144200000211
Figure FDA00024739144200000212
为
Figure FDA00024739144200000213
除以
Figure FDA00024739144200000214
的商，并基于Doodson规格化的引潮位展开即将
Figure FDA00024739144200000215
表达为如下的“潮波”形式：
Figure FDA00024739144200000216
其中：
Figure FDA00024739144200000217
为无量纲的数值，ki为整数，可由ki组合得到Doodson编码，τ、s、h、p、N′、ps为Doodson定义的日月天文辐角参数；
第三步，构建球面天文学中黄道坐标系、赤道坐标系和时角坐标系间的函数表达式，基于ELP/MPP02月球历表和Newcomb太阳历表中，分别将月球、太阳在地心天球中的黄经λ、黄纬β、地心距R的三角函数表达式
Figure FDA00024739144200000218
其中，λ、β对应正弦，R对应余弦，D、F、l、l′为Delaunay天文辐角参数，得到Delaunay天文辐角参数与Doodson天文辐角参数之间的关系为：
Figure FDA0002473914420000031
记
Figure FDA0002473914420000032
(对于月球
Figure FDA0002473914420000033
对于太阳
Figure FDA0002473914420000034
)，
Figure FDA00024739144200000312
则有：
Figure FDA0002473914420000035
Figure FDA0002473914420000036
Figure FDA0002473914420000037
则可基于ELP/MPP02和Newcomb历表，经过简单转换后，将xJ、βJ、cJ/RJ三者均表达为形如
Figure FDA0002473914420000038
的三角函数级数的形式，其中振幅均为无量纲的数值，且xJ、βJ对应正弦，cJ/RJ对应余弦；
由此得到球面天文学中黄道坐标系、赤道坐标系和时角坐标系间的关系式为：
Figure FDA0002473914420000039
其中，由于无法基于天文历表直接计算cosmHJ，在计算过程中，需要将cosmHJ按照倍角公式展开，并与cosδJ相乘后得到cosζJ的各次幂，然后再基于
Figure FDA00024739144200000310
进行展开计算；
第四步，数据结构与算法设计，在将
Figure FDA00024739144200000311
展开为潮波的过程，核心是sinδJ、cosζJ、cJ/RJ以及三者各次幂的计算，同时由于：
sinδJ、cosζJ是由sinβJ、sinλJ、cosβJ、cosλJ计算得到；
sinλJ、cosλJ是由sinxJ、cosxJ计算得到；
sinβJ、cosβJ、sinxJ、cosxJ则是由βJ、xJ以及二者的各次幂计算得到；
同时，由于在天文历表中，βJ、xJ、cJ/RJ均可表达由三角函数级数
Figure FDA0002473914420000041
进行表述，因此sinδJ、cosζJ、sinλJ、cosλJ、sinβJ、cosβJ、sinxJ、cosxJ的计算可归结为三角函数乘法和加法的计算工作，因此在进行计算时，首先定义一个数组[A k1 k2 k3 k4 k5 k6 q]，用来代表三角函数A sin/cos(k1τ+k2s+k3h+k4p+k5N′+k6ps)，
其中数组中最后一个元素q为正弦或余弦的标志，当为正弦时q＝1，余弦时q＝-1则基于该数据结构，两个三角函数相加的算法可定义为上述两个数组的并列，而两个三角函数相乘的算法可根据三角函数积化和差公式表达为两个数组之和，于计算过程中产生大量的数组，为了提高计算效率，在三角函数相加的算法中还需要进行三角函数同类项的合并工作，且需考虑到sin(-θ)＝-sinθ、cos(-θ)＝cosθ这两种特殊情况，其中的乘法算法部分可根据三角函数积化和差公式得到，并由此获得引潮位Doodson规格化展开的计算流程；
第五步，“伪波”与“滤波”处理，在第四步的计算过程中，由于截断的原因，将会产生“伪波”现象，且不同周期的潮波体现在cos mH项中，H表示天体的时角，周期为一日，因此，在计算过程中，需要进行“滤波”处理，即对于包含cos mH项的展开式而言，在展开计算过程中，对于τ的系数不等于m的项需要删除掉，将产生的“伪波”予以消除；
第六步，计算结果与精度分析，欲使引潮位展开表达到10-11ms-2(1ngal)的精度水平，即测站点在球坐标系下的径向引潮力|gr,S|≥10-11ms-2，由于：
Figure FDA0002473914420000042
其中V(t)为引潮位计算函数模型；
鉴于ELP/MPP02月球历表和Newcomb太阳历表的截断精度，需将月球、太阳的引潮位分别展开至5阶、3阶，振幅
Figure FDA0002473914420000043
绝对值的截断阈值为10-7；
因此在展开计算时，理论公式中相关天文、测地常数均采用IERS 2010规范推荐值，并在将计算过程中产生的“伪波”进行“滤波”处理后，最终得到一个展开式的引潮位展开表，
然后将Doodson、XI89、RATGP95、HW95以及第四步和第五步中的主要潮波项的振幅大小进行了对比，由于HW95原始展开表理论公式没有采用Doodson规格化，因此需要将其转换为Doodson规格化形式，HW95展开表理论公式为：
Figure FDA0002473914420000051
由于在6阶以内
Figure FDA0002473914420000052
各个阶次极值绝对值均在1～3.606之间，可直接以式
Figure FDA0002473914420000053
为基础进行展开，在HW95展开表中，潮波系数
Figure FDA0002473914420000054
的单位为m2s-2，
由于展开方法、天体历表、展开阶次以及振幅截断阈值的不同，对展开表中主要潮波振幅的比较并不能反映各个展开表的精度水平，对于引潮位展开表的精度评定，目前国际上通用做法是采用Wenzel提出的基于引潮力基准序列的精度评定方法，由于：
Figure FDA0002473914420000055
为在更大时间范围内评价各个引潮位展开表的精度水平，根据JPL最新发布的DE431历表，计算若干时间段内基准序列BFDE431，然后基于Doodson、XI89、RATGP95、HW95展开表以及本文得到的展开表，分别计算各展开表相应的法向引潮力序列值，并求得与BFDE431基准序列间的差值序列，统计结果列表汇总中，并对结果进行汇总比较，并当Doodson的计算值与XI89、RATGP95、HW95中的XI89数值接近时证明计算结果精度达到要求，为否则再次返回到第四部进行数据运算。
2.根据权利要求1所述的一种引潮位的Doodson规格化展开及其精度评定方法，其特征在于，所述的第六步中，计算若干时间段内基准序列BFDE431时，时间段选择单位为年，且时间段的时间跨度不小于10年，在选定的时间跨度范围内，计算时的数据的时间间隔为1小时。
3.根据权利要求1所述的一种引潮位的Doodson规格化展开及其精度评定方法，其特征在于，所述的第六步中，计算若干时间段内基准序列BFDE431的具体计算方法为：
1)基于DE431历表计算得到各天体在地心天球参考系中的直角坐标；
2)通过基于春分点的岁差章动转换方法或是基于无旋转原点的转换方法[10]，得到各天体在国际地球参考系中的直角坐标，并转换为球面坐标(αJ,δJ,RJ)；
3)由(αJ,δJ)计算cosZJ的值，并基于勒让德函数递推算法计算Pn(x)各阶的值[11]；
4)根据式(19)计算得到径向引潮力gr,S的具体数值，并转换为法向引潮力gr，相关天文常数采用DE431历表头文件中提供的数值。
CN201710238447.6A 2017-04-13 2017-04-13 一种引潮位的Doodson规格化展开及其精度评定方法 Expired - Fee Related CN106960111B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710238447.6A CN106960111B (zh) 2017-04-13 2017-04-13 一种引潮位的Doodson规格化展开及其精度评定方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710238447.6A CN106960111B (zh) 2017-04-13 2017-04-13 一种引潮位的Doodson规格化展开及其精度评定方法

Publications (2)

Publication Number Publication Date
CN106960111A CN106960111A (zh) 2017-07-18
CN106960111B true CN106960111B (zh) 2020-07-31

Family

ID=59483440

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710238447.6A Expired - Fee Related CN106960111B (zh) 2017-04-13 2017-04-13 一种引潮位的Doodson规格化展开及其精度评定方法

Country Status (1)

Country Link
CN (1) CN106960111B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110727914B (zh) * 2019-09-30 2022-12-02 中国人民解放军战略支援部队信息工程大学 一种基于向量运算的垂线偏差单点计算方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105203104A (zh) * 2015-09-16 2015-12-30 北京航空航天大学 一种适用于高精度惯导系统的重力场建模方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105203104A (zh) * 2015-09-16 2015-12-30 北京航空航天大学 一种适用于高精度惯导系统的重力场建模方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
不同规格化的引潮位展开及其转换;郗钦文;《地球物理学报》;20070131;第50卷(第1期);182-194 *
引潮位展开的不同规格化形式及其转换;雷伟伟 等;《大地测量与地球动力学》;20161215;第36卷(第12期);1105-1108 *
引潮位的Doodson 规格化展开及其精度评定;雷伟伟 等;《测绘科学技术学报》;20170503;第34卷(第1期);5-9 *
引潮位的计算机数值解析展开研究;闫朋远;《中国优秀硕士学位论文全文数据库 基础科学辑》;20130315;A012-1 *
精密引潮力的计算及其影响因素分析;雷伟伟 等;《地球物理学进展》;20161015;第31卷(第5期);1917-1926 *

Also Published As

Publication number Publication date
CN106960111A (zh) 2017-07-18

Similar Documents

Publication Publication Date Title
Haigh et al. Global influences of the 18.61 year nodal cycle and 8.85 year cycle of lunar perigee on high tidal levels
Seo et al. Terrestrial water mass load changes from Gravity Recovery and Climate Experiment (GRACE)
Petrie et al. Higher‐order ionospheric effects on the GPS reference frame and velocities
Sjöberg et al. Gravity inversion and integration
Konopliv et al. High‐resolution lunar gravity fields from the GRAIL primary and extended missions
Swinbank et al. Compilation of wind data for the Upper Atmosphere Research Satellite (UARS) reference atmosphere project
Johnson et al. Ocean bottom pressure seasonal cycles and decadal trends from GRACE Release‐05: Ocean circulation implications
Fang et al. Empirical cotidal charts of the Bohai, Yellow, and East China Seas from 10 years of TOPEX/Poseidon altimetry
Crétaux et al. Seasonal and interannual geocenter motion from SLR and DORIS measurements: Comparison with surface loading data
Vinogradova et al. Relation between sea level and bottom pressure and the vertical dependence of oceanic variability
Hill et al. Combination of geodetic observations and models for glacial isostatic adjustment fields in Fennoscandia
Stopa et al. Sea state trends and variability: Consistency between models, altimeters, buoys, and seismic data (1979–2016)
Müller et al. Toward an internal gravity wave spectrum in global ocean models
Guo et al. Improvements in the monthly gravity field solutions through modeling the colored noise in the GRACE data
Zhu et al. Assimilation of coastal acoustic tomography data using an unstructured triangular grid ocean model for water with complex coastlines and islands
Mendillo et al. Simultaneous ionospheric variability on Earth and Mars
Knudsen et al. Correcting GRACE gravity fields for ocean tide effects
Shim et al. CEDAR‐GEM challenge for systematic assessment of ionosphere/thermosphere models in predicting TEC during the 2006 December storm event
Li et al. Ellipsoidal correction in GRACE surface mass change estimation
Hopkins et al. Storm impact on morphological evolution of a sandy inlet
Ghobadi‐Far et al. A transfer function between line‐of‐sight gravity difference and GRACE intersatellite ranging data and an application to hydrological surface mass variation
CN106960111B (zh) 一种引潮位的Doodson规格化展开及其精度评定方法
Huang et al. Seasonal behavior of the semidiurnal and diurnal tides, and mean flows at 95 km, based on measurements from the High Resolution Doppler Imager (HRDI) on the Upper Atmosphere Research Satellite (UARS)
Valle‐Levinson et al. Solar activity and lunar precessions influence extreme sea‐level variability in the US Atlantic and Gulf of Mexico coasts
Kleinherenbrink et al. Trends and interannual variability of mass and steric sea level in the T ropical A sian S eas

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20200731