一种预先确定星载激光测高仪的足印位置的方法与流程

本发明涉及星载激光测高仪在轨几何检校领域，特别涉及一种预先确定星载激光测高仪的足印位置的方法，应用于星载激光测高仪在轨几何检校过程中的激光地面探测器布设位置的预报。

背景技术：

在轨几何检校试验的成功开展是星载激光测高仪获得高精度测高数据的必要条件之一。星载激光测高仪足印位置的预报是在轨几何检校的第一步，也是极为关键的一步。预报的精度不仅关系着整个检校试验外业布设工作的难易程度，甚至决定着试验的成败。

目前，仅有美国ICESat-GLAS(Ice Cloud Elevation Satellite Geoscience Laser Altimeter System)开展过相关的研究和工程实践，其激光脉冲发射重频为40Hz，激光沿轨相邻两个足印中心间距为170米左右，地面足印平均大小为65米左右，沿轨布设大于170米的探测器阵列，必然可以捕获到光斑。由于其硬件的优势，使其在足印预报方法上并不需要开展过多的研究。但是其高重频是根据卫星的使用目的而设计的，即对冰川的厚度进行高精度测量；而对于光学测绘卫星而言，星上本身携带了高分辨率相机、姿态测控系统和轨道确定系统等多个载荷，激光测高仪的脉冲出射频率必定要控制在一定的范围内，以保证整星的能源分配问题，且用于符合光学测绘的激光测高数据在实际应用中也不需要过高的激光出射频率。低频率的脉冲出射为星载激光测高仪的在轨几何检校增加了试验难度，最突出地体现在激光足印的预报上。由于国际上暂时没有涉及该方面技术的测绘卫星在轨运行，激光足印预报方面的研究几乎属于空白。

因此，本领域需要一种能够准确预报星载激光测高仪足印位置的方法。

技术实现要素：

为此，本发明提出一种预先确定星载激光测高仪的足印位置的方法，弥补了目前激光测高仪高精度足印位置预报方法的空白，并解决了现有简单预报技术的限制和缺陷导致的一个或多个问题。

本发明另外的优点、目的和特性，一部分将在下面的说明书中得到阐明，而另一部分对于本领域的普通技术人员通过对下面的说明的考察将是明显的或从本发明的实施中学到。通过在文字的说明书和权利要求书及附图中特别地指出的结构可实现和获得本发明目的和优点。

本发明提供了一种预先确定星载激光测高仪的足印位置的方法，其特征在于，所述方法包括以下步骤：

步骤1，计算激光测高仪在卫星本体坐标系中的指向M_{b_l}，步骤1具体包括以下子步骤：

步骤1.1，预估激光测高仪与平台之间的安装角(α_r,β_r)；

步骤1.2，基于步骤1.1预估的安装角，计算激光测高仪在卫星本体坐标系中的指向M_{b_l}；

步骤2，将激光测高仪的指向从卫星本体坐标系旋转至轨道坐标系M_{o_l}；

步骤3，将激光测高仪的指向从轨道坐标系旋转至地球固定坐标系，获得激光测高仪在地球固定坐标系下的指向M_{T_l}；

步骤4，通过计算地球固定坐标系下的激光测高仪的指向与地球椭球的交点，确定星载激光测高仪的足印位置，其中，所述步骤4具体包括以下子步骤：

步骤4.1，根据下述公式计算激光指向与地球椭球相交的地球固定坐标系下的地面点坐标(X₁,Y₁,Z₁)：

其中，M_T＝[X₁ Y₁ Z₁]^T，X₁、Y₁和Z₁为待求的激光指向与地球椭球相交的地球固定坐标系下的地面点坐标，μ为待求的比例系数，M_{T_l}为激光测高仪在地球固定坐标系下的指向，P₀为卫星激光质心在地球固定坐标系的坐标值，a和b分别为地球椭球长半轴和短半轴长度，h为地面大地高程；

步骤4.2，根据下述公式将地球固定坐标系下的地面点坐标(X₁,Y₁,Z₁)转换成大地经纬度坐标系下的地面点坐标(L,B,H)：

