CN104864853B - 一种高分辨率三线阵卫星沿轨方向姿态颤振的探测方法 - Google Patents

Abstract

本发明涉及一种高分辨率三线阵卫星沿轨方向姿态颤振的探测方法，包括以下步骤：1)三线阵影像同名点匹配，并进行匹配粗差剔除；2)获取三线阵影像沿轨方向投影偏差，并建立颤振估计函数模型；3)三线阵影像沿轨方向姿态颤振估计。与现有技术相比，本发明具有方法先进、计算精确、应用范围广等优点。

Description

一种高分辨率三线阵卫星沿轨方向姿态颤振的探测方法

技术领域

本发明涉及一种卫星姿态颤振的探测方法，尤其是涉及一种高分辨率三线阵卫星沿轨方向姿态颤振的探测方法。

背景技术

平台姿态的稳定性，是决定高分辨率卫星几何定位精度的关键因素之一，而卫星平台颤振是高分辨率卫星在轨运行普遍存在的复杂现象，对成像几何精度的影响是一个尚未解决的科学难题。许多遥感卫星受到姿态颤振的影响导致几何定位精度大幅下降，平台颤振是国内外遥感卫星在几何定位方面共同面临的挑战(Iwasaki,2011)。

解决姿态颤振最常规的方法是在卫星上安装具有更高测量精度和更高采样频率的姿态传感器，但无疑增加了卫星的制造成本，最重要的是，对于已在轨运行的卫星，该方法则是无法实现的(Lehner and Müller,2003；Schwind et al.,2009)。另一种姿态颤振探测方法是，利用密集的地面控制点改正卫星姿态颤振(Shin et al.,1997)。Gwinner etal.(2010)开发了连续摄影测量平差(Sequential Photogrammetric Adjustment)方法，利用激光高度计数据生成的DEM(数字高程模型)作为地面控制数据，用于改正火星快车的HRSC影像受到的颤振影响，结果显示颤振改正后生成的DEM中不再包含周期性的误差。这些方法依赖于外部地面控制数据，而控制数据的获取需要较高的人力和经济成本。Tong etal.(2015c)提出一种基于三线阵影像的垂轨方向颤振探测方法，该方法能够有效探测三线阵传感器在垂轨方向上受到的姿态颤振影响，但是无法探测沿轨方向姿态颤振。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种方法先进、计算精确、应用范围广的高分辨率三线阵卫星沿轨方向姿态颤振的探测方法。

本发明的目的可以通过以下技术方案来实现：

一种高分辨率三线阵卫星沿轨方向姿态颤振的探测方法，包括以下步骤：

1)三线阵影像同名点匹配，并进行匹配粗差剔除；

2)获取三线阵影像沿轨方向投影偏差，并建立颤振估计函数模型；

3)三线阵影像沿轨方向姿态颤振估计。

所述的步骤1)具体包括以下步骤：

11)根据SIFT匹配算法，获取初始的同名点，采用相对定向的方法，消除三视影像之间相对的系统性偏差，使得同名摄影光线相交于同一物方点；

12)将下视影像作为参考影像，以N×N像素大小的空间格网点作为待匹配点，利用基于几何约束的互相关方法在三视影像上进行同名点的匹配；

13)通过设定相关系数阈值剔除部分粗差点，结合相对定向参数，通过前方交会获取三线阵同名点对应的物方点坐标，并将物方点分别投影至三个影像的像方空间，计算相应的匹配点与投影点之间的像方偏差d′，统计像方偏差的平均值μ和标准差σ，剔除粗差点公式为：

所述的步骤2)具体包括以下步骤：

21)利用同名点对三线阵影像进行相对定向处理，结合相对定向参数，将前、后视影像中受到颤振影响的同名点，通过前方交会获取受到颤振影响的地面点G′的坐标，并将地面点G′后向投影至下视影像的像方空间，获得地面点坐标的投影点；

22)计算投影点与相应的下视匹配像点之间的偏差，即三线阵影像的系统性误差和颤振误差共同引起的视差d′_x123：

