CN106959454B - 一种基于数字域tdi及连续复线阵成像模式的颤振反演方法 - Google Patents

Abstract

一种基于数字域TDI及连续复线阵成像模式的颤振反演方法，其包含以下步骤：S1、对地物点P(x,y,z)进行多次拍摄，得到图像I₁(t)和I₂(t+Δt)，I₁(t)和I₂(t+Δt)的成像目标相同但二者仍存在颤振引起的像移；从中获取相移序列；S2、根据S1的结果拟合振动曲线，分别提取垂轨和沿轨方向各谐波分量A_xi(ω_i)和A_yi(ω_i)，从而反演出不同频率ω_i下最大颤振角度；S3、根据S1的结果对I₂(t+Δt)进行校正，得到I’₂(t+Δt)，对I₁(t)和I’₂(t+Δt)进行能量叠加，复原出的二级积分图像I(t)。这种方法不依赖高精度的定姿设备以及密集的地面控制点，通过分析颤振的变化规律，反演出卫星颤振的相关参数。

Description

一种基于数字域TDI及连续复线阵成像模式的颤振反演方法

技术领域

本发明涉及航天遥感遥测领域，具体涉及一种基于数字域TDI及连续复线阵成像模式的颤振反演方法。

背景技术

当探测器以推扫方式成像时，由于数字域时间延迟积分(TDI)技术的特殊工作方式，要求探测器沿TDI方向上的行周期与目标相对位置应严格同步，即探测器对目标的采样应该满足位置匹配条件，否则无法保证累积在一起的数据是像元对同一目标所成的像，导致图像模糊。卫星上活动部件的各种运动均会引起卫星的振动，相较于一般卫星，微纳遥感卫星体积小、质量轻，在轨运行稳定性不足，尚不能像大卫星一样配备执行机构的隔振平台。当振动的频率和幅值超过某一临界值就会变得发散，从而发生颤振，导致在曝光时探测器与目标的位置失配，如图1所示。因此，微纳卫星在轨摄像期间的卫星颤振反演成为提高数字域TDI成像图像质量的关键。

目前星上的星敏感器和GPS等测量部件的数据更新率在10Hz左右，无法依靠其实现更高频率的颤振检测，学者们已提出利用图像实现高频颤振的检测方法。英国剑桥大学利用面阵凝视相机采集的图像序列获取颤振参数。长光所提出了在焦平面空闲空间安装快速面阵CCD，通过一种改进的灰度投影算法计算颤振参数的方法。武汉大学进行了在轨颤动检测的可行性研究，利用小型高速面阵CCD相机实时获取与主光轴CCD同光轴或平行的图像，并与上一采样周期所得图像或标准图像进行相关运算，获得高速探测CCD多次曝光时间内光轴的偏移量。但是这两种方法需要在卫星TDICCD焦平面位置加装快速面阵CCD设备，因此不具有广泛的适用性。以色列本古里安大学利用CCD交错排列的方式获取振动参数，由于国内遥感相机焦平面很少采用这种结构因而受到了限制。日本东京大学利用光学卫星多光谱相机的特性针对TERRA卫星定量检测了110Hz以下的卫星颤振。中国科学院电子学研究所针对资源三号以及天绘一号，利用光学卫星多光谱相机以及多片CCD拼接的全色相机的特性实现了卫星在0.105Hz-4Hz的颤振检测，并根据检测结果实现了颤振图像补偿。同济大学通过不同波段的多光谱影像，以及不同视角的全色三线阵立体影像间的密集匹配精确量测和获取影像视差差异值，反演卫星姿态颤振的规律和特性。针对资源三号实现了0.65Hz的颤振探测。综上所述，在不加装快速面阵CCD设备的情况下，现有的方法针对多光谱相机、全色三线阵立体相机以及多片拼接的全色相机实现了卫星颤振检测，并且只在110Hz以下进行了验证，更高频率的颤振探测还处于空白区域。

成像时刻的载荷高频颤振主要是由于星上执行部件微动造成的周期性姿态不稳，由于数字域TDI是多个线阵分时成像的，对同一目标多次曝光，分时成像，因此这种不稳定性在每个线阵的重叠区域会有相应的错位表现。