其中，e₁和e₂分别为第一偏心率和第二偏心率，a和b分别为地球椭球长半轴和短半轴长度，X₁、Y₁、Z₁为通过步骤4.1计算出的激光指向与地球椭球相交的地球固定坐标系下的地面点坐标。

本发明还提供了一种计算机程序，其中，当执行所述计算机程序时，执行本发明任一实施例所述的方法。

本发明针对星载激光测高仪足印位置预报领域的空白，提供了一种预先确定星载激光足印位置的方法，该方法结合了传统的光学相机成像的原理，并充分顾及了激光测高仪自身的测高机理，预报结果能够较好地符合卫星真实的飞行过境情况，具有较高的预报精度，并且极大限度地减少不必要的人力物力，节省了大量的试验开销。排除激光指向角、卫星平台姿态预报数据和卫星质心轨道预报数据精度损失的情况下，该方法的星载激光测高仪足印位置预报精度能够控制在5米。

附图说明

图1为根据本发明实施例的确定星载激光测高仪的足印位置的方法的流程图。

具体实施方式

下面参照附图对本发明进行更全面的描述，其中说明本发明的示例性

实施例。

根据本发明的实施例，可以随机抽取星载激光测高仪进行足印位置预报模型的构建；优选的，本发明的实施案例以资源三号02星的激光测高仪的足印位置预报模型构建为例，但并不限于此。

如图1所示，激光测高仪预报模型的构建和光学严密几何成像模型的构建存在同样的过程，涉及到各个坐标系之间的旋转，通过各坐标系的旋转最终将星上指向传输至地面点，本发明所提出的星载激光测高仪足印位置预报模型构建方法具体包括以下步骤：

步骤1，计算激光测高仪在卫星本体坐标系中的指向M_{b_l}，步骤1具体包括以下子步骤：

步骤1.1，预估激光测高仪与平台之间的安装角(α_r,β_r)；

由于卫星发射升空过程中的剧烈震荡和在轨后太空环境与地面实验室的差异较大，会引起激光测高仪与平台之间的安装关系发生变化，导致激光测高仪的安装角发生变化，因此需要计算真实的激光测高仪与平台之间安装角。

参考资源三号02星激光测高仪设计指标，获取激光测高仪的实验室标定的安装角α＝90.86229°，β＝0.861985°。以实验室标定安装角为基础，采用基于地形匹配的安装角预估方法计算真实的指向角。这里采用金字塔计算策略，逐步缩小计算区域和遍历间隔，将时间限定在一定的区间内，保证计算的精度。将金字塔分为3层来进行计算，第一层金字塔的区域间隔为2°，遍历间隔为0.1°，第二层金字塔的区域间隔为0.4°，遍历间隔为1′，第三层金字塔的区域间隔为0.02°，遍历间隔为1″，具体的计算公式如下式所示：

其中，(α_r,β_r)为预估的安装角(其中，α_r为激光指向与本体坐标系X轴之间的夹角，β_r为激光指向与本体坐标系Z轴之间的夹角)，表示从1层至3层逐层计算，α_i、β_i表示第i层遍历的安装角，表示括号内第i层取最小值时对应的α_i、β_i，并将其赋值给α₀、β₀，当i＝1首次计算时，α₀、β₀的值取实验室测量值90.86229°、0.861985°，q(i)为第i层遍历间隔，p(i)为第i层的区域间隔，表示对括号内进行遍历，行遍历序号m，列遍历序号n，E(α_i,β_i)为α_i、β_i指向下计算某轨的高程值与已知高程值之间的标准差，具体的高程值计算公式可参照激光测高仪的严格几何测高模型，已知高程采用公开全球范围的AW3D30数据。

以资源三号02星的第382轨激光数据为例，计算过程如下：首先进行第一层金字塔遍历，区域间隔为2°，遍历间隔为0.1°，

经过第一层的遍历，可以获得第一层的预估结果α₀、β₀分别为89.9°、0.0°，下面进行第二层金字塔的遍历，区域间隔为0.4°，遍历间隔为1′，

