-
受地球自转和轨道摄动的影响, 卫星的星下点轨迹在地球表面会形成交叉点。利用卫星测高交叉点处海面高不符值进行平差计算可以有效削弱卫星径向轨道误差、海面时变误差以及系统误差等,提高海面高数据的精度。因而交叉点数据在卫星大地测量和海洋学领域得到广泛应用[1, 2],例如,测高仪时标偏差校准[3]、重力场模型精度评估[4, 5]、冰盖高程变化监测[6, 7]、海洋重力场反演[8]、大洋环流[9, 10]、海洋潮汐[11, 12]等研究。
由于卫星轨道的离散采样,交叉点的位置需要通过内插计算得到[2]。对于单星交叉点的计算相对简单,可根据卫星轨道六参数和地球自转建立模型 (这里只考虑卫星地球二体运动并且地球视为匀质椭球体) 计算出理论交叉点[13-15],再将其作为初值内插计算精确交叉点[6]。对于双星交叉点而言,其理论交叉点的计算较为复杂,利用数值法直接计算交叉点位置较为方便。即对卫星地面轨迹进行曲线拟合,再通过计算曲线交点来确定交叉点,目前用的最多是二次多项式拟合[8, 15-16]。二项式拟合法的计算效率取决于初值的选取,并需求解二次方程组,存在解算失败导致交叉点缺失的可能, 在高纬度地区这一现象尤为严重。一些学者在二项式拟合法基础上提出了分区拟合或分段拟合的改进方法,但是并没有解决本质问题,同时还引入了如何最优分区或分段的新问题[17, 18]。为了克服上述不足,本文提出了一种快速定位交叉点的矩形收缩法,能够高效准确地搜索到每条轨迹上距离交叉点最近的数据采样点,最后采用求两直线交点的方法计算得到精确的交叉点位置。
-
数据预处理包括经度跳变的处理和交叉点判断两项内容。卫星运行有顺行和逆行之分,顺行卫星轨道倾角小于90°,经度随时间增加;逆行卫星轨道倾角大于90°,经度随时间减小。在卫星测高中,地面轨迹根据星下点的纬度变化分为升弧段和降弧段。因此,理论上星下点的经纬度都是连续单调变化的。但是由于经度的取值范围是[0°, 360°), 如果卫星轨迹穿过本初子午线,其星下点经度在0°(360°) 附近会产生跳变。以前处理经度跳变是通过加或减360°的方法使经度变化连续,这一做法的不便之处是需要根据求交叉点的另一条轨迹情况来决定加或减。此外,当两条轨迹存在两个交叉点时,这一方法可能导致只能求得一个交叉点。本文直接将跨越本初子午线的卫星轨迹分成两个弧段。
为了避免无谓的计算,事先判断两条弧段是否存在交叉点是非常有必要的。对于同一个卫星而言,只有升弧段和降弧段之间可能存在交叉点。对于两个倾角不同的卫星,升 (降) 弧段与升 (降) 弧段之间也可能存在交叉点。现有文献中通常只讨论了升-降弧段之间的交叉点存在条件[18],而很少述及同升或同降弧段之间交叉点的存在条件。若顾及轨道运行方向和弧段的升降,两个相交弧段存在10种可能的组合 (图 1),其中 (1) 和 (2) 两种形式即是单星交叉点的组合。这些组合可以归纳为“同升同降”(组合 (1)~(4)) 和“一升一降”(组合 (5)~(10)) 两类,它们存在交叉点时必须满足以下条件:
(1) 式中,λFirst1、λLast1、λFirst2、λLast2分别表示两条弧段第一点和最后一点的经度。上述条件只是存在交叉点的必要而非充分条件,使用时还应使两条弧段的纬度范围保持一致。
-
矩形收缩法的主要思想是先根据两条弧段的端点坐标确定一个初始矩形,然后按照一定的规则不断缩小矩形,直到其中一条弧段在矩形范围内的数据点小于两个时停止收缩,最后利用矩形内的数据点内插出交叉点概略位置。其具体步骤如下。
1) 设两条弧段4个端点的坐标为 (λ1, φ1)、(λ2, φ2)、(λ3, φ3)、(λ4, φ4),分别取4个端点经度和纬度的两个中值为初始矩形的经度、纬度范围 (λ0L~λ0H, φ0L~φ0H)。矩形的两条经线分别与两条弧段各有一个交点,对应的纬度差为ΔφL、ΔφH,利用这两个纬度差可以计算出一个收缩系数s:
(2) 2) 根据收缩系数缩小矩形的经度范围,收缩的规则为:
(3) (4) 式中,L0=c×(λ0H-λ0L),c为 (0, 1) 区间内的任意常数。理论上c值越大,收缩速度越快,但是由于数据是离散采样的,c值过于接近1可能导致无法进行迭代收缩。试验表明,c取0.8~0.95时收敛速度较快,且迭代次数仅相差1~2次,故在本文算例中取为0.9。
3) 由上一步的结果可以截取两条新的子弧段,重复前两个步骤,直到其中一条子弧段的点数小于2为止。
4) 对两条子弧段进行直线拟合,计算两直线的交点作为概略交叉点。
图 2给出了矩形收缩法的示意图,实线矩形给出的是同一卫星两条弧段 (点划线) 的概略交叉点搜索情形,虚线矩形是不同倾角的不同卫星两条弧段的情况。可以看出,同一卫星的两条弧段相对其交叉点具有对称性,收缩系数为0.5。对于不同卫星的两条弧段,收缩系数则与两颗卫星的轨道倾角相关,同时也会随交叉点的纬度变化而变化。不管哪种情况,矩形收缩一般只需2~3次迭代就可完成。
通过以上迭代过程得到的概略交叉点已经非常接近真实交叉点位置。为了进一步得到精确的交叉点位置,只需分别在两条弧段上按经度大小查找出离概略交叉点最近的两个数据点,并进行线性内插即可。理论上,交叉点应位于最近两点的中间,可以根据这一条件来迭代计算交叉点的精确位置[18]。但是实测数据可能存在因数据缺失而导致距离最近的两点位于交叉点同一侧的情形,这种情况下如果采用迭代计算则无法收敛。大量试算表明,利用矩形收缩法得到的概略交叉点求精确交叉点仅需一次计算就能完成,无需迭代计算。如果最后查找到的实测数据点离概略交叉点很远时,内插出的交叉点数据意义不大,一般不再计算该交叉点。
-
为了验证算法的有效性,本文进行了3个算例的计算。
算例1利用T/P卫星的相关参数模拟了一个周期的卫星轨道数据,并计算其交叉点。之所以采用模拟数据,是因为可以得到交叉点的理论值,便于对算法的精度进行评估。在球形地球的假设下,圆轨道卫星的地面轨迹可以表示为时间的函数[14]:
(5) (6) 式中,φ、λ分别是轨迹点的纬度和经度;i是轨道倾角,T/P卫星为66°;Ω是相对轨道面的地球自转角速度;λ0为地面轨迹与赤道交点的经度值,这里采用T/P卫星GDR (geophysical data record) 数据手册中给出的数值[19];Δt是相对过赤道时刻的时间差,满足|Δt|≤π/2;j取1为升弧段,-1为降弧段。图 3中的灰色曲线即为利用上述方法模拟的T/P卫星轨道,升弧段和降弧段各127条。与GDR数据一样,每条弧段有3 360个采样点。它们的交叉点可以通过解析式精确计算得到[14], 共有16 002个理论交叉点。
图 3 T/P卫星模拟轨道 (灰色) 和HY-2卫星第91号轨道 (黑色)
Figure 3. Simulated T/P Tracks (in Gray) and No. 91 Track of HY-2(in Black)
算例2选取了海洋二号 (HY-2) 卫星IGDR (interim geophysical data record) 数据中的一条轨道,图 3中黑色曲线为第14周期第91号升弧段,计算其与模拟T/P数据的交叉点。HY-2的轨道倾角[20]约99.34°,为逆行轨道,所以算例2可以评估算法用于计算两颗不同倾角卫星之间交叉点的性能,同时也能检验算法处理实际卫星轨道数据的能力。
算例3选用Envisat卫星第90周期的实测数据,目的是测试算法在计算近极轨道在两极附近的交叉点的适用性。
-
利用矩形收缩法、二项式拟合法和分段拟合法分别计算了两个算例的交叉点,其中分段拟合法采用了文献[10]给出的针对Jason-1卫星最优的type-12分段方式。三种方法均采用距离交叉点最近的两点来计算精确交叉点,所以得到的交叉点精度是一致的。交叉点计算精度利用算例1的结果与理论交叉点进行比较得出,这里仅给出矩形收缩法计算结果的精度 (表 1)。可以看出,低纬度的交叉点计算精度较高,纬度越大交叉点的精度越差。对于模拟T/P轨道,交叉点经度的计算精度已经达到机器误差水平 (计算过程中数据只保存到小数点后10位),纬度误差稍大,高纬地区的交叉点精度达0.000 01°,足以满足各种应用对交叉点的精度要求。
表 1 矩形收缩法的交叉点计算精度
Table 1. Accuracy of Crossovers Computed by Rectangle Shrinking Algorithm
纬度范围 |φ| < 30° 30°≤|φ| < 60° |φ|≥60° 全球 交叉点个数 2 540 6 604 6 858 16 002 经度标准差/(°) 0 7.0×10-12 1.2×10-10 8.0×10-11 纬度标准差/(°) 2.4×10-6 8.1×10-6 1.1×10-5 8.9×10-6 尽管3种方法都能得到同样精度的结果,但是在效率上有着明显差别。效率包括所用的时间和计算得到的交叉点数量两个方面,表 2给出了3种方法的表现。3种交叉点计算方法的程序均由Visual Fortran 6.6软件编写和编译,并在电脑 (CPU主频2.66 GHz、内存5 GB,64位Windows 7操作系统) 上运行。表 2中统计的计算时间为各程序运行10次的平均值。
表 2 3种计算方法的效率比较
Table 2. Comparison of the Computational Efficiency of Three Methods
本文方法 二项式拟合法 分段拟合法 算例1 交叉点个数 16 002(100%) 16 002(100%) 15 748(98.4%) 计算时间/s 34.854 35.265 35.403 算例2 交叉点个数 127(100%) 113(89%) 81(63.8%) 计算时间/s 0.522 0.554 0.560 算例1中,本文提出的矩形收缩法计算时间最短,并计算出了全部交叉点。而分段拟合法用的时间最长,且丢失了254个交叉点,约占总数的1.6%。二项式拟合法的效率取决于初始选取的拟合数据,若采用整条弧段直接进行拟合不仅计算时间长,还会导致大量交叉点无法计算得到[10]。本文以两弧段4个端点经度的平均值为参考,选取经度与平均值差值最小的8个数据点作为初始拟合数据,极大改善了二项式拟合法的计算效率,计算得出了全部交叉点,计算时间位于矩形收缩法和分段拟合法之间。
对分段拟合法计算结果缺失的交叉点进行分析发现,这些交叉点分布在南北纬60°附近。显然,这是由于对弧段进行人为分段引起的。以图 4所示的第1号升弧和第32号降弧为例,其交叉点位于北纬60°附近。由于分段拟合时恰好以60°纬线为分界将弧段分开,得到的子弧段不满足存在交叉点的第2个条件,从而被误判为不存在交叉点。因此,在采用分段拟合或类似分区方法时,应根据交叉点分布合理选择分段或分区方式,分段节点不能在交叉点纬度带附近。在交叉点分布未知时,可使相邻两个分段或分区之间保留一定的重合区域,以避免类似的误判结果。
图 4 T/P轨迹第1号升弧段和第32号降弧段的交叉点
Figure 4. Crossover of No. 1 Ascending Pass and No. 32 Descending Pass of T/P
算例2的结果与算例1类似,矩形收缩法以最快的速度计算出了全部交叉点,二项式拟合法解算出89%的交叉点,而分段拟合法的解算成功率仅有63.8%。图 5给出了3种方法计算出的交叉点分布情况,右上角的子图是北纬64°以上高纬度地区的交叉点解算结果。图 5中小圆点、空心圆和竖线段分别是矩形收缩法、二项式拟合法以及分段拟合法计算的交叉点,虚线为分段拟合法采用的分段节点纬线。结果显示,二项式拟合法未能成功解算出的14个交叉点均位于大于64°的高纬度地区,是由于拟合的二项式联立求解失败造成的。分段拟合法除了分段节点附近的交叉点缺失外,在高纬度地区仍有大量交叉点不能成功解算,表明该方法并不能解决二项式拟合法的缺点,在实测卫星数据点缺失较严重的情况下反而会增加解算不成功的可能。
算例3采用本文方法计算了Envisat卫星轨道在南极附近的交叉点,图 6给出了部分区域的结果。可以看出,即使在接近极点的区域,本文方法仍能成功计算得到高精度的交叉点,验证了算法对近极轨道的适用性。
-
卫星测高技术发展早期通常只处理单一卫星数据,因此交叉点被认为只是升弧段与降弧段的交点。随着测高卫星的增多,多种卫星数据的联合处理成为一种趋势。对于不同倾角的两个卫星而言,升弧段与升弧段、降弧段与降弧段也可能存在交叉点。为了提高交叉点计算方法的通用性和效率,本文扩展了交叉点存在的判断条件,可用于判断任意两条弧段是否有交叉点。本文应用纯几何原理,提出了一种快速搜索概略交叉点进而计算精确交叉点的矩形收缩算法。该算法既不需要卫星轨道的设计参数,也不要进行复杂的多项式拟合及二次方程组求解的运算。两个算例的结果表明,本文提出的矩形收缩法不仅计算速度快、精度高,而且没有交叉点遗漏现象。对于单一卫星的交叉点计算,若能给定合适的初值,二项式拟合法不失为一种比较好的方法。当计算不同卫星的轨道交叉点时,由于初值难于确定,所以不宜使用二项式拟合法。同时,不建议采用分段拟合法计算交叉点,因为其本质仍是二项式拟合,不能完全克服二项式拟合法的缺点,计算速度也没有提高, 反而可能因为分段造成更多数量的交叉点丢失。分段拟合法的另一个不足是不同卫星需要选择不同的最优分段类型,即最优分段类型不具有通用性。
Fast Numerical Algorithm for the Calculation of Altimetric Crossovers from Satellite Ground Tracks
-
摘要: 计算交叉点是卫星测高数据处理中的重要基础性工作。扩展了交叉点存在的判断条件,可用于判断任意两条卫星地面轨迹是否有交叉点。提出了一种快速计算交叉点的数值算法--矩形收缩算法。采用一个周期的Topex/Poseidon(T/P)卫星模拟轨道和一条海洋二号(HY-2)卫星实际轨迹设计了两个算例,以验证算法的精度和效率。结果表明矩形收缩法可以快速、高精度地计算出全部交叉点。以Envisat数据为例验证了算法计算近极轨道两极交叉点的适用性。该方法不仅可以计算单一卫星轨迹的交叉点,也可计算两个不同倾角卫星的轨迹交叉点,具有很强的通用性。Abstract: Numerical calculation of crossovers is important groundwork for satellite altimeter data processing. The criteria for judging the existence of crossover are extended to be applicable to any two satellite ground tracks. A fast numerical algorithm named rectangle shrinking method is proposed in order to improve the computational efficiency. Based on a cycle of simulated orbit data of Topex/Poseidon (T/P) and a pass of observations from HY-2, two experiments are performed to assess the precision and efficiency of the algorithm. The results demonstrate that the rectangle shrinking method can rapidly work out all of crossovers with high precision. Another test on Envisat data validated the applicability to near polarorbit. The proposed algorithm has strong universality, not only solving the crossovers from tracks of a single satellite, but also two satellites with different inclinations.
-
Key words:
- satellite altimetry /
- crossover /
- satellite ground track /
- fast numerical algorithm
-
表 1 矩形收缩法的交叉点计算精度
Table 1. Accuracy of Crossovers Computed by Rectangle Shrinking Algorithm
纬度范围 |φ| < 30° 30°≤|φ| < 60° |φ|≥60° 全球 交叉点个数 2 540 6 604 6 858 16 002 经度标准差/(°) 0 7.0×10-12 1.2×10-10 8.0×10-11 纬度标准差/(°) 2.4×10-6 8.1×10-6 1.1×10-5 8.9×10-6 表 2 3种计算方法的效率比较
Table 2. Comparison of the Computational Efficiency of Three Methods
本文方法 二项式拟合法 分段拟合法 算例1 交叉点个数 16 002(100%) 16 002(100%) 15 748(98.4%) 计算时间/s 34.854 35.265 35.403 算例2 交叉点个数 127(100%) 113(89%) 81(63.8%) 计算时间/s 0.522 0.554 0.560 -
[1] Kim M C.Theory of Satellite Ground-track Crossovers[J].Journal of Geodesy, 1997, 71:749-767 doi: 10.1007/s001900050141 [2] 翟国君, 黄谟涛, 谢锡君, 等.卫星测高数据处理的理论与方法M].北京:测绘出版社, 2000 Zhai Guojun, Huang Motao, Xie Xijun, et al. The Theory and Method of Satellite Altimeter Data Processing[M].Beijing:Surveying and Mapping Press, 2000 [3] Schut B E, Tapley B D, Shum C K. Evaluation of the SEASAT Altimeter Time Tag Bias[J]. Journal of Geophysical Research, 1982, 87(C5):3239-3245 doi: 10.1029/JC087iC05p03239 [4] Klokoník J, Wagner C A, Kostelecky J. Spectral Accuracy of JGM-3 from Satellite Crossover Altimetry[J]. Journal of Geodesy, 1999, 73:38-46 doi: 10.1007%2Fs001900050229 [5] Klokoník J, Wagner C A, Kostelecky J, et al. Testing the Test: Accuracy Assessment of Gravity Field Models by Independent Satellite Crossover Altimetry, with Example for EIGENIS[C]. Satellite Altimetry for Geodesy, Geophysics and Oceanography, International Association of Geodesy Symposia, Nice, France, 2004 [6] 史红岭, 陆洋, 鲍李峰, 等.利用ICESat交叉点分析探测恩德比地冰盖近年高程变化[J].武汉大学学报·信息科学版, 2009, 34(4):440-443 http://ch.whu.edu.cn/CN/Y2009/V34/I4/440 Shi Hongling, Lu Yang, Bao Lifeng, et al. Recent Elevation Change Detection of Enderby Land Ice Sheet Using ICESat Crossover Analysis[J].Geomatics and Information Science of Wuhan University, 2009, 34(4):440-443 http://ch.whu.edu.cn/CN/Y2009/V34/I4/440 [7] 杨元德, 鄂栋臣, 王泽明, 等.利用Envisat数据探测中山站至Dome A条带区冰盖高程变化[J].武汉大学学报·信息科学版, 2013, 38(4):383-385 http://ch.whu.edu.cn/CN/Y2013/V38/I4/383 Yang Yuande, E Dongchen, Wang Zemin, et al. Elevation Change from Zhongshan Station to Dome A Using Envisat Data[J]. Geomatics and Information Science of Wuhan University, 2013, 38(4):383-385 http://ch.whu.edu.cn/CN/Y2013/V38/I4/383 [8] 黄谟涛, 王瑞, 翟国君, 等.多代卫星测高数据联合平差及重力场反演[J].武汉大学学报·信息科学版, 2007, 32(11):988-993 http://ch.whu.edu.cn/CN/abstract/abstract2021.shtml Huang Motao, Wang Rui, Zhai Guojun, et al. Integrated Data Processing for Multi-Satellite Missions and Recovery of Marine Gravity Field[J]. Geomatics and Information Science of Wuhan University, 2007, 32(11):988-993 http://ch.whu.edu.cn/CN/abstract/abstract2021.shtml [9] Hwang W, Shih C, Guo Jinyun, et al. Zonal and Meridional Ocean Currents at Topex/Poseidon and JASON-1 Crossovers Around Taiwan: Error Analysis and Limitation[J]. Terrestrial, Atmospheric and Oceanic Sciences, 2008, 19(1/2): 151-162 https://www.researchgate.net/profile/Cheinway_Hwang/publication/250211930_Zonal_and_Meridional_Ocean_Currents_at_TOPEXPoseidon_and_JASON1_Crossovers_around_Taiwan_Error_Analysis_and_Limitation/links/56a9874c08aef6e05df2b966.pdf?origin=publication_detail [10] 郭金运, 常晓涛, 黄金维, 等.基于卫星测高交叉点的海洋表面地转流速度[J].地球物理学报, 2010, 53(11): 2582-2589 http://www.cnki.com.cn/Article/CJFDTOTAL-DQWX201011009.htm Guo Jinyun, Chang Xiaotao, Hwang C W, et al. Oceanic Surface Geostrophic Velocities Determined with Satellite Altimetric Crossover Method[J].Chinese Journal Geophysics, 2010, 53(11): 2582-2589 http://www.cnki.com.cn/Article/CJFDTOTAL-DQWX201011009.htm [11] Wagner C A, Tai C K, Kuhn J M. Improved M2 Ocean Tide from Topex/Poseidon Altimetry[J]. Journal of Geophysical Research, 1994, 99(C12):24853-24865 doi: 10.1029/94JC01347 [12] 许军, 暴景阳, 刘雁春, 等.以相邻点海面高度差为观测量的沿迹调和分析新方法[J].武汉大学学报·信息科学版, 2006, 31(11):1003-1006 http://ch.whu.edu.cn/CN/Y2006/V31/I11/1003 Xu Jun, Bao Jingyang, Liu Yanchun, et al. Harmonic Tidal Analysis Along T/P Tracks Taking Difference of Observed Sea Surface Heights at Adjacent Points as Observations[J]. Geomatics and Information Science of Wuhan University, 2006, 31(11):1003-1006 http://ch.whu.edu.cn/CN/Y2006/V31/I11/1003 [13] Hagar H. Engineering Memorandum 315-29: Computation of Circular Orbit Ground Trace Intersections[R]. Jet Propulsion Laboratory, Pasadena, CA, 1977 [14] Tai C K. An Efficient Algorithm for Computing the Crossovers in Satellite Altimetry[J].Journal of Atmospheric and Oceanic Techmology, 1989, 6:13-18 doi: 10.1175/1520-0426(1989)006<0013:AEAFCT>2.0.CO;2 [15] 王广运, 王海瑛, 许国昌.卫星测高原理[M].北京:科学出版社, 1995:207-208 Wang Guangyun, Wang Haiying, Xu Guochang. The Theory of Satellite Altimetry[M]. Beijing:Science Press, 1995:207-208 [16] 姜卫平.卫星测高技术在大地测量学中的应用[D].武汉:武汉大学, 2001:25-27 Jiang Weiping. The Application of Satellite Altimetry in Geodesy[D]. Wuhan:Wuhan University, 2001:25-27 [17] Spoecker N T. Crossover Generation[R]. Technical Report ERS-D-ADP-32601, DGFI Dept. 1, Munich, Germany, 1991 [18] 周晓光, 苗洪利, 王云海, 等.卫星地面轨迹分段拟合确定交叉点的方法研究[J].测绘学报, 2012, 41(6):811-815 http://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201206006.htm Zhou Xiaoguang, Miao Hongli, Wang Yunhai, et al. Study on the Determination of Crossovers by Piecewise Fitting of Satellite Ground Track[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(6):811-815. http://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201206006.htm [19] AVISO/Altimetry. AVISO User Handbook for Merged Topex/Poseidon Products[R]. AVI-NT-02-101, Edition 3.0, 1996 [20] 赵罡, 周旭华, 吴斌.海洋二号卫星SLR精密定轨[J].科学通报, 2012, 57(36):3475-3483 http://www.cnki.com.cn/Article/CJFDTOTAL-KXTB201236006.htm Zhao Gang, Zhou Xuhua, Wu Bin. Precise Orbit Determination of Haiyang-2 Using Satellite Laser Ranging[J]. Chin Sci Bull, 2012, 57(36):3475-3483 http://www.cnki.com.cn/Article/CJFDTOTAL-KXTB201236006.htm -