发明内容

本发明的目的在于提供一种基于数字域TDI及连续复线阵成像模式的颤振反演方法，适用于全色单线阵相机影像的高频颤振，这种方法不依赖高精度的定姿设备以及密集地面控制点，可以直接从影像上进行1000Hz以上的颤振检测，通过分析颤振的变化规律最终可以反演出卫星颤振的相关参数。

一种基于数字域TDI及连续复线阵成像模式的颤振反演方法，其特征是，包含以下步骤：

S1、对地物点P(x,y,z)进行拍摄，以间隔时间为Δt的两行单级积分图像为一组，反复多次向前推进，从而累加成行数为M的图像I₁(t)和I₂(t+Δt)，I₁(t)和I₂(t+Δt)的成像目标相同但二者仍存在颤振引起的像移；从中获取相移序列：

其中(x,y,z)为地物点在大地测量坐标系中的三维坐标；u和u’分别为地物点在单级积分图像I₁和I₂上的横坐标，v和v’分别为地物点在单级积分图像I₁和I₂上的纵坐标，0<u≤N，0<v≤M，0<u’≤N，0<v’≤M，M、N分别为上述图像的行数和列数；

S2、根据S1的结果拟合振动曲线，分别提取垂轨和沿轨方向各谐波分量A_xi(ω_i)和A_yi(ω_i)，从而反演出不同频率ω_i下最大颤振角度：

其中，A_xi和A_yi分别为垂轨和沿轨方向第i个谐波分量的幅值，θ_pi为频率ω_i下俯仰轴角振幅，θ_ri为频率ω_i下滚动轴角振幅，f为相机焦距；

S3、根据S1的结果对I₂(t+Δt)进行校正，得到I’₂(t+Δt)，对I₁(t)和I’₂(t+Δt)进行能量叠加，复原出的二级积分图像：

I(t)＝I₁(t)+I’₂(t+Δt)

至此，整幅图像全部处理完毕，得到去除颤振的清晰图像。

上述的一种基于数字域TDI及连续复线阵成像模式的颤振反演方法，其中，所述的步骤S1具体包含：

S11、令第一线阵成像开始时间为t_s，第一线阵在t_s时刻输出的单级积分图像为I₁(t_s)，在经过时间Δt后，第一线阵在t_s+Δt时刻继续输出第二行单级积分图像I₁(t_s+Δt)；第二线阵在t_s+Δt时刻输出第一行单级积分图像I₂(t_s+Δt)；

S12、令t＝t_s+k·Δt，其中k为当前图像累积的行数，0<(k+1)≤M，令k逐渐递增，逐行累积图像I₁(t)，得到累积了(k+1)行的图像，尺寸为(k+1)·N；逐行累积图像I₂(t+Δt)，得到累积了k行的图像，尺寸为k·N；

S13、利用特征提取算子，对图像I₁和I₂从第1行至第M行逐行进行同名点匹配，实现像素级的像元匹配；

S14、在S13实现像素级同名点匹配的基础上，对图像I₁和I₂进行亚像元级同名点匹配，Δt时间内产生的像移为(Δu,Δv)，则：

M_x(u,v)和M_y(u,v)分别为[t,t+Δt](t＝t_s+v·Δt)时间段内垂轨方向和沿轨方向的像移量序列；

S15、利用平滑滤波对垂轨方向像移序列M_x(u,v)和沿轨方向像移序列M_y(u,v)剔除粗差和飞点。

上述的一种基于数字域TDI及连续复线阵成像模式的颤振反演方法，其中，所述的步骤S2具体包含：

S21、将像移序列减去第一线阵和第二线阵之间的固定位移差，得到颤振像移序列；

S22、通过垂轨方向颤振像移序列和沿轨方向颤振像移序列中的大量散点，进行振动曲线拟合；

S23、分别对垂轨、沿轨两个方向的离散颤振像移序列进行离散傅里叶变换，计算其频谱和功率谱密度，最终提取颤振信号各谐波分量的幅值A_xi、A_yi和相对应谐波分量的频率ω_i，其中，i代表第i个谐波分量；

S24、利用提取出的颤振各谐波分量的频率ω_i和幅值A_xi、A_yi反演出不同频率下的最大颤振角度：

