第３４卷第４期 ２０１４年７月 海洋测 绘 Ｖｏ１．３４．Ｎｏ．４ ＨＹＤＲ０ＧＲＡＰＨＩＣ ＳＵＲＶＥＹＩＮＧ Ａ ＮＤ ＣＨＡＲＴＩＮＧ Ｊｕ１．，２０１４ ＤＯＩ：１０．３９６９／ｊ．ｉｓｓｎ．１６７１—３Ｏ ．２０１４．０４．００２ 超高阶地球重力场模型的高程异常优化算法 曲政豪，曾 添，周 强 （解放军信息工程大学地理空间信息学院，河南郑州４５００５２） 摘要：通过分析比较标准前向列推、标准前向行推、跨阶次递推和Ｂｅｌｉｋｏｖ列推４种缔合勒让德递推算法的精 度、稳定性以及计算速度，提出了选取Ｂｅｌｉｋｏｖ列推法来解算超高阶重力场模型高程异常；研究探讨了基于严密球 谐级数展开、保留泰勒级数展开的零阶项和保留至泰勒级数展开一阶项计算模型高程异常的三种算法，并进行了 实验计算分析。结果表明，保留至泰勒级数一阶项的模型高程异常既能保证计算速度也能达到足够的精度，可满 足大区域高分辨率高程异常建模的需求。 关键词：缔合勒让德函数；超高阶地球重力场模型；高程异常；泰勒级数展开 中图分类号：Ｐ２２３ 文献标志码：Ａ 文章编号：１６７１＿３Ｏ４Ｊ４（２０１４）０４—０００５—０４ 解超高阶的缔合勒让德函数，为此需要优化计算方 法，使得精度、稳定性以及效率有所提高，这样才能 基于超高阶地球重力场模型和已知某一块格网 在较短的时间内获得更稳定更精确的结果。由于直 的高程数据求解格网模型高程异常，当利用严密球 谐级数展开模型时，需要对每个格网点求解缔合勒 让德函数和地心向径，运算量庞大，花费时间长，非 接计算法所需时间太长，因此这种算法不予考虑。 下面通过比较现有的几种缔合勒让德函数递推计算 方法…：标准前向列推法、标准前向行推法、跨阶次 递推法和Ｂｅｌｉｋｏｖ列推法，从计算时间和计算精度两 方面分析以比较不同方法的优劣。 ２．１ 算法模型 常不实用。而如果不考虑地形数据的影响，即只计 算每个纬度圈里的缔合勒让德函数，虽然计算速度 提高很多，但求解的精度不是很理想，无法满足高程 异常计算精度的需要。那么有没有一种既使精度损 标准前向列推法可表示为 ２ ： Ｐ （ｃｏｓ０）＝口 ｃｏｓＯ Ｐ 一１． （ｃｏｓＯ）一 ｂ Ｐ ． 失不大，又能提高程序运行时间的方法呢？本文就 此问题展开了探讨，一个重要思想就是把地球重力 场模型高程异常看成以地心向径为自变量的函数， 在高程等于零处按泰勒公式进行展开，并保留其一 阶项。这样求解勒让德函数时就不用代人每个格网 点的高程，大大节省了运算时间，同时又较好地提高 了精度。由于计算模型高程异常时，对不同纬度圈 其缔合勒让德函数都必须重新求解，每一次求解运 （ｃｏｓ０），（Ｖｎ≤２）（１） 标准前向行推法 ： （ｃ０ ＝√ （ｇ￣ｍＣＯｔＯ ｌ（ｃ伽）一 ｈ Ｐ ． ＋ （ｃｏｓ０）） （２） 算量都较大，为此本文也探讨了标准前向列推、标准 前向行推、跨阶次、Ｂｅｌｉｋｏｖ列推４种不同勒让德递 推算法的运算速度和精度稳定问题，最后综合取舍 当ｍ＝ｎ时，采用下式计算 ］： Ｐ （ｃｏｓ０）＝Ｃｍｓｉｎ０ Ｐ ＿ｌ＿ 一１（ｃｏｓＯ），（Ｖｍ＞１） （３） 给出了较为适宜超高阶模型计算的一种算法，即 Ｂｅｌｉｋｏｖ列推法 跨阶次递推法 ： Ｐ， （ｃｏｓ０）＝ａ Ｐ 一２． ｃｏｓ０）＋／３ Ｐ 一２． 一２ ｃｏｓ０）一 Ｐ ２ ｃｏｓ０） （４） ２ 不同缔合勒让德递推算法的分析 为了计算超高阶重力场模型高程异常，必须求 收稿日期：２０１３—１２—１７；修回日期：２０１４—０４－２５ 基金项目：国家自然科学基金（４１２７４０２９）；国家８６３计划（２０１３ＡＡ１２２５０２）。 作者简介：曲政豪（１９９２一），男，河南新密人，硕士研究生，主要从事大地测量研究。 ６ 海洋测 绘 第３４卷 对于Ｂｅｌｉｋｏｖ列推法，其递推式为 ｃｏｓ０）一 Ｐ棚（ｃｏｓ０）＝ｔＰ ｏ．： 误差也比较稳定，但跨阶次递推法的计算时间较长， 当要计算大范围数据时，选择Ｂｅｌｉｋｏｖ列推法较为 适宜。 。詈 ， （ｃ。ｓ ），ｍ＝０） Ｐ （ｃｏｓＯ）＝ｔＰ 一１ ｃｏｓ０）一 “标准『ＤｅｌｉＫＯ ｖ瓤 Ｉ Ｕ砸偿 Ｉ 踌阶次递推法 标准向前列推法Ｉ『［丢 （ｃｏｓ ）一 ｎ－ｌ，ｍ－ｌ（ｃｏｓ ）］， 一一 （ｍ＞０） （５） 再将其转化为相应的正常化系数，式中Ⅳ 表 达参见文献［１］。 一～ ～ ———————————一Ｐ （ｃｏｓ０）＝￣／２ｎ＋１ Ｎ Ｐ （ｃｏｓ０） （６） ２．２ 四种算法计算速度与精度的比较分析 时问比较时，采取计算２ １６０阶的缔合勒让德 函数所需时间，所得结果见表１。 表１ ４种递推法计算时间比较 单位：ｓ 方法 时间 标准前向列推法 ０．２８５ 标准前向行推法 ０．４２５ 跨阶次递推法 ０．１６１ Ｂｅｌｉｋｏｖ列推法 Ｏ．０４６ 进行精度评估 ： Ｊ（ ｃｏｓ －（２凡＋１）ｆ （７） 分别取 ＝３。、５４。进行编程计算，计算结果见 图１、２。 ×１０。 １．８ ——１．６ 标准向前列推法１ １ｔ４ ； ～一跨阶次递推法ＪＢｅ１ ０ｖ列推法１巨 一｛一…～…一 …一＿ ： 一． “ １．２ ｉ＿ １ 蜒Ｏ ８ Ｏ・６ 二 ：＿＿＿＿：一卜… 一 ０．４ ／ Ｉ ； ．． ０．２ ０ ：： ： 兰兰 ０ ５００ １０００ １５００ ２０００ ２５００ 阶数 图１ ＝３。时４种递推法的比较图示 从图１、２可以看出，随着纬度增大，各种递推算 法误差减小。当纬度较大时，标准向前行推法一旦 超过１４００阶，误差急剧增大；而跨阶次递推法和 Ｂｅｌｉｋｏｖ列推法在纬度较大时依然保持较好的精度， ・，一 ０ ５００ １０００ １５００ ２０００ ２５００ 阶数 图２ ＝５４。时４种递推法的比较图示 ３模型高程异常计算的优化算法 ３．１理论分析 假设（ｒ ， ，Ａ ）为第ｉ个纬度圈第 个网格中 点（曰 ， ， ）处的地心向径和经纬度，高程异常为 可表示为下式 ＇ｍ ： （Ｂ ， ，Ｈ ）＝∑（Ｅ ｉ ｃｏｓｍ ＋Ｆ ｓｉｎｍＡ，） 塞（ Ｃ—ｎ＊Ｐ ｓｉｎ （８） 差（ Ｓ （ｓｉ ） 式中，ｙ 为不同纬度下正常重力值，可根据索米里 安公式 求解。 由上面的分析知，由于顾及到网格平均高 ，， 同一纬度圈ｒ 提不等的，这样对每个点都必须重新 循环求解Ｐ （ｓｉｎｑ￣）、ｒ，运算量大。如果假定高程 Ｈ ＝０，此时高程异常为 （Ｂ ， ，０），式（８）简 化为： 靠（Ｂ ， ，０）＝∑（Ｅ ｉ ｃｏｓｍＡＪ＋Ｆ ｓｉｎｍＡ ） Ｎ Ｅ 。 ｎ ） （９） Ｎ ＝ ｎ 。） 利用式（９）计算 ，每一个纬度圈的缔合勒让 德函数只需计算一次，大大提高了运算速度。但计 算 时，由于没有考虑格网高程影响，显然对精度 影响很大。由分析可知， 是关于ｒｉｊ的函数，计算 实质就是将 在 ＝０处按泰勒公式展开并取至零 １６ 汉：武汉大学，２００７． 海洋测 绘 第３４卷 Ｐｒｅｅｉｓｅ Ｏｒｂｉｔｓ ｏｆ Ｌｏｗ Ｅａｒｔｈ Ｏｒｂｉｔｅｒｓ Ｕｓｉｎｇ ｔｈｅ ＧｌｏｂＮ ［２］ 唐卫明，刘 智．ＧＰＳ载波相位平滑伪距精度分析与 Ｐｏｓｉｔｉｏｎｇ Ｓｙｓｔｅｍ［Ｄ］．ＰＨ ｔｈｅｓｉｓ．Ａｓｔｒｏｎｏｍｉｃａｌ Ｉｎｓｔｉｔｕｔｅ Ｕｎｉｖｅｒｓｉｔｙ ｏｆ Ｂｅｒｎｅ，Ｓｅｈｗｉｔｚｅｒｌａｎｄ． 应用探讨［Ｊ］．测绘信息与工程，２００５，３０（３）：３７—３９． ［３］ 牛作鹏，张书毕，张丽．星历精度对ＧＰＳ单点定位的 涌．ＧＰＳ广播星历误差及其对导 ［８］ Ｓ Ｃ Ｈａｎ，Ｊ Ｈ Ｋｗｏｎ，Ｃ Ｊｅｋｅｌｉ．Ａｃｃｕｒａｔｅ ａｂｓｏｌｕｔｅ ＧＰＳ 影响【Ｊ］．海洋测绘，２０１０，３０（４）：５７—５９． ［４］ 帅平，陈定昌，江航定位精度的影响［Ｊ］．数据采集与处理，２００４，１９ （１）：１０７一ｌ１０． ｐｏｓｉｔｉｏｎｉｎｇ ｔｈｒｏｕｇｈ ｓａｔｅｌｌｉｔｅ ｃｌｏｃｋ ｅｒｒｏｒ ｅｓｔｉｍａｔｉｏｎ［Ｊ］． Ｊｏｕｒｎａｌ ｏｆ Ｇｅｏｄｅｓｙ，２００１（７５）：３３—４３． ［９］ 张全德．海岛（礁）测绘工程［Ｒ］．中国地理信息应用 报告，２０ｌ０． ［５］ 刘朋姣，陈义．基于ＩＧＵ星历的实时卫星钟差加密 ［１０］ 孙成星，许家琨．浙江南鼎星岛全野外数字测图关键 方法研究［Ｊ］．全球定位系统，２０１０（６）：３５—３９． 技术应用研究［Ｊ］．浙江测绘，２０１０（３）：２—１１． 鄂栋臣，詹必伟，姜卫平，等．应用ＧＡＭＩＴ／ＧＬ０ＢＫ软 件进行高精度ＧＰＳ数据处理［Ｊ］．极地研究，２００５，１７ （３）：１７３．１８２． ［６］ 叶世榕．ＧＰＳ非差相位精密度单点定位理论与实现 ［Ｄ］．武汉：武汉大学，２００２． ［７］ Ｂｏｅｋ Ｈ（２００３）．Ｅｆｆｉｅｉｅｎｔ Ｍｅｔｈｏｄｓ ｆｏｒ Ｄｅｔｅｒｍｉｎｉｎｇ Ｒｅｓｅａｒｃｈ ｏｎ Ｃａｌｃｕｌａｔｉｎｇ Ｓｃｈｅｍｅ Ｏｐｔｉｍｉｚａｔｉｏｎ ｏｆ Ｐｒｅｃｉｓｅ Ｐｏｉｎｔ Ｐｏｓｉｔｉｏｎｉｎｇ ｉｎ Ｓｅａｃｏａｓｔ Ｃｏｎｔｒｏｌ Ｓｕｒｖｅｙｓ ＺＨＡＮＧ Ｈｕａ，ＴＡＮＧ Ｊｕｎ，ＳＯＮＧ Ｇｕｏｄａ，ＣＡＯ Ｊｕｎｈｏｎｇ （Ｎａｖｙ Ｐｒｅｓｓ，Ｔｉａｎｊｉｎ ３００４５０，Ｃｈｉｎａ） Ａｂｓｔｒａｃｔ：Ａｉｍｉｎｇ ａｔ ｔｈｅ ｄｒａｗｂａｃｋｓ ｉｎ ｓｅａｃｏａｓｔ ｃｏｎｔｒｏｌ ｓｕｒｖｅｙｓ，ｔｈｉｓ ｐａｐｅｒ ｓｔｕｄｉｅｓ ｓｅｖｅｒａｌ ｏｐｔｉｍｉｚｉｎｇ ｃａｌｃｕｌａｔｉｎｇ ｍｅｔｈｏｄｓ ｏｆ ｐｒｅｃｉｓｅ ｐｏｉｎｔ ｐｏｓｉｔｉｏｎｉｎｇ ｗｈｉｃｈ ｉｎｃｌｕｄｅ ｄｅｌｅｔｉｎｇ ｏｂｓｅｒｖａｔｉｏｎ ｔｉｍｅ ｉｎ ｐｏｏｒ ｃｏｎｄｉｔｉｏｎｓ，ｓｈｏｒｔｅｎｉｎｇ ｏｂｓｅｒｖａｔｉｏｎ ｔｉｍｅ，ｓｕｂｓｔｉｔｕｔｉｎｇ ｐｒｅｃｉｓｅ ｅｐｈｅｍｅｒｉｓ ｗｉｈ ｑｕｉｔｃｋ ｅｐｈｅｍｅｒｉｓ ａｎｄ ｕｓｉｎｇ ｑｕｉｃｋ ｅｐｈｅｍｅｒｉｓ ｃａｌｃｕｌａｔｉｏｎ ｔｏ ｓｈｏ￣ｅｎ ｏｂｓｅｒｖａｔｉｏｎ ｔｉｍｅ，ａｎｄ ｐｒｏｐｏｓｅｓ ｔｗｏ ｏｐｔｉｍｉｚｅｄ ｃａｌｃｕｌａｔｉｎｇ ｓｃｈｅｍｅｓ ｂａｓｅｄ ｏｎ ａｎａｌｙｚｉｎｇ ｔｈｅ ａｄｖａｎｔａｇｅｓ ａｎｄ ｄｉｓａｄｖａｎｔａｇｅｓ ｏｆ ｔｈｅ ａｂｏｖｅ ｍｅｔｈｏｄｓ．Ｅｘｐｅｒｉｍｅｎｔａｌ ｒｅｓｕｌｔ ｓｈｏｗｓ ｔｈａｔ ｔｈｅ ｐｒｏｐｏｓｅｄ ｓｃｈｅｍｅｓ ｃａｎ ｒｅｄｕｃｅ ｔｈｅ ｗｏｒｋｅｒ’ｓ ｌｏａｄ ａｎｄ ｍｅｅｔ ｔｈｅ ｎｅｅｄｓ ｏｆ ｓｅａｃｏａｓｔ ｃｏｎｔｒｏｌ ｓｕｒｖｅｙｓ． Ｋｅｙ ｗｏｒｄｓ：ｓｅａｃｏａｓｔ ｃｏｎｔｒｏｌ ｓｕｒｖｅｙ；ｐｒｅｃｉｓｅ ｐｏｉｎｔ ｐｏｓｉｔｉｏｎｉｎｇ（ＰＰＰ）；ｑｕｉｃｋ ｅｐｈｅｍｅｒｉｓ ◆…Ｉ◆ｌｌＩ◆ｉｊ◆… ●…◆…ｌ◆ｌｌ◆…Ｉ◆ｌｌ｛◆…ｌ◆…◆…Ｉ◆ｌｌ】◆…Ｉ◆ｌ ｌｊ◆…Ｉ◆ｌｌ】◆…Ｉ４．Ｐ…◆…Ｉ◆ｌｌ●…ｌ◆ｌｌ◆¨ｌ◆ｌｌ◆¨ｌ◆ｌ ４．ｌ１¨ｌ◆ｌ ４．１¨Ｉｌ◆…◆ｉｌｌ◆…◆¨ｌ◆¨◆１¨ｌ◆¨◆¨ｌ◆ｉｌｌ◆…◆…ｂｏ．１ｌ◆…◆…◆…◆ （上接第８页） Ｏｐｔｉｍｉｚａｔｉｏｎ Ａｌｇｏｒｉｔｈｍ ｆｏｒ Ｈｅｉｇｈｔ Ａｎｏｍａｌｙ ｏｆ Ｕｌｔｒａ—ｈｉｇｈ Ｄｅｇｒｅｅ ａｎｄ Ｏｒｄｅｒ Ｅａｒｔｈ Ｇｒａｖｉｔｙ Ｆｉｅｌｄ Ｍｏｄｅｌ Ｑｕ Ｚｈｅｎｇｈａｏ，ＺＥＮＧ Ｔｉａｎ，ＺＨＯＵ Ｑｉａｎｇ （Ｉｎｓｔｉｔｕｔｅ ｏｆ Ｇｅｏ－ｓｐａｔｉａｌ Ｉｎｆｏｒｍａｔｉｏｎ，Ｉｎｆｏｒｍａｔｉｏｎ Ｅｎｇｉｎｅｅｒｉｎｇ Ｕｎｉｖｅｒｓｉｔｙ，Ｚｈｅｎｇｚｈｏｕ ４５００５２，Ｃｈｉｎａ） Ａｂｓｔｒａｃｔ：Ｆｏｕｒ ｒｅｃｕｒｓｉｏｎ ｍｅｔｈｏｄｓ ｏｆ ｃｏｍｐｕｔｉｎｇ ｔｈｅ ｕｌｔｒａ—ｈｉｇｈ ｄｅｇｒｅｅ ａｎｄ ｏｒｄｅｒ ｆｕｌｌｙ ｎｏｒｍａｌｉｚｅｄ ａｓｓｏｃｉａｔｅｄ Ｌｅｇｅｎｄｒｅ ｆｕｎｃｔｉｏｎｓ ｓｕｃｈ ａｓ ｓｔａｎｄａｒｄ ｆｏｒｗａｒｄ ｃｏｌｕｍｎ ｍｅｔｈｏｄ，ｓｔａｎｄａｒｄ ｆｏｒｗａｒｄ ｒｏｗ ｍｅｔｈｏｄ，ｒｅｃｕｒｓｉｖｅ ｍｅｔｈｏｄ ｂｅｔｗｅｅｎ ｅｖｅｒｙ ｏｔｈｅｒ ｏｒｄｅｒ ａｎｄ ｄｅｇｒｅｅ ａｎｄ Ｂｅｌｉｋｏｖ ｃｏｌｕｍｎ ｍｅｔｈｏｄ，ｗｅｒｅ ａｎａｌｙｚｅｄ ａｎｄ ｃｏｍｐａｒｅｄ ｔｏ ｔｈｅ ａｄｖａｎｔａｇｅｓ ａｎｄ ｄｉｓａｄｖａｎｔａｇｅｓ ｆｒｏｍ ｔｈｒｅｅ ａｓｐｅｃｔｓ ｏｆ ｎｕｍｅｒｉｃａｌ ｐｒｅｃｉｓｉｏｎ，ｎｕｍｅｒｉｃａｌ ｓｔａｂｉｌｉｔｙ ａｎｄ ｎｕｍｅｒｉｃａｌ ｅｆｉｃｉｅｎｃｙ ｉｎ ｔｆｈｉｓ ｐａｐｅｒ．Ｔｈｅ Ｂｅｌｉｋｏｖ ｃｏｌｕｍｎ ｍｅｔｈｏｄ ｗａｓ ｃｈｏｓｅｎ ｔｏ ｃａｌｃｕｌａｔｅ ｔｈｅ ｈｅｉｇｈｔ ａｎｏｍａｌｙ ｏｆ ｔｈｅ ｅａｒｔｈ’ｓ ｇｒａｖｉｔｙ ｆｉｅｌｄ ｍｏｄｅ１．Ｔｈｅｎ ｔｈｒｅｅ ｄｉｆｆｅｒｅｎｔ ａｌｇｏｒｉｔｈｍｓ ｆｏｒ ｃａｌｃｕｌａｔｉｎｇ ｍｏｄｅｌ ｈｅｉｇｈｔ ａｎｏｍａｌｙ ｂａｓｅｄ ｏｎ ｓｔｒｉｃｔ ｓｐｈｅｒｉｃａｌ ｈａｒｍｏｎｉｃ ｅｘｐａｎｓｉｏｎ ｆｏｒｍｕｌａ，ｒｅｔａｉｎｉｎｇ ｔｈｅ ｆｏｒｍｕｌａ’ｓ ｆｕｎｃｔｉｏｎ ｏｆ Ｔａｙｌｏｒ ｓｅｒｉｅｓ ａｂｏｕｔ ｇｅｏｃｅｎｔｒｉｃ ｒａｄｉｕｓ ｖｅｃｔｏｒ ｕｎｆｏｌｄｅｄ ｔｏ ｚｅｒｏ ｏｒｄｅｒ ｔｅｒｍ ａｎｄ ｒｅｔａｉｎｉｎｇ ｔｈｅ ｆｏｒｍｕｌａ’ｓ ｆｕｎｃｔｉｏｎ ｏｆ Ｔａｙｌｏｒ ｓｅｒｉｅｓ ｕｎｆｏｌｄｅｄ ｔｏ ｏｎｅ ｏｒｄｅｒ ｔｅｒｍ ｗｅｒｅ ｄｉｓｃｕｓｓｅｄ，ａｎｄ ｓｏｍｅ ｅｘｐｅｒｉｍｅｎｔｓ ｗｅｒｅ ｃｏｎｄｕｃｔｅｄ．Ｔｈｅ ｒｅｓｕｌｔ ｓｈｏｗｅｄ ｔｈａｔ ｕｎｆｏｌｄｉｎｇ ｔｏ ｆｉｒｓｔ ｏｒｄｅｒ Ｔａｙｌｏｒ ｓｅｒｉｅｓ ｏｆ ｆｏｒｍｕｌａ ａｂｏｕｔ ｇｅｏｃｅｎｔｒｉｃ ｒａｄｉｕｓ ｖｅｃｔｏｒ ｔｏ ｃａｌｃｕｌａｔｅ ｈｅｉｇｈｔ ａｎｏｍａｌｙ ｏｆ ｔｈｅ ｅａｒｔｈ’ｓ ｇｒａｖｉｔｙ ｆｉｅｌｄ ｍｏｄｅｌ ｃｏｕｌｄ ｂｏｔｈ ｅｎｓｕｒｅ ｇｏｏｄ ｃｏｍｐｕｔｉｎｇ ｓｐｅｅｄ ａｎｄ ｒｅａｃｈ ｅｎｏｕｇｈ ｐｒｅｃｉｓｉｏｎ，ｍｅｅｔｉｎｇ ｔｈｅ ｒｅｑｕｉｒｅｍｅｎｔ ｏｆ ｂｕｉｌｄｉｎｇ ｈｅ ｌａｔｒｇｅ ａｒｅａ ｈｉｇｈ—ｒｅｓｏｌｕｔｉｏｎ ｈｅｉｇｈｔ ａｎｏｍａｌｙ ｍｏｄｅ１． Ｋｅｙ ｗｏｒｄｓ：ａｓｓｏｃｉａｔｅｄ Ｌｅｇｅｎｄｒｅ ｆｕｎｃｔｉｏｎ；ｕｌｔｒａ—ｈｉｇｈ ｄｅｇｒｅｅ ａｎｄ ｏｒｄｅｒ ｅａｒｔｈ ｇｒａｖｉｔｙ ｆｉｅｌｄ ｍｏｄｅｌ；ｈｅｉｇｈｔ ａｎｏｍａｌｙ；Ｔａｙｌｏｒ ｓｅｒｉｅｓ ｅｘｐａｎｓｉｏｎ