其中，Δx₁为姿态角Pitch在下视影像上造成的像方偏差，Δx₂为姿态角Pitch在前视影像上造成的像方偏差，Δx₃为姿态角Pitch在后视影像上造成的像方偏差，κ₁为下视影像的分辨率，κ₂为前视影像的分辨率，κ₃为后视影像的分辨率，x₁和y₁为下视影像上同名点的像方坐标，x₂和y₂为前视影像上同名点的像方坐标，x₃和y₃为后视影像上同名点的像方坐标，a₁,b₁,c₁,a₂,b₂,c₂,a₃,b₃,c₃分别为下、前、后视影像系统误差的仿射模型系数，f为焦距，t₁，t₂和t₃分别为下、前、后视影像曝光时刻，β₂为前视传感器视线矢量与下视传感器视线矢量的夹角，β₃为后视传感器视线矢量与下视传感器视线矢量的夹角，J_Pitch(t)为轨道姿态角Pitch颤振函数模型，t为卫星的飞行时间，A_i为第i个谐波分量的振幅，ω_i为第i个谐波分量的角频率，为第i个谐波分量的振幅，角频率和初始相位值；

23)利用基于仿射补偿模型的相对定向，消除三线阵影像之间相对的系统性误差，得到投影偏差d_x123：

24)建立颤振估计函数模型为：

其中，V_x为同名点行坐标的残差向量。

所述的步骤3)具体包括以下步骤：

根据颤振估计函数模型，建立最优化目标函数，并根据最优化估计算法，估计颤振模型中的参数，使得目标函数最小，得到姿态角Pitch颤振函数模型J_Pitch(t)的参数结果，包括谐波分量的振幅、角频率和初始相位值，消除由于颤振造成的行坐标残差向量，

所述的目标函数为：

与现有技术相比，本发明具有以下优点：

1)方法先进：在投影偏差中存在地形信息和颤振信息，无法明确区分二者，本方法有效分离地形视差和颤振视差，建立了颤振估计函数模型，解决了沿轨方向颤振探测受地形影响的难题；

2)计算精确：本发明采用全局最优化估计算法，颤振模型参数可获得最优值，且估计过程中无需设定初始值；

3)应用范围广：本发明能够应用于具有较大基高比的三线阵立体影像的颤振探测。

附图说明

图1为姿态角Pitch颤振探测原理图。

图2为沿轨方向颤振探测流程图。

图3为A景影像缩略图，其中，图(3a)为A景下视影像subCCD1～3影像缩略图，图(3b)为A景前视影像subCCD1～3影像缩略图，图(3c)为A景后视影像subCCD1～3影像缩略图。

图4为B景影像缩略图，其中，图(4a)为B景下视影像subCCD1～3影像缩略图，图(4b)为B景前视影像subCCD1～3影像缩略图，图(4c)为B景后视影像subCCD1～3影像缩略图。

图5为A和B景立体像对之间受地形影响的像方投影偏差，其中，图(5a)为A景沿轨方向下视与前视的视差，图(5b)为A景沿轨方向下视与后视的视差，图(5c)为B景沿轨方向下视与前视的视差，图(5d)为B景沿轨方向下视与后视的视差。

图6为A、B景三线阵影像之间沿轨方向消除地形影响后的投影偏差图，其中，图(6a)为A景原始投影偏差图，图(6b)为B景原始投影偏差图，图(6c)为A景相对定向后投影偏差图，图(6d)为B景相对定向后投影偏差图，图(6e)为A景颤振补偿后投影偏差图，图(6f)为B景颤振补偿后投影偏差图。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。

实施例：

影响在沿轨方向上的像方偏差主要是由于姿态角Pitch引起，假设姿态测量值为卫星不受颤振影响的理想轨道参数，实际姿态值为受颤振影响的轨道参数，姿态角Pitch颤振探测原理如图1所示。

图中，不受颤振影响的理想轨道下，地面点G在下、前和后视影像像方空间上分别成像于理想像点和设下、前和后视传感器的曝光时刻分别为t₁，t₂和t₃。而在实际姿态下，姿态角Pitch发生变化，理想像点分别偏移至受颤振影响的像点p₁，p₂和p₃，则下、前、后视影像上，受颤振影响像点与理想像点之间的距离即为姿态角Pitch变化造成的像方偏差分别为-|Δx₁|，-|Δx₂|和|Δx₃|，其对应的物方空间的偏差为：

其中，和分别为像方对应在物方空间上的偏差，κ₁，κ₂和κ₃分别为下、前和后视影像的分辨率。