其中θ_ri为频率ω_i下滚动轴角振幅、θ_pi为频率ω_i下俯仰轴角振幅，f为相机焦距。

上述的一种基于数字域TDI及连续复线阵成像模式的颤振反演方法，其中，所述的步骤S3具体包含：

S31、根据S1检测出的帧间像移量，对前后两帧互相关图像进行像移补偿，以I₁(t)为基准行图像，结合M_x(u,v)和M_y(u,v)对I₂(t+Δt)逐行进行像移校正，得到I’₂(t+Δt)：

I’₂(t+Δt)＝F(I₂(t+Δt))

其中，F(·)为像移校正函数；

S32、逐行对I₁(t)和I’₂(t+Δt)进行能量累加，复原出消除颤振混迭的能量更强的二级积分图像I(t)：

I(t)＝I₁(t)+I’₂(t+Δt)

至此，整幅图像全部处理完毕，得到去除颤振的清晰图像。

本发明提供了一种在数字域TDI及连续复线阵的成像模式下的颤振检测方法，该方法通过对分时成像的互相关全色影像进行逐行逐列的亚像元匹配，从而逐像元地求取相邻阵列之间的行、列方向偏移，最后减去相邻线阵间的固定偏移量即可定量地反算出每个像元的颤振偏移量。

附图说明

图1是颤振引起的探测器行间错位示意图。

图2是数字域TDI及连续复线阵成像模式示意图。

图3是基于数字域TDI及连续复线阵成像模式的颤振反演流程图。

具体实施方式

以下结合附图，通过详细说明一个较佳的具体实施例，对本发明做进一步阐述。

如图2所示，星上照相机在数字域TDI成像模式下对同一地标以固定的间隔时间Δt连续拍摄多幅图像，多行累加后得到如图2中图像3所示的多级积分图像。数字域TDI及连续复线阵的成像模式可以分为两类功能区域：TDI技术成像阵列区域(如图2中图像3及以上部分所示)和复线阵成像区域(如图2中图像1和图像2所示)。TDI技术成像阵列区域输出一幅多级积分的目标图像，各行像素依次标记为第一行像素p₁，第二行像素p₂，……，第j行像素p_j。在经过时间t_i后，各行图像分别输出多行累加后的多次曝光图像，例如，第j行像素在t_i时刻输出的图像I_3j(t_i)＝p_j(t_i)+p_j-1(t_i-1)……+p₁(t_i-j+1)。但该图像在累积多级能量的同时也会累积0～t_i时间段内的颤振信息而导致图像模糊，本发明对这一部分图像不进行处理。复线阵成像区域由2个独立成像的线阵组成，各行像素依次标记为q₁和q₂，分别输出只含一级颤振信息的单行图像I₁(t_i)＝q₁(t_i)和I₂(t_i)＝q₂(t_i)。本发明利用复线阵成像区域的图像实现颤振反演。

S1、对地物点P(x,y,z)进行拍摄，以间隔时间为Δt的两行单级积分图像为一组，反复多次向前推进，从而累加成行数为M的图像I₁(t)和I₂(t+Δt)，I₁(t)和I₂(t+Δt)的成像目标相同但二者仍存在颤振引起的像移；从中获取相移序列：

其中(x,y,z)为地物点在大地测量坐标系中的三维坐标；u和u’分别为地物点在单级积分图像I₁和I₂上的横坐标，v和v’分别为地物点在单级积分图像I₁和I₂上的纵坐标，0<u≤N，0<v≤M，0<u’≤N，0<v’≤M，M、N分别为上述图像的行数和列数。

具体可分为如下5步：

S11、令第一线阵成像开始时间为t_s，第一线阵在t_s时刻输出的单级积分图像为I₁(t_s)，此图像尺寸为1·N，N为探测器的列数(每一个线阵在单级积分时间内只能输出一行图像)，即I₁(t_s)为I₁(t)累积的第一行图像，如图2中图像1所示。在经过时间Δt后，第一线阵在t_s+Δt时刻继续输出第二行单级积分图像I₁(t_s+Δt)；第二线阵在t_s+Δt时刻输出第一行单级积分图像I₂(t_s+Δt)，I₂(t_s+Δt)为I₂(t+Δt)累积的第一行图像，如图2中图像2所示，此图像尺寸亦为1·N。I₁(t_s)与I₂(t_s+Δt)之间在列与行方向上均有重叠。