经过第二层的遍历，可以获得第二层的预估结果α₀、β₀分别为89.945°、0.046°，下面进行第三层金字塔的遍历，区域间隔为0.02°，遍历间隔为1″，

通过上式可以求得预估的安装角α_r＝89.94445°，β_r＝0.04692°

步骤1.2，基于步骤1.1预估的安装角，根据下述公式计算激光测高仪在卫星本体坐标系中的指向M_{b_l}：

式中，M_{b_l}为步骤1待求的激光测高仪在卫星本体坐标系中的指向，α_r、β_r为步骤1.1求得的激光测高仪的安装角。

步骤1.1求得的激光测高仪安装角预估值代入到上式可以求得预估指向

步骤2，将激光测高仪的指向从卫星本体坐标系旋转至轨道坐标系M_{o_l}。

根据光束的物理传播路径，传输的第一步是将激光测高仪的指向从卫星本体坐标系旋转至轨道坐标系，步骤2具体包括以下子步骤：

步骤2.1，估计卫星平台的本体坐标系相对于轨道坐标系的姿态值数值单位为弧度，其中，roll、pitch、yaw分别为卫星本体坐标系相对于轨道坐标系的三轴欧拉角：横滚角、俯仰角和偏航角，即，roll为卫星本体坐标系相对于轨道坐标系的横滚角，pitch为卫星本体坐标系相对于轨道坐标系的俯仰角，yaw为卫星本体坐标系相对于轨道坐标系的偏航角。

一般情况下，由于卫星姿控系统的高精度实时控制，平台的本体坐标系相对于轨道坐标系的姿态变化较小，因此这里采用了取零值的方法，即三轴姿态均为零

步骤2.2，根据下述公式计算本体坐标系相对于轨道坐标系的旋转矩阵M_{o_b}；

a₁＝cos(pitch)*cos(yaw)

a₂＝-cos(pitch)*sin(yaw)

a₃＝sin(pitch)

b₁＝-sin(roll)*sin(pitch)*cos(yaw)+cos(roll)*sin(yaw)

b₂＝sin(roll)*sin(pitch)*sin(yaw)+cos(roll)*cos(yaw)

b₃＝sin(roll)*cos(pitch)

c₁＝-cos(roll)*sin(pitch)*cos(yaw)-sin(roll)*sin(yaw)

c₂＝cos(roll)*sin(pitch)*sin(yaw)-sin(roll)*cos(yaw)

c₃＝cos(roll)*cos(pitch)

将步骤2.1得到的卫星平台的本体坐标系相对于轨道坐标系的姿态值(即，)代入到M_{o_b}，可得本体坐标系相对于轨道坐标系的旋转矩阵

本实施例的步骤2.1优选的采用了的方法，即代入到步骤2.2计算。

步骤2.3，根据下述公式将激光测高仪的指向由本体坐标系转换到轨道坐标系；

M_{o_l}＝M_{o_b}*M_{b_l}

其中，M_{o_l}为待求的激光测高仪在轨道坐标系下的指向，M_{o_b}为步骤2.2计算的本体坐标系相对于轨道坐标系的旋转矩阵，M_{b_l}为步骤1.2计算的激光测高仪的指向。

将步骤1.2和步骤2.2计算得到的激光测高仪载荷与平台之间的预估指向与本体坐标系相对于轨道坐标系的旋转矩阵相乘，可以得到激光测高仪的指向由本体坐标系转换到轨道坐标系的向量

步骤3，将激光测高仪指向从轨道坐标系旋转至地球固定坐标系。

根据光束的物理传播路径，传输的第二步是将激光测高仪指向从轨道坐标系旋转至地球固定坐标系，步骤3具体包括以下子步骤：

步骤3.1，计算卫星激光质心在地球固定坐标系的坐标值P₀＝[X₀ Y₀ Z₀]^T和位移速度值V₀＝[V_X0 V_Y0 V_Z0]^T；