在图1中，物方点G′为前、后视影像受到颤振影响的像点在物方上的交会点，则物方点G′中既包含了地形信息，也包含了前、后视影像分别受到的姿态角Pitch颤振信息，物方点G′在下视影像上的后向投影点p′₂₃与下视影像上受颤振影响的像点之间的投影偏差为|d_x123|，则：

由于前、后视传感器与下视传感器之间的夹角相等，所以在物方空间上：

则：

投影偏差值|d|_x123可以表示为：

将式(1)第一式，式(4)代入式(5)得：

考虑像方偏差Δx₁，Δx₂和Δx₃的正负号，三线阵影像之间投影偏差d_x123可以表示为：

从上式中可以看出，三线阵影像之间的投影偏差d_x123是由三线阵影像对应的像方偏差共同决定，也即三线阵影像受到的姿态颤振影响决定了投影偏差d_x123，可见投影偏差中不包含任何地形信息。

沿轨方向颤振探测采用与垂轨方向颤振探测相同的同名点，利用前、后卫星平台颤振是一种复杂的平台振动现象，简单的谐波函数无法建模平台颤振，因此，根据傅里叶分析原理，将平台颤振的波形分解为多个简单谐波的叠加(Schwind et al.,2009；Tong etal.,2014；)，如下表示：

式中，J(t)表示随时间变化的卫星平台姿态颤振函数，t为卫星的飞行时间，A_i，ω_i和分别为第i个谐波分量的振幅，角频率和初始相位值。

由于姿态角Pitch变化导致的下视、前视和后视像方偏差Δx₁，Δx₂和Δx₃可表示为：

式中，f表示焦距，Δβ₁，Δβ₂和Δβ₃分别为在t₁，t₂和t₃时刻，下、前和后视传感器对应的姿态角Pitch颤振值。

由于卫星平台姿态颤振随时间变化，导致三线阵影像之间的投影偏差d_x123也随时间发生变化，并且三线阵影像在相同的轨道上拍摄，其对同一地物的成像时间间隔较短。沿轨方向像方偏差仅受姿态角Pitch的影响，故下、前和后视影像受到的卫星平台姿态角Pitch颤振的影响具有很强的连续性和相关性。因此，前、下和后视影像受到的周期性颤振影响可以用同一个颤振模型来建模。故结合式(8)，式(9)改写为考虑曝光时刻的一般形式：

式中，J_Pitch(t)为该段轨道姿态角Pitch颤振函数模型。与垂轨方向类似，沿轨方向像方也受到其它误差的影响，为保证颤振估计模型的一般性，在式(10)中加入仿射模型表述系统性误差。

式中，x₁,y₁,x₂,y₂,x₃,y₃分别为下、前、后视影像上同名点的像方坐标，a₁,b₁,c₁,a₂,b₂,c₂,a₃,b₃,c₃分别为表示下、前、后视影像系统误差的仿射模型系数。

考虑系统性误差的影响，同时结合式(7)和式(11)，三线阵影像的系统性误差和颤振误差共同引起的投影偏差d′_x123可重写成考虑时间变化和系统性误差的一般形式：

式中，t₁，t₂和t₃分别为下、前和后视传感器对相同地面点的成像时刻，d′_x123(t₁,t₂,t₃)为t₁，t₂和t₃时刻三线阵影像的系统性误差和颤振误差共同引起的投影偏差，κ₁，κ₂和κ₃分别为下、前和后视影像的分辨率。

从式(12)中可以看出，a₁,a₂,a₃之间，b₁,b₂,b₃之间，以及c₁,c₂,c₃之间具有相关性。因此，对系统性偏差系数进行合并，可得：

采用基于仿射模型的相对定向消除三线阵影像之间相对的系统性误差，残余的误差主要是由平台姿态颤振引起的周期性误差，则投影偏差d_x123重新表示为：

将式(14)左侧移项至右侧得到颤振估计函数模型：

式中，V_x为同名点行坐标的残差向量。

那么，目标函数为：

基于三线阵影像探测卫星平台沿轨方向姿态颤振主要包含以下主要步骤：

1)根据SIFT匹配算法，获取少量的同名特征点，采用相对定向的方法，消除三线阵影像之间相对的系统性偏差，使得同名摄影光线相交于同一物方点。