S12、令t＝t_s+k·Δt，其中k为当前图像累积的行数，0<(k+1)≤M，令k逐渐递增，逐行累积图像I₁(t)，得到累积了(k+1)行的图像，尺寸为(k+1)·N；逐行累积图像I₂(t+Δt)，得到累积了k行的图像，尺寸为k·N。设该系统Δt＝280μs，为保证频谱分析的频率分辨率Δf<1，根据Δf＝1/(Δt·M)，可知M>1/Δt，取M>3571。

S13、利用特征提取算子，对图像I₁和I₂从第1行至第M行逐行进行同名点匹配，实现像素级的像元匹配。

S14、在S13实现像素级同名点匹配的基础上，对图像I₁和I₂进行亚像元级同名点匹配。假设第一线阵在t时刻对地物点P(x,y,z)成像后，P(x,y,z)在图像I₁上的坐标为(u,v)，经过时间Δt后，第二线阵对同一地物点P(x,y,z)曝光成像，P(x,y,z)在图像I₂上的坐标为(u’,v’)，Δt时间内产生的像移为(Δu,Δv)，则：

其中(x,y,z)为地物点在大地测量坐标系的三维坐标；u、u’分别为图像I₁(t)和图像I₂(t+Δt)上的横坐标，v、v’分别为图像I₁(t)和图像I₂(t+Δt)上的纵坐标；N＝3072、M＝3571；0<u≤3072，0<v≤3571；0<u’≤3072，0<v’≤3571；M_x(u,v)和M_y(u,v)分别为[t,t+Δt](t＝t_s+v·Δt)时间段内垂轨方向和沿轨方向的像移量序列。

S15、利用平滑滤波对垂轨方向像移序列M_x(u,v)和沿轨方向像移序列M_y(u,v)剔除粗差和飞点。

S2、根据S1的结果式(1)拟合振动曲线，分别提取垂轨和沿轨方向各谐波分量A_xi(ω_i)和A_yi(ω_i)，从而反演出不同频率ω_i下最大颤振角度：

其中，A_xi和A_yi分别为垂轨和沿轨方向第i个谐波分量的幅值，θ_pi为频率ω_i下俯仰轴角振幅，θ_ri为频率ω_i下滚动轴角振幅，f为相机焦距。

平台的颤振可分为滚动轴、俯仰轴和偏航轴三个方向的振动。下面就三个方向分别进行分析。

设星体坐标系xyz，θ、ψ分别为滚动角、俯仰角和偏航角。设θ_r、θ_p、θ_y分别为三个方向上的最大角振幅，相机焦距为f。

滚动轴方向

绕滚动轴方向的颤振会在垂直于卫星运动方向，即y方向上产生像移。角振幅为θ_r，则像移的振幅为A_y1＝f×tanθ_r。

俯仰轴方向

俯仰轴方向的颤振会在沿卫星运动方向，即x方向产生像移。若角振幅为θ_p，则像移的振幅为A_x1＝f×tanθ_p。

偏航轴方向

偏航轴方向的颤振会同时在运动方向和垂直方向，即x方向和y方向上产生像移。若CCD像元大小为a，角振幅为θ_y，则x方向的像移振幅A_x2＝a×sinθ_y，y方向的像移振幅为A_y2＝a×(1-cosθ_y)。

根据上述分析，当卫星的三个轴同时存在同一个频率的颤振时，在一个积分时间内，垂直TDI积分方向的像移振幅应为：A_y＝A_y1+A_y2，沿TDI积分方向的像移振幅为：A_x＝A_x1+A_x2。综上可得像移振幅公式：

A_x＝ftanθ_p+asinθ_y(5)

A_y＝ftanθ_r+a(1-cosθ_y)(6)

然而在实际应用中，相机的焦距f一般是米量级，CCD像元大小a一般是微米量级，两者相差近6个数量级，而颤振引起的各个轴的角度振幅一般在同一量级，由此可见由偏航角引起的像移可以忽略不计，因此可以将式(5)近似为：