根据本发明的优选实施例，步骤3.1的具体计算过程如下：

式中，为指定时刻的惯性系卫星位置关于时间的二阶导数即加速度，为待求量t为某一指定时刻，r为指定时刻的惯性系卫星位置，v为指定时刻的惯性系卫星速度，C_D为大气阻力系数，C_R为太阳光压系数，a_iR为作用在卫星轨道径向上的第i段经验加速度，ai_T为作用在卫星轨道切向上的第i段经验加速度，a_iN为作用在卫星轨道法向上的第i段经验加速度，e_R(t)为指定时刻卫星轨道径向方向单位向量，e_T(t)为指定时刻卫星轨道切向方向单位向量，e_N(t)为指定时刻卫星轨道法向方向单位向量，t_i为第i段经验加速度参数开始时刻，t_i+1为第i+1段经验加速度开始时刻，也即第i段经验加速度结束时刻。

参考时刻卫星位置速度、大气阻力系数C_D、太阳光压系数C_R、经验加速度参数a_iR、a_iT、a_iN等来自于预报区间之前的精密轨道拟合。

通过上式可求得新的加速度因此可以利用下式求待估参数(P₀,V₀)，已知初值位置和速度P₀'，V₀'，C_D，C_R及加速度，即

式中，C_D和C_R作为常值估计。

根据待预报地区，粗略估计卫星过境时间t，利用轨道预报计算公式，结合卫星过境前24小时的长时段轨道数据作为初始值，对过境的轨道参数进行预报。由于该预报一个渐进递推的长时段处理，过程涉及的数据较多，这里仅给出一例作为参考，24小时前的轨道数据初值为：P₀'＝[-728959.1 6827631.7 473238.5]^T和V₀'＝[1543.3 -367.0 7528.4]^T。结合其他可获取的参数，可计算得到过境预报区域的卫星的激光质心在地球固定坐标系的坐标值P₀＝[-1855244.6 4669501.6 4693461.4]^T和位移速度值V₀＝[-287.4 5397.1 -5468.8]^T。

步骤3.2，根据下述公式计算轨道坐标系相对于地球固定坐标系的旋转矩阵M_{T_o}；

其中，P₀＝[X₀ Y₀ Z₀]^T为卫星激光质心在地球固定坐标系的坐标值，V₀＝[V_X0 V_Y0V_Z0]^T为卫星激光质心在地球固定坐标系的位移速度值，A₁,A₂,A₃,B₁,B₂,B₃,C₁,C₂,C₃为旋转矩阵M_{T_o}的各个元素。

将步骤3.1得到的过境预报区域的卫星的激光质心在地球固定坐标系的坐标值P和位移速度值V代入到步骤3.2的公式，可以求得轨道坐标系相对于地球固定坐标系的旋转矩阵

步骤3.3，根据下述公式将激光测高仪的指向由轨道坐标系转换到地球固定坐标系；

M_{T_l}＝M_{T_o}*M_{o_l}

式中，M_{o_l}表示激光测高仪在轨道坐标系下的指向，M_{T_o}表示轨道坐标系相对于地球固定坐标系的旋转矩阵，M_{T_l}表示激光测高仪在地球固定坐标系下的指向。

将步骤2.3和步骤3.2计算得到的激光测高仪的指向由本体坐标系转换到轨道坐标系的旋转矩阵与轨道坐标系相对于地球固定坐标系的旋转矩阵相乘，可以得到激光测高仪的指向由本体坐标系转换到地球固定坐标系的向量

步骤4，通过计算地球固定坐标系下的激光测高仪的指向与地球椭球的交点，确定星载激光测高仪的足印位置。

根据光束的物理传播路径，在激光到达地面后，需要求地球固定坐标系下的激光测高仪的指向与地球椭球的交点，所述交点即本发明要确定的星载激光测高仪的足印位置，即激光探测器的布设位置的坐标。步骤4具体包括以下子步骤：

步骤4.1，根据下述公式计算激光指向与地球椭球相交的地球固定坐标系下的地面点坐标(X₁,Y₁,Z₁)：