2)将下视影像作为参考影像，以N×N像素大小的空间格网点作为待匹配点，利用基于几何约束的互相关方法在三视影像上进行同名点的匹配，获取更多的同名点。

3)匹配粗差剔除。通过设定相关系数阈值剔除部分粗差点；结合相对定向参数，通过前方交会计算三线阵同名点对应的物方点坐标，并将物方点分别投影至三个影像的像方空间，计算相应的匹配点与投影点之间的像方偏差d′，统计像方偏差的平均值μ和标准差σ，按照下式剔除粗差点；

4)三线阵影像沿轨方向投影偏差计算。利用同名点对三线阵影像进行相对定向处理，结合相对定向参数，利用前、后视影像的同名点，通过前方交会计算地面点坐标，并将地面点后向投影至下视影像的像方空间，计算投影点与相应的下视匹配像点之间的偏差，即三线阵影像的系统性误差和颤振误差共同引起的视差d′_x123，利用基于仿射补偿模型的相对定向，消除三线阵影像之间相对的系统性误差，得到投影偏差d_x123；

5)三线阵影像沿轨方向姿态颤振估计，根据式(15)建立颤振估计函数模型，利用最优化估计算法，估计颤振模型中的参数，使得目标方程(16)最小。

具体的流程图如图2所示。

本发明所使用的国产某测绘卫星三线阵连续两景影像，均为0级影像数据，仅经过辐射校正处理，且未经过任何的几何校正和几何拼接，各片CCD阵列对地推扫式成像，形成各自单独的影像文件，如图3和4所示。

视影像之间的同名点，通过前方交会计算地面点坐标，并将地面点投影至下视影像的像方空间，并根据基于三线阵影像的卫星平台沿轨方向姿态颤振探测和估计理论，同时考虑到影像数据量较少(两景三线阵影像)，依据式(14)采用两个累加的谐波函数建模沿轨方向姿态颤振，并利用遗传算法估计颤振模型中的参数，使得目标方程(16)最小。A、B两景影像的颤振估计结果见表1，姿态颤振的频率约为0.6～0.7Hz之间，A景下、前和后视影像受到的最大颤振值分别为1.55，0.69和1.66像素，B景下、前和后视影像受到的最大颤振值分别为1.15，1.21和1.65像素。

三线阵立体影像沿轨方向上投影偏差受到地形起伏影响，消除地形之前，两景影像下视与前视、下视与后视之间的投影偏差如图5所示。从图中可以看出，影像之间的投影偏差受到明显的地形影响，A景下视与前视影像之间的投影偏差在-330～-10像素之间，下视与后视影像之间的偏差视差在10～330像素之间，B景下视与前视影像之间的投影偏差在-300～0像素之间，下视与后视影像之间的偏差视差在0～300像素之间。可见，地形起伏对沿轨方向的颤振探测的影响很大，故沿轨方向颤振探测必须首先消除地形。

表1 A、B景影像沿轨方向颤振估计值

在没有外部DEM参与的情况下，利用三线阵影像两两构成立体像对，抵消地形起伏的影响。根据颤振探测理论，在下视影像像方空间上，图6给出了A、B两景的原始像方投影偏差、相对定向后沿轨方向像方投影偏差和颤振补偿后像方投影偏差的平面图，投影偏差统计结果见表2。从图6(a)和(d)可以看出，A、B两景影像在无偏差补偿之前，均存在明显的系统性和周期性投影偏差，但是无任何地形信息。值得注意的是，偏差图中不存在任何地形信息，进一步验证了本文提出的基于三线阵影像的沿轨方向颤振探测方法是一种地形无关的探测方法。同时，可以看出像方偏差明显偏离0像素，说明影像中存在明显的系统性误差，另外，投影偏差图中存在条纹形状的误差，这主要是由于周期性的姿态颤振导致。理论上，同一行影像受到的颤振影像是相同的，但是在偏差图中，条纹左右方向上存在一定程度的倾斜，这主要是由于系统误差导致三线阵CCD对同一地物成像时不平行导致。周期性条纹沿着影像行方向变化，这也说明周期性的像方偏差与卫星的飞行时间相关。如表2所示，考虑地形影响后，A景和B景影像像方偏差的均方根误差分别为8.33像素和8.48像素，而受地形影响时，最大投影偏差高达350像素和300像素，这表明地形引起的视差被有效抵消。