A_x＝ftanθ_p (7)

同理，将式(6)近似为：

A_y＝ftanθ_r (8)

进一步，由于颤振引起的各个轴的角度振幅一般在角秒级，因此tanθ≈θ，将式(7)和(8)近似为：

A_x＝f×θ_p (9)

A_y＝f×θ_r (10)

由此可得到：

θ_p＝A_x/f (11)

θ_r＝A_y/f (12)

通过上述分析，当通过对像移量的频谱分析提取出的颤振分量的频率和振幅后，即可以通过上述像移振幅和角度的关系式(11)和(12)反演出不同频率的颤振角度。

具体过程可分为如下4步：

S21、将像移序列减去第一线阵和第二线阵之间的固定位移差后，即可得到颤振像移序列。本实施例中第一线阵和第二线阵之间在垂轨方向的位移可视为零，因此仅需将沿轨方向像移序列M_y(u,v)减去第一线阵和第二线阵之间的固定行位移差后，即可得到颤振像移序列M_y’(u,v)。

S22、通过垂轨方向像移序列M_x(u,v)和沿轨方向像移序列M_y’(u,v)中的大量散点，进行振动曲线拟合。

S23、分别对离散序列M_x(u,v)和M_y’(u,v)进行离散傅立叶变换，计算其频谱和功率密度谱，最终提取颤振信号各谐波分量的幅值A_xi、A_yi和相对应谐波分量的频率ω_i。其中，i代表第i个谐波分量，A_xi为垂轨方向第i个谐波分量的幅值，A_yi为沿轨方向第i个谐波分量的幅值。

根据奈奎斯特采样定理，可测的最大带宽为采样频率的一半，设该系统时间差Δt＝280μs，则测量带宽为1/2×Δt＝1785Hz。即可探测到1785Hz以下的各个频率的颤振信息。

S24、利用提取出的颤振各谐波分量的频率ω_i和幅值A_xi、A_yi，反演出不同频率ω_i下的最大颤振角度。

其中θ_ri为频率ω_i下滚动轴角振幅、θ_pi为频率ω_i下俯仰轴角振幅，f为相机焦距。

S3、根据S1的结果对I₂(t+Δt)进行校正，得到I’₂(t+Δt)，对I₁(t)和I’₂(t+Δt)进行能量叠加，复原出二级积分图像：

I(t)＝I₁(t)+I’₂(t+Δt) (3)。

具体过程可分为如下2步：

S31、根据S1检测出的帧间像移量，对前后两帧互相关图像进行像移补偿。以I₁(t)为基准行图像，结合M_x(u,v)和M_y(u,v)对I₂(t+Δt)逐行进行像移校正，得到I’₂(t+Δt)，

I’₂(t+Δt)＝F(I₂(t+Δt)) (4)

其中，F(·)为像移校正函数。

S32、逐行对I₁(t)和I’₂(t+Δt)进行能量累加，复原出消除颤振混迭的能量更强的二级积分图像I(t)：

I(t)＝I₁(t)+I’₂(t+Δt) (3)。

至此，整幅图像全部处理完毕，得到去除颤振的清晰图像。

尽管本发明的内容已经通过上述优选实施例作了详细介绍，但应当认识到上述的描述不应被认为是对本发明的限制。在本领域技术人员阅读了上述内容后，对于本发明的多种修改和替代都将是显而易见的。因此，本发明的保护范围应由所附的权利要求来限定。

Claims (4)

1.一种基于数字域TDI及连续复线阵成像模式的颤振反演方法，其特征在于，包含以下步骤：

S1、对地物点P(x,y,z)进行拍摄，以间隔时间为Δt的两行单级积分图像为一组，反复多次向前推进，从而累加成行数为M的图像I₁(t)和I₂(t+Δt)，I₁(t)和I₂(t+Δt)的成像目标相同但二者仍存在颤振引起的像移；从中获取相移序列：

其中(x,y,z)为地物点在大地测量坐标系中的三维坐标；u和u’分别为地物点在单级积分图像I₁和I₂上的横坐标，v和v’分别为地物点在单级积分图像I₁和I₂上的纵坐标，0＜u≤N，0＜v≤M，0＜u’≤N，0＜v’≤M，N为单级积分图像I₁和I₂的列数；