其中，M_T＝[X₁ Y₁ Z₁]^T，X₁、Y₁和Z₁为待求的激光指向与地球椭球相交的地球固定坐标系下的地面点坐标，μ为待求的比例系数，^MT_l为激光测高仪在地球固定坐标系下的指向，P₀为卫星激光质心在地球固定坐标系的坐标值，a和b分别为地球椭球长半轴和短半轴长度(a＝6378137.0米，b＝6356752.3米)，h为地面大地高程，可以通过其他辅助高程数据获取，初始值给0。

其中，M_T＝P₀+μ*M_{T_l}可以表示为：或表示为将X₁、Y₁、Z₁代入到公式中，可以得到下述方程

式中，X₀、Y₀和Z₀为步骤3.1计算得到卫星激光质心在地球固定坐标系的坐标值，M_{T_l}为步骤3.3计算得到的值。

解上述一元二次方程，得到比例系数μ，因为方程会有两个解，最终选择绝对值小的u＝506437.3为最终解，根据M_T＝P₀+μ*M_T__l从而可以计算出M_T＝[-1718742.3 4325848.3 4347414.8]^T。

步骤4.2，根据下述公式将地球固定坐标系下的地面点坐标(X₁,Y₁,Z₁)转换成大地经纬度坐标系下的地面点坐标(L,B,H)；

(X₁,Y₁,Z₁)对于我们来说，没法直接使用，需要转换为我们经常用的经纬度格式(L,B,H)，这个也是从天上到地面的必要转换。

其中，e₁和e₂分别为第一偏心率和第二偏心率，a和b分别为地球椭球长半轴和短半轴长度，X₁、Y₁、Z₁为通过步骤4.1计算出的激光指向与地球椭球相交的地球固定坐标系下的地面点坐标。

根据上述公式，可以求得地面坐标(L,B,H)＝(111.66887,43.23643,1079.99)。

根据本发明的优选实施例，为了得到更精确的星载激光测高仪的足印位置，步骤4还可以包括：

步骤4.3，通过迭代计算，获得精确的激光探测器的布设位置；步骤4.3具体包括：将步骤4.2计算得到的地面点坐标(L,B)代入到高程数据中，重新内插得到更精确的高程值H'，将计算得到的H'代入到步骤4.1和步骤4.2进行迭代运算，以得到精确的星载激光测高仪的足印位置。

具体来说，将步骤4.2得到的地面点坐标(L,B)代入到高程数据中，重新内插得到更精确的高程值H'，内插公式采用双线性内插：

H'＝(1-Δx)(1-Δy)H₁+(1-Δx)ΔyH₂+Δx(1-Δy)H₃+ΔxΔyH₄

式中，H₁、H₂、H₃和H₄为H'附近x、y方向上相邻点的高程值，Δx为H'距H₃和H₄的x方向距离差，Δy为H'距H₁和H₂的y方向距离差。

将计算得到的H'代入到步骤4.1和步骤4.2，重复计算，此过程一般迭代三次，最终得到精确的星载激光测高仪的足印位置(L,B,H)，即，足印中心点的坐标。

将步骤4.2计算得到的地面点坐标(111.66887,43.23643)代入到我们使用的AW3D30数据(公开数据)中，采用上述内插公式，可以内插得到新的高程值H＝951.21；将H＝951.21代入到步骤4.1和步骤4.2得到新的地面点坐标(111.66879,43.23638，951.21),将(111.66879,43.23638)代入到AW3D30数据中，可以内插得到新的高程值H＝950.14；将H＝950.14代入到步骤4.1和步骤4.2得到新的地面点坐标(111.66879,43.23638，950.08)，该位置即为待铺设探测器阵的中心点。

本发明结合了传统的光学相机成像的原理，并充分顾及了激光测高仪自身的测高机理，预报结果能够较好地符合卫星真实的飞行过境情况，具有较高的预报精度。

以上内容仅为本发明的较佳实施例，对于本领域的普通技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，本说明书内容不应理解为对本发明的限制。