为消除系统偏差影响，基于像方仿射补偿模型对三线阵影像进行相对定向处理，并重新计算三线阵影像之间的像方投影偏差，结果如图6(b)和(e)所示。可以看出，像方仿射模型有效的消除了三线阵影像之间的系统性偏差，使得周期性的投影偏差更加清晰明显，且处于相同影像行上的颤振值基本相等。需要特别注意的是，图中条纹的颜色沿行方向逐渐变化，说明投影偏差的峰值逐渐变小，因此实验中采用两个累加的谐波来表述变化的姿态颤振。如表2所示，相对定向后，A景和B景影像像方偏差的平均误差基本为0像素，均方根误差分别减小至1.35个像素和1.08个像素。

基于相对定向后的投影偏差，根据式(15)给出的沿轨方向像方偏差估计模型，采用遗传算法估计三线阵影像沿轨方向颤振的绝对分量。在像方空间上，利用估计的颤振值并进行颤振补偿，从图6(c)和(f)中可以看出，系统性和周期性像方偏差均被有效补偿，像方偏差几乎为0，表示立体像对之间的相对系统性偏差已经被补偿，并且周期性的条纹基本消除，反映了本文提出的颤振估计方法能够有效的补偿沿轨方向周期性偏差，并且偏差图中无任何的地形信息，再次证明基于三线阵影像的沿轨方向姿态颤振探测和估计方法，能够消除地形起伏的影响。如表2所示，颤振补偿后A景和B景影像的像方偏差均方根误差大幅减小，分别减小至0.34像素和0.32像素。

表2 A、B景三线阵影像沿轨方向像方偏差结果统计(单位：pixel)

本发明以国产某测绘卫星姿态颤振作为研究对象，基于本文提出的沿轨方向姿态颤振探测和估计方法，消除了地形起伏的影像，利用三线阵影像之间的投影偏差，消除了地形起伏的影像，解算了由于姿态颤振造成的三线阵影像周期性像方偏差。实验中利用严格传感器模型进行姿态颤振探测，结果显示在沿轨方向上存在姿态颤振，且受颤振影响较大，以2.1米的下视影像分辨率为参照标准，沿轨方向上颤振值最大为1.66像素，振幅大小随时间变化，但颤振的频率变化不大，集中在0.6～0.7Hz之间。颤振补偿后，A景影像的投影偏差从1.35像素减小至0.34像素，B景影像的投影偏差从1.08像素减小至0.32像素。

Claims (1)

1.一种高分辨率三线阵卫星沿轨方向姿态颤振的探测方法，其特征在于，包括以下步骤：

1)三线阵影像同名点匹配，并进行匹配粗差剔除，所述的步骤1)具体包括以下步骤：

11)根据SIFT匹配算法，获取初始的同名点，采用相对定向的方法，消除三视影像之间相对的系统性偏差，使得同名摄影光线相交于同一物方点；

12)将下视影像作为参考影像，以N×N像素大小的空间格网点作为待匹配点，利用基于几何约束的互相关方法在三视影像上进行同名点的匹配；

13)通过设定相关系数阈值剔除部分粗差点，结合相对定向参数，通过前方交会获取三线阵同名点对应的物方点坐标，并将物方点分别投影至三个影像的像方空间，计算相应的匹配点与投影点之间的像方偏差d′，统计像方偏差的平均值μ和标准差σ，剔除粗差点公式为：

$d^{\&prime;}=\left\{\begin{array}{cc}Inliers& (\mathrm{\&mu;}-3\mathrm{\&sigma;}<d^{\&prime;}<\mathrm{\&mu;}+3\mathrm{\&sigma;})\\ Outliers& \left(Others\right)\end{array}\right.;$

2)获取三线阵影像沿轨方向投影偏差，并建立颤振估计函数模型，所述的步骤2)具体包括以下步骤：

21)利用同名点对三线阵影像进行相对定向处理，结合相对定向参数，将前、后视影像中受到颤振影响的同名点，通过前方交会获取受到颤振影响的地面点G′的坐标，并将地面点G′后向投影至下视影像的像方空间，获得地面点坐标的投影点；