S2、根据S1的结果拟合振动曲线，分别提取垂轨和沿轨方向各谐波分量A_xi(ω_i)和A_yi(ω_i)，从而反演出不同频率ω_i下最大颤振角度：

其中，A_xi和A_yi分别为垂轨和沿轨方向第i个谐波分量的幅值，θ_pi为频率ω_i下俯仰轴角振幅，θ_ri为频率ω_i下滚动轴角振幅，f为相机焦距；

S3、根据S1的结果对I₂(t+Δt)进行校正，得到I’₂(t+Δt)，对I₁(t)和I’₂(t+Δt)进行能量叠加，复原出二级积分图像：

I(t)＝I₁(t)+I’₂(t+Δt) (3)。

2.如权利要求1所述的一种基于数字域TDI及连续复线阵成像模式的颤振反演方法，其特征在于，所述的步骤S1具体包含：

S11、令第一线阵成像开始时间为t_s，第一线阵在t_s时刻输出的单级积分图像为I₁(t_s)，在经过时间Δt后，第一线阵在t_s+Δt时刻继续输出第二行单级积分图像I₁(t_s+Δt)；第二线阵在t_s+Δt时刻输出第一行单级积分图像I₂(t_s+Δt)；

S12、令t＝t_s+k·Δt，其中k为当前图像累积的行数，0＜(k+1)≤M，令k逐渐递增，逐行累积图像I₁(t)，得到累积了(k+1)行的图像，尺寸为(k+1)·N；逐行累积图像I₂(t+Δt)，得到累积了k行的图像，尺寸为k·N；

S13、利用特征提取算子，对图像I₁(t)和I₂(t+Δt)从第1行至第M行逐行进行同名点匹配，实现像素级的像元匹配；

S14、在S13实现像素级同名点匹配的基础上，对图像I₁(t)和I₂(t+Δt)进行亚像元级同名点匹配，Δt时间内产生的像移为(Δu,Δv)，则整幅图的移相序列为：

M_x(u,v)和M_y(u,v)分别为[t,t+Δt](t＝t_s+v·Δt)时间段内垂轨方向和沿轨方向的像移量序列；

S15、利用平滑滤波对垂轨方向像移序列M_x(u,v)和沿轨方向像移序列M_y(u,v)剔除粗差和飞点。

3.如权利要求1所述的一种基于数字域TDI及连续复线阵成像模式的颤振反演方法，其特征在于，所述的步骤S2具体包含：

S21、将像移序列减去第一线阵和第二线阵之间的固定位移差，得到颤振像移序列；

S22、通过垂轨方向颤振像移序列和沿轨方向颤振像移序列中的大量散点，进行振动曲线拟合；

S23、分别对垂轨、沿轨两个方向的离散颤振像移序列进行离散傅里叶变换，计算其频谱和功率谱密度，最终提取颤振信号各谐波分量的幅值A_xi、A_yi和相对应谐波分量的频率ω_i，其中，i代表第i个谐波分量；

S24、利用提取出的颤振各谐波分量的频率ω_i和幅值A_xi、A_yi反演出不同频率下的最大颤振角度：

其中θ_ri为频率ω_i下滚动轴角振幅、θ_pi为频率ω_i下俯仰轴角振幅，f为相机焦距。

4.如权利要求1所述的一种基于数字域TDI及连续复线阵成像模式的颤振反演方法，其特征在于，所述的步骤S3具体包含：

S31、根据S1检测出的帧间像移量，对前后两帧互相关图像进行像移补偿，以I₁(t)为基准行图像，结合M_x(u,v)和M_y(u,v)对I₂(t+Δt)逐行进行像移校正，得到I’₂(t+Δt)：

I’₂(t+Δt)＝F(I₂(t+Δt)) (4)

其中，F(·)为像移校正函数；

S32、逐行对I₁(t)和I’₂(t+Δt)进行能量累加，复原出消除颤振混迭的能量更强的二级积分图像I(t)：

I(t)＝I₁(t)+I’₂(t+Δt) (3)。