2. 中国科学院遥感应用研究所,北京 100101
2. Institute of Remote Sensing Applications, Chinese Academy of Sciences, Beijing 100101, China
1 引 言
月球测绘是月球探测工程和月球研究中的基础工作。中国第一颗月球探测卫星嫦娥一号(CE-1)获取了月面CCD影像和激光高度计(laser altimeter,LAM)数据,为研究月球地质特征、月球起源和演变提供基础的地形地貌信息。其中,激光高度计测量卫星到月表星下点的距离,用于制作月球全球表面数字高程模型(DEM)。文献[1]对CE-1两个月正飞阶段的激光观测数据进行了转换和处理,获取了第一个月球全球地形模型CLTM-s01,文献[2]制作了空间分辨率3 km的全月DEM模型,文献[3]也利用CE-1卫星激光测高数据进行了月球DEM的制作,并与Clementine、ULCN2005、CLTM-s01进行了比较。在采用CE-1激光测高数据生成DEM时,局部区域尤其在比较平坦区域的DEM存在条带噪声,表明由于轨道与测量的误差,多个轨道的高度计数据间以及与CE-1获取的CCD影像之间存在着误差和不一致性[4, 5]。为获取更高精度的月球DEM,为其他科学研究提供地形数据支持,对月球激光高度计数据进行误差分析和平差处理研究是非常有必要的。
在国内外对地观测海平面模型和重力场模型研究中,交叉点分析常用于消除轨道的径向误差[6, 7, 8, 9, 10, 11, 12]。在美国宇航局的火星MOLA(Mars orbiter laser altimeter)数据、月球LOLA(lunar orbiter laser altimeter)数据的处理中,交叉点分析也取得了较好的效果[13, 14, 15, 16],日本SELENE计划中,激光测高数据LALT的交叉点被用来对卫星轨道解算进行约束[17]。CE-1卫星采用200 km的高极轨道,飞行姿态稳定,有利于卫星精密定轨,其交叉点反映了不同轨在同一地形点上的差异,不论这些误差具体来源于轨道、测距、时间、环境或者仪器,如果能减小这些误差,即可改进地形精度。本文采用交叉点分析和平差的方法对CE-1激光高度计数据进行误差处理以生成高精度的DEM。由于月球卫星的轨道并不同于对地观测卫星的轨道,其轨道的稳定性和测控精度与地面卫星均有不同[18, 19, 20, 21],交叉点的计算和处理也会不同。
2 CE-1交叉点与交叉点不符值 2.1 交叉点和交叉点不符值卫星激光高度测量中,由于其周期运行,会在不同时间点经过目标天体同一位置,这个位置即为交叉点。由于轨道误差以及测量过程的误差,轨道在两次不同时间点上对同一地点产生的高度测量值出现偏差,这个偏差称为交叉点不符值(图 1)。假设t和t′为轨道在交叉点处的测量时间,h(t)和h(t′)为产生的2个高度测量值,则
即为交叉点不符值。 2.2 CE-1交叉点和不符值计算从中国科学院国家天文台获取CE-1运行期间的1397轨有效激光高度计2B级数据,去掉其中4条测高点少于50个的轨道,采用其中1393轨数据用于交叉点计算与分析。每轨数据约有5000~7000个有效测高点,连续激光点时间间隔为1 s,共约912万个激光高度测量点。计算交叉点数据之前,先对测高数据进行粗差剔除,假设相邻两点的高差为Δh,相邻两点的直线距离为Δd,当Δh/Δd>tan 60°时,该点被判定为粗差点,在本研究中,共去除31 077个粗差点。
由于月球引力场的特征和差异,月球卫星轨道不同于地球卫星冻结轨道[22, 23],不能利用轨道运行的周期性来求解交叉点的位置。本研究采用逐一线段检测的方式,判断每一轨上每两个连续的LAM轨迹点组成的线段是否与另一轨任两个连续的LAM轨迹点组成的线段相交,若线段相交,则为交叉点,并求解交点坐标,然后根据求出的交点坐标与线段端点的坐标,线性解算交点处在交叉两轨(i,j)上的时间值ti、tj。在每条轨道上选取交叉点左右邻近各3个点,共6个点处的LAM高度值,采用Akima插值[24]方法解算出交点处相应的高度值ai、aj。去除粗差点后的测高数据点,由于覆盖全月面,在交叉点计算中数据量非常大,因此采用了分块的策略分别计算每块区域内的交叉点。
3 CE-1交叉点特点分析 3.1 全球空间分布特性经过计算,全月面共产生1 620 056个交叉点,由于激光测高过程中有数据缺失的情况,为避免交叉点处参与插值的LAM点时间跨度过大而造成插值后交叉点处的高程出现较大的偏差,对相交线段端点处的LAM点进行了时间检核,保留了其中插值范围在7 s以内(对应月面距离约9800 m)的1 463 107个点,其不符值直方图分布如图 2。
CE-1卫星精密轨道由统一S波段(USB)测距测速和甚长基线干涉(VLBI)测量数据联合确定,轨道径向误差30 m[1, 18, 21]。而对不符值幅度300 m以内的1 150 346个点进行统计,其RMS值为120.05 m,可见交叉点不符值是轨道、测距、时间、环境或者仪器所造成的误差的综合反映。由于月球的自转,交叉点可以出现在任何纬度地带。CE-1采用的是极圆轨道,轨道倾角为88.2°[2],星下点轨迹在中低纬度地区与经线近似平行,随着轨道的不断运行,后期运行轨迹与前期轨迹的交叉点不断增多,因此,CE-1激光高度计产生的交叉点在中低纬度比较稀疏,而在两极区域分布非常密集。其中,93.1%的交叉点分布在南北纬80°以上的两极地区,如图 3(见文末)。对交叉点不符值的统计显示,97.5%不符值超过300 m的点分布在南北两极70°~90°范围以内。因此,对交叉点进行平差处理时,需要先对交叉点数据进行粗差剔除以避免交叉点不符值误差过大而造成高程纠正出现较大的偏差。
3.2 时间特性交叉点随时间变化的规律是确定平差模型的重要依据。图 4显示了一个地球日(2008年5月17日,第2180~2191轨)时段内交叉点不符值幅度的变化(a)和交叉点纬度随时间变化的分布(b),每轨交叉点的数量在两极位置出现高峰值。随着轨道的南北运行,不符值整体幅度呈现周期性变化,周期为2,即一轨两次变化,且相邻轨道上交叉点不符值的分布和变化规律相似,因此,可以采用基于时间参量的函数来描述不符值大小的变化。由于沿轨交叉点不符值幅度的变化类似于正余弦曲线,本文依据轨道运行周期,采用三角函数为不符值拟合的基本模型,并在三角函数模型基础上简化成多项式的形式进行交叉点不符值平差。
图 5为第0247、第0555、第2500、第2900轨上的交叉点不符值随时间的变化,其中轨道运行周期约为2 h 12 min,参考时间起点为每一轨数据第一个测高点的时间值。图 5显示了每轨变化规律在相似的基础上仍然存在差异,因此,可以对每轨数据采用相同形式但不同参数的模型进行处理,并通过时间参量的设置对相邻轨道不符值的变化进行约束。
4 CE-1 LAM交叉点平差模型本文参考了火星激光高度计数据中交叉点的处理方法[13],并根据不符值的定义和CE-1交叉点不符值沿轨分布特性,采用基于时间函数的模型来描述交叉点不符值的变化。为减小交叉点处测高点的偏差,假设f(t)和f(t′)分别为交叉点处两个不同时间点上高度值的修正量,则可以建立t和t′时刻的不符值误差方程
通过使Δd(t)TΔd(t)最小达到交叉点处高度测量的修正。其中,f(t)和f(t′)是基于交叉点处时间的函数,其形式为
即位于第j轨上的第i个交叉点,Gi为时间参数构成的系数矩阵,pj=[p0 p1 … pn]T为模型未知参数矩阵。考虑到不同模型形式对平差结果的影响,采用了4种不同形式的模型对试验区进行了试验,并对平差效果进行了对比,4种模型形式如下:模型1
模型2
模型3
模型4
由式(3)可知,每轨有(n+1)个待解参数pj,在每一个交叉点处,式(2)需要解算2(n+1)个未知系数。由于每一轨上的交叉点数据都可以建立式(2)的误差方程,因此有足够的观测值来解算所有的未知系数。对于所有的交叉点数据,G为一个非常庞大的稀疏矩阵,且值比较近似,解算pj的过程中会出现奇异矩阵。因此,设置pj的初始值为0,即假设所有轨道的高度测量均没有误差,并通过协方差逆阵Cpp-1对相邻轨道交叉点不符值的相似变化进行约束并迭代求解最终的pj值[25],迭代过程为
根据解算出每一轨的pj值,式(2)得以求解,然后依据激光高度测量点处的时间值获取相对应的f(t)值,则每一个测高点的纠正为
式中,t为该点时刻;h(t)0为最初获取的高度值;h(t)为改正之后的高度值。
5 试验结果及分析选择CE-1在月球(0°N~60°N,50°W~0°W)范围内的DEM进行交叉点平差试验,该区交叉点处理前生成的DEM有比较明显的条带现象(如图 6(a),见文末)。采用2阶多项式模型对该区产生的2593个交叉点进行了平差试验。结果显示,不符值残差由平差前的62.1 m降低至平差后的36.8 m。图 7为交叉点平差前后不符值的直方图分布,表 1显示了交叉点平差前后不符值分布范围的变化,与直方图分布显示的结果一致。其中,98.7%的交叉点不符值降低到了100 m以下,并且大部分处于30 m以下,而平差之前,仅有50%是在30 m以下。图 6(a)、(b)为运用该区的399 629个激光测高点生成的10 km分辨率的DEM效果,(a)为直接采用激光高度点进行克里金插值后的结果,(b)为采用2阶多项式进行交叉点平差纠正之后的激光高度点克里金插值生成的DEM。可以看出,(b)的条带现象较(a)有明显的改善,表明交叉点平差能明显改进LAM数据的精度和DEM的质量。
不符值范围 | >100 m | 50~100 m | 30~50 m | 10~30 m | 0~10 m |
平差前 | 9.99% | 20.32% | 19.09% | 30.24% | 20.36% |
平差后 | 1.27% | 4.17% | 7.98% | 32.66% | 53.91% |
为验证交叉点平差后整体精度的变化,将经过交叉点平差前后的DEM与美国LRO任务中LOLA产生的DEM进行了对比。LOLA采用5束激光进行测距,其获取的DEM平面精度可达5 m,高程方向约1 m[14],是目前发布的月球最高精度的DEM。在月球(0°N~60°N,50°W~0°W)范围内,LOLA约有2.69亿个高程点。将CE-1和LOLA的DEM统一到月固坐标系下,并分别采样成10 km分辨率等级,通过迭代最邻近点(iterative closest point,ICP)算法[26]对两个DEM进行整体配准,对两者X、Y、Z 3个方向的差异进行统计,结果如表 2。可看出,交叉点平差前后的DEM与LOLA的DEM的差异基本保持不变,平差前后的DEM均能与LOLA的DEM进行很好的配准,表明交叉点平差并未对整体地形产生变形。图 6为该区交叉点平差前后DEM与LOLA的DEM显示。同时,为了解交叉点平差对细节地形的纠正,以该区哥白尼撞击坑(8.17°N~11.23°N,21.63°W~18.57°W)(红框区域)的DEM为例,按CE-1点距1400 m的分辨率与LOLA进行对比(表 3),可看出,平差后的DEM与LOLA的DEM在XYZ方向上的差异明显缩小,且与LOLA的DEM进行配准时,Z方向的差异由原来的1.55 m降至-0.91 m,表明平差后的DEM与LOLA更为接近。需要说明的是,常用的DEM比较方法是在X、Y格网完全一致的情况下只统计Z方向的误差;而ICP配准属于刚体变换,使待配准DEM产生整体的三维旋转和平移,配准后的两个DEM的X、Y格网并不完全一致;为了避免再次重采样产生进一步的误差,采用最近距离的方法统计了X、Y、Z 3个方向的误差。
m | ||||||||
X | Y | Z | X | Y | Z | |||
平差前 | 平均值 | 88.20 | -46.23 | -65.41 | ICP | 0 | 0 | 0 |
DEM | 中误差 | 73.72 | 49.84 | 69.07 | 配准后 | 75.15 | 48.31 | 68.79 |
平差后 | 平均值 | 88.35 | -46.42 | -65.53 | ICP | 0 | 0 | 0 |
DEM | 中误差 | 73.27 | 49.93 | 69.09 | 配准后 | 74.79 | 48.29 | 68.83 |
m | ||||||||
X | Y | Z | X | Y | Z | |||
平差前 | 平均值 | 448.14 | -172.03 | -164.12 | ICP | -10.47 | -18.67 | 1.55 |
DEM | 中误差 | 209.42 | 388.34 | 453.76 | 配准后 | 209.17 | 375.67 | 342.44 |
平差后 | 平均值 | 94.65 | -43.74 | -83.67 | ICP | -9.74 | -19.81 | -0.91 |
DEM | 中误差 | 205.89 | 382.46 | 452.63 | 配准后 | 207.77 | 376.83 | 343.28 |
为分析基于时间模型的纠正效果和不同模型形式对结果的影响。表 4列出了采用4种不同形式的模型进行平差后的数值结果。从表中可以看出,4种模型均达到了减小不符值的效果,交叉点高程不符值中误差(RMS=)均降至了37 m以内。为显示4种模型的不同效果。截取了条带现象比较明显的一小块局部区域(18°N~44°N,33°W~19°W)进行DEM显示(图 8,见文末),4种模型都不同程度地减弱了原有DEM的条带现象,但是基于时间的2阶和3阶多项式模型产生的DEM效果更好。
不符值参数 | 中误差 | 均值 | 中值 | 最小值 | 最大值 |
平差前 | 62.11 | 2.35 | 1.68 | -298.50 | 293.82 |
模型1 | 36.30 | 0.37 | -0.07 | -242.04 | 270.13 |
模型2 | 34.35 | 0.44 | -0.23 | -245.51 | 230.11 |
模型3 | 36.83 | 0.14 | -0.10 | -240.27 | 268.92 |
模型4 | 32.84 | 0.34 | 0.012 | -229.01 | 259.10 |
在4种模型中,模型1、2、3每轨解算3个未知参数,模型4每轨解算4个,数值解算效果稍好。模型1和模型3的结果接近,但是在运用模型1的三角函数模型解算所得的参数进行式(9)高程纠正时,易出现局部区域的地形振荡现象而产生新的条带。其可能原因为三角函数模型比模型3的多项式模型变化的幅度大,在每轨交叉点位置处有较好的拟合,但当轨道覆盖的区域较长时,会在个别轨上出现较大的改正幅度而导致新的条带。模型2为三角函数模型的多项式近似表示,含2阶、4阶多项式和常数项,先对时间值进行归一化处理,t∈[-0.5,0.5],f(t)的值变化范围与模型1相同,其结果在模型1基础上有较好的改善,但不能完全消除这种局部新条带现象。模型3和模型4是模型2的变体,模型4的效果比模型3效果稍好,但每轨需要多解算一个未知数,计算量比模型3大,需要保证所有轨道都有至少4个交叉点参与计算。由于两者产生的DEM效果非常接近,当局部区域交叉点个数较少时,而模型3解算未知数要少,则更适用。
6 结 论本文根据CE-1激光高度计数据的特点,探讨了全月球范围内交叉点位置和高程计算的方法,对交叉点时空分布特性进行了分析,发现了交叉点在两极密集分布,中低纬度处密度较低,沿轨不符值随时间规律变化且相邻轨道变化规律相似。采用4种基于时间的函数模型在局部试验区域(0°N~60°N,50°W~0°W)内进行交叉点平差试验,均能使交叉点高程不符值中误差显著减小(从平差前的62.1 m降至37 m以内),在此基础上产生的DEM条带现象明显减弱,表明这些方法是可行的。
交叉点不符值的产生是由很多误差因素造成的,是轨道、测距、时间、环境或者仪器所造成的误差综合反映。在运用二阶多项式模型用于极区(60°N~90°N,0°W~40°W)和全月正面范围内的交叉点处理时,不符值均方根误差降低了9~14 m左右(分别为99.44 m至85.71 m,117.5 m至109.6 m),DEM效果有所改善,但不够显著。在交叉点产生的过程中,是根据沿轨邻近点的高度值插值产生的,且由于月球的转动和卫星轨道测控的精度,在CE-1卫星运行过程中,很多时段尤其是背面运行时段没有测控数据[2, 18, 21],卫星位置姿态的精度为交叉点分析带来了不确定性。因此,全月范围内的交叉点平差模型和方式还需要进一步考虑轨道精度在其中的影响和不符值模型的适用范围。在后期的工作中,将在依据时间建立不符值模型的基础上,进一步研究各种误差源的影响,并考虑交叉角度和地形因素对不符值的影响,在局部区域交叉点平差的基础上,建立全月交叉点整体平差模型。
[1] | PING Jinsong,HUANG Qian,YAN Jianguo,et al. Lunar Topographic Model CLTM-s01 from Chang’E-1 Laser Altimeter[J]. Science in China Series G:Physics,Mechanics & Astronomy, 2009, 52(7): 1105-1114. |
[2] | LI Chunlai,REN Xin,LIU Jianjun,et al. Laser Altimetry Data of Chang’E-1 and the Global Lunar DEM Model[J]. Science China:Earth Science,2010,40(3):281-293. (李春来, 任鑫, 刘建军, 等. 嫦娥一号激光测距数据及全月球DEM模型[J]. 中国科学: 地球科学, 2010, 40(3): 281-293.) |
[3] | CAI Zhanchuan,LIANG Yanyan,LI Jian,et al. Digital Elevation Model of the Moon from the Chang’E-1 Laser Altimeter[J]. Progress in Geophys, 2010, 25(4): 1153-1160. (蔡占川, 梁延研, 李坚, 等. 基于嫦娥一号卫星激光测高数据的月球数字高程模型[J]. 地球物理学进展, 2010, 25(4): 1153-1160.) |
[4] | ZHAO Shuangming,LI Deren,MOU Lingli. Inconsistency Analysis of CE-1 Stereo Camera Images and Laser Altimeter Data[J]. Acta Geodaetica et Cartographica Sinica,2011, 40(6): 751-755. (赵双明, 李德仁, 牟伶俐. CE-1立体相机影像与激光高度计数据不一致性分析[J]. 测绘学报, 2011, 40(6): 751-755.) |
[5] | DI Kaichang,HU Wenmin,LIU Yiliang,et al. Coregistration of Chang’E-1 Stereo Images and Laser Altimeter Data with Crossover Adjustment and Image Sensor Model Refinement[J]. Advances in Space Research, 2012, 50(12): 1615-1628. |
[6] | NISHIIMURA C E,FORSYTH D W. Improvements in Navigation Using Sea Beam Crossing Errors[J]. Marine Geophysical Research,1988,9(4): 333-352. |
[7] | TAI C K,FU L L. On Crossover Adjustment in Satellite Altimetry and Its Oceanographic Implications[J]. Journal of Geophysical Research,1986,91(C2): 2549-2554. |
[8] | KIM M. Theory of Satellite Ground-track Crossovers[J]. Journal of Geodesy, 1997, 71(12): 749-767. |
[9] | LUTHCKE S,MCCARTHY J J,PAVLIS D E,et al. Spaceborne Laser-altimeter-pointing Bias Calibration from Range Residual Analysis[J]. Journal of Spacecraft and Rockets,2000,37(3):374-384. |
[10] | LI Jiancheng,CHEN Junyong. Determination of Gravity Anomalies over the South China Sea by Combination of TOPEX/Poseidon,ERS2 and Geosat Altimeter Data[J]. Acta Geodaetica et Cartographica Sinica, 2001, 30(3): 192-202. (李建成,陈俊勇. 联合 TOPEX/Poseidon,ERS2 和 Geosat 卫星测高资料确定中国近海重力异常[J]. 测绘学报,2001,30(3): 192-202.) |
[11] | 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. (黄谟涛,王瑞,翟国君,等. 多代卫星测高数据联合平差及重力场反演[J]. 武汉大学学报:信息科学版,2007,32(11): 988-993.) |
[12] | ZHAI Guojun. The Theory and Method of Satellite Altimeter Data Processing[M]. Beijing:Surveying and Mapping Press, 2000. (翟国君. 卫星测高数据处理的理论与方法[M]. 北京:测绘出版社, 2000.) |
[13] | NEUMANN G A,ROWLANDS D D,LEMOINE F G,et al. Crossover Analysis of Mars Orbiter Laser Altimeter Data[J]. Journal of Geophysical Research, 2001, 106(E10): 23753-23768. |
[14] | MAZATICO E,NEUMANN G A,ROWLANDS D D,et al. Geodetic Constraints from Multi-beam Laser Altimeter Crossovers[J]. Journal of Geodesy,2010,84(6): 343-354. |
[15] | MAZATICO E,ROWLANDS D D,NEUMANN G A,et al. Selenodesy with LRO:Radio Tracking and Altimetric Crossovers to Improve Orbit Knowledge and Gravity Field Estimation[C]//Proceedings of the 42nd Lunar and Planetary Science Conference.The Woodlands: [s.n.], 2011: 2215-2216. |
[16] | SMITH D E,ZUBER M T,JACKSON G B,et al. The Lunar Orbiter Laser Altimeter Investigation on the Lunar Reconnaissance Orbiter Mission[J]. Space Science Reviews, 2010, 150(1): 209-241. |
[17] | GOOSSENS S,MATSUMOTO K,ROWLANDS D D,et al. Orbit Determination of the SELENE Satellites Using Multi-satellite Data Types and Evaluation of SELENE Gravity Field Models[J]. Journal of Geodesy, 2011, 85(8): 487-504. |
[18] | QIAO Shubo,LI Jinling,SUN Fuping. Application Analysis of Lunar Exploration Satellite Positioning by VLBI Technique[J]. Acta Geodaetica et Cartographica Sinica,2007, 36(3): 262-268. (乔书波,李金岭,孙付平. VLBI 在探月卫星定位中的应用分析[J]. 测绘学报,2007, 36(3): 262-268.) |
[19] | YAN Jianguo,PING Jinsong,LI Fei,et al. Chang’E-1 Precision Orbit Determination and Lunar Gravity Field Solution[J]. Advances in Space Research,2010,46(1): 50-57. |
[20] | HUANG Qian,PING Jingsong,YAN Jianguo,et al. Chang’E-1 Laser Altimetry Data Processing[J]. Advances in Geosciences: Planetary Science, 2010, 19: 137-149. |
[21] | CHEN Ming,TANG Geshi,CAO Jianfeng,et al. Precision Orbit Determination of CE-1 Lunar Satellite[J]. Geomatics and Information Science of Wuhan University,2011,36(2): 212-217. (陈明,唐歌实,曹建峰,等. 嫦娥一号绕月探测卫星精密定轨实现[J]. 武汉大学学报: 信息科学版,2011,36(2): 212-217.) |
[22] | HUANG Yong. Orbit Determination of the First Chinese Lunar Exploration Spacecraft CE-1[D]. Shanghai:Shanghai Astronomical Observatory of Chinese Academy of Sciences,2006. (黄勇. “嫦娥一号” 探月飞行器的轨道计算研究[D]. 上海:中国科学院上海天文台,2006.) |
[23] | LIU Lin,LIU Shiyuan,WANG Yanrong. The Frozen Orbit of the Orbiter of the Major Planet (or Moon)[J]. Journal of Spacecraft TT&C Technology, 2003, 22(2): 19-24. (刘林, 刘世元, 王彦荣. 关于大行星(或月球)轨道器的冻结轨道[J]. 飞行器测控学报, 2003, 22(2): 19-24.) |
[24] | AKIMA H. A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures[J]. Journal of the ACM, 1970, 17(4): 589-602. |
[25] | TARANTOLA A,VALETTE B. Generalized Nonlinear Inverse Problems Solved Using the Least Squares Criterion[J]. Reviews of Geophysics and Space Physics, 1982, 20(2): 219-232. |
[26] | BESL P J,MCKAY N D. A Method for Registration of 3-D Shapes[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1992, 14(2): 239-256. |