22)计算投影点与相应的下视匹配像点之间的偏差，即三线阵影像的系统性误差和颤振误差共同引起的视差d′_x123：

$d_{x123}^{\&prime;}=-{\mathrm{\&Delta;x}}_1+\frac{1}{2}(\frac{{\mathrm{\&kappa;}}_2}{{\mathrm{\&kappa;}}_1}{\mathrm{\&Delta;x}}_2+\frac{{\mathrm{\&kappa;}}_3}{{\mathrm{\&kappa;}}_1}{\mathrm{\&Delta;x}}_3)$

$\left\{\begin{array}{c}{\mathrm{\&Delta;x}}_1=a_1+b_1{x}_1+c_1{y}_1-J_{Pitch}\left(t_1\right)f\\ {\mathrm{\&Delta;x}}_2=a_2+b_2{x}_2+c_2{y}_2-J_{Pitch}\left(t_2\right)f/{\cos }^2{\mathrm{\&beta;}}_2\\ {\mathrm{\&Delta;x}}_3=a_3+b_3{x}_3+c_3{y}_3-J_{Pitch}\left(t_3\right)f/{\cos }^2{\mathrm{\&beta;}}_3\end{array}\right.$

其中，Δx₁为姿态角Pitch在下视影像上造成的像方偏差，Δx₂为姿态角Pitch在前视影像上造成的像方偏差，Δx₃为姿态角Pitch在后视影像上造成的像方偏差，κ₁为下视影像的分辨率，κ₂为前视影像的分辨率，κ₃为后视影像的分辨率，x₁和y₁为下视影像上同名点的像方坐标，x₂和y₂为前视影像上同名点的像方坐标，x₃和y₃为后视影像上同名点的像方坐标，a₁,b₁,c₁,a₂,b₂,c₂,a₃,b₃,c₃分别为下、前、后视影像系统误差的仿射模型系数，f为焦距，t₁，t₂和t₃分别为下、前、后视影像曝光时刻，β₂为前视传感器视线矢量与下视传感器视线矢量的夹角，β₃为后视传感器视线矢量与下视传感器视线矢量的夹角，J_Pitch(t)为轨道姿态角Pitch颤振函数模型，t为卫星的飞行时间，A_i为第i个谐波分量的振幅，ω_i为第i个谐波分量的角频率，为第i个谐波分量的初始相位值；

23)利用基于仿射补偿模型的相对定向，消除三线阵影像之间相对的系统性误差，得到投影偏差d_x123：

$d_{x123}(t_1,t_2,t_3)=J_{Pitch}\left(t_1\right)f-\frac{1}{2}\left(\frac{{\mathrm{\&kappa;}}_2}{{\mathrm{\&kappa;}}_1}{J}_{Pitch}\right(t_2)f/{\cos }^2{\mathrm{\&beta;}}_2+\frac{{\mathrm{\&kappa;}}_3}{{\mathrm{\&kappa;}}_1}{J}_{Pitch}(t_3)f/{\cos }^2{\mathrm{\&beta;}}_3);$

24)建立颤振估计函数模型为：

$V_x=J_{Pitch}\left(t_1\right)f-\frac{1}{2}\left(\frac{{\mathrm{\&kappa;}}_2}{{\mathrm{\&kappa;}}_1}{J}_{Pitch}\right(t_2)f/{\cos }^2{\mathrm{\&beta;}}_2+\frac{{\mathrm{\&kappa;}}_3}{{\mathrm{\&kappa;}}_1}{J}_{Phch}(t_3)f/{\cos }^2{\mathrm{\&beta;}}_3)-d_{x123}(t_1,t_2,t_3)$

其中，V_x为同名点行坐标的残差向量；

3)三线阵影像沿轨方向姿态颤振估计，所述的步骤3)具体包括以下步骤：

根据颤振估计函数模型，建立最优化目标函数，并根据最优化估计算法，估计颤振模型中的参数，使得目标函数最小，得到姿态角Pitch颤振函数模型J_Pitch(t)的参数结果，包括谐波分量的振幅、角频率和初始相位值，消除由于颤振造成的行坐标残差向量，

所述的目标函数为：

$\begin{array}{cc}min& {\mathrm{\&Sigma;V}}_x^T{V}_x\end{array}.$