CN107966137B - 一种基于tdiccd拼接区图像的卫星平台颤振探测方法 - Google Patents

Abstract

本发明公开了一种TDICCD拼接区图像的卫星平台颤振探测方法，包括以下步骤：1)重叠图像的获取，利用时间延迟积分电荷耦合器件拼接区获得重叠图像。2)相对成像位置差的计算，对两幅重叠图像进行精确密集匹配处理，获得同名匹配点，计算同一目标在两幅重叠图像中的相对成像位置差；3)卫星平台颤振的估计，根据相对成像位置差，估计卫星平台颤振。本发明能有效解决颤振探测准确度低、仅能探测若干孤立频率点的问题，为卫星平台提供一种高精度、较宽频段内全覆盖的颤振探测方法，从而提高卫星姿态探测能力和遥感图像信息提取的准确度。

Description

一种基于TDICCD拼接区图像的卫星平台颤振探测方法

技术领域

本发明涉及航天遥感成像领域，具体涉及一种基于时间延迟积分电荷耦合器件(Time-delayed and Integration Charge Coupled Device TDICCD)拼接区图像的卫星平台颤振探测方法。

背景技术

卫星姿态信息对于在轨运行阶段的卫星姿态稳定控制和地面后期图像信息的提取都至关重要，颤振作为卫星运行过程中的普遍现象，对卫星姿态数据的精确度有着直接的影响。在高性能卫星姿态传感器的支撑下，欧美等遥感科技强国能够对一定范围内的颤振进行有效的探测，从而获得高精度的卫星姿态数据，然而，目前国内的传感技术尚不成熟，且国外在高性能姿态传感器市场进行了销售封锁，我国当前采用的星载姿态传感器存在测量带宽低、精度差等问题，无法对卫星平台颤振进行有效的测量，制约了我国高分辨率空间相机的成像质量和后期图像信息提取的准确度。如何在不依赖于高精度卫星姿态传感设备的情况下，对高分辨率遥感卫星平台颤振进行有效探测，成为我国空间遥感成像领域亟待解决的关键问题之一。

近年来，基于遥感图像的颤振探测方法受到国内外关注，作为颤振探测的首要条件，具有重叠区域的遥感图像的获取成为国内外相关学者关注的焦点。英国剑桥大学利用面阵凝视相机获得重叠图像，长光所和武汉大学提出在焦面上增加辅助面阵成像传感器的方式拍摄重叠图像，日本东京大学和同济大学利用多光谱相机不同谱段的成像传感器获得重叠景物图像。虽然重叠图像获取方式多种多样，但颤振探测算法却大致相同，通常采用频谱分析技术求解颤振的幅值、频率和相位，进而利用三角函数拟合颤振曲线，该算法能够对能量集中的孤立频率点处的主要颤振分量进行提取，而对于能量较低、占据频带较宽的次要颤振分量则较难探测，此外，由于频谱分析技术中难以避免的能量泄漏现象引起的颤振计算误差是不能忽视的。因此，如何解决现有技术中颤振探测准确度低、仅能探测若干孤立频率点的问题，提出一种全新的高精度、较宽频段内全覆盖的颤振探测方法，成为当前我国高分辨率遥感卫星的迫切任务。

发明内容

本发明的目的是针对上述现有技术存在的缺陷，提出一种基于TDICCD拼接区图像的卫星平台颤振探测方法，以期能有效解决颤振探测准确度低、仅能探测若干孤立频率点的问题，为卫星平台提供一种高精度、较宽频段内全覆盖的颤振探测方法，从而提高卫星姿态探测能力和遥感图像信息提取的准确度。

本发明为解决技术问题采用如下技术方案：

本发明一种基于TDICCD拼接区图像的卫星平台颤振探测方法的特点是按如下步骤进行：

步骤1、重叠图像的获取：

步骤1.1、在卫星平台上搭载有高分辨率空间相机，所述高分辨率空间相机的焦平面上存在多片交错拼接排布的TDICCD，任意选取两片相邻的TDICCD，记为TDICCD₁和TDICCD₂，在沿TDI方向上，TDICCD₁和TDICCD₂的间隔行数记为L，在垂直于TDI方向上，TDICCD₁和TDICCD₂的交接处存在若干列重叠像元形成拼接区，记为OA₁和OA₂，所述拼接区OA₁和OA₂内的重叠像元列数记为C；

所述高分辨率空间相机任意一次拍摄任务的起始时刻记为0时刻，拍摄持续时长记为T_w，其中T_w∈Q，T_w＞0，Q表示有理数集合；

利用所述拼接区OA₁和OA₂在[0,T_w]期间对探测区域内的所有目标进行两次分时拍摄，获得图像I₁＝(I₁ ¹,I₁ ²,…,I₁ ^r,…,I₁ ^R)^T和I₂＝(I₂ ¹,I₂ ²,…,I₂ ^r,…,I₂ ^R)^T，其中为所述图像I₁和I₂的总行数，T_r为每行图像的曝光时长，I₁ ^r和I₂ ^r分别为所述图像I₁和I₂的第r行图像，T_r∈Q，T_r＞0，1≤r≤R；

步骤1.2、提取所述图像I₁和I₂中的重叠部分(I₁ ¹,I₁ ²,…,I₁ ^R-L)^T和(I₂ ^L+1,I₂ ^L+2,…,I₂ ^R)^T，分别记为重叠图像O₁和O₂，所述重叠图像O₁和O₂的总列数为C，总行数为记为M＝R-L；

步骤2、相对成像位置差的计算：

步骤2.1、利用图像配准方法对所述重叠图像O₁和O₂进行逐行匹配处理，获得同名配准点集合P＝{(p₁ ¹,p₂ ¹),(p₁ ²,p₂ ²),…,(p₁ ^m,p₂ ^m),…,(p₁ ^M,p₂ ^M)}，其中(p₁ ^m,p₂ ^m)为所述重叠图像O₁和O₂的第m行图像配准后得到的一对同名配准点，1≤m≤M，并有： 和分别为点p₁ ^m在所述重叠图像O₁内的行坐标和列坐标，和分别为点p₂ ^m在所述重叠图像O₂内的行坐标和列坐标，

步骤2.2、分别计算所述同名配准点集合P中每对同名配准点内两点的行坐标差值和列坐标差值，获得沿TDI方向上的相对成像位置差集合S^||＝{s^||(1),s^||(2)，…,s^||(m)，…,s^||(M)}和垂直于TDI方向上的相对成像位置差集合S^⊥＝{s^⊥(1),s^⊥(2)，…,s^⊥(m)，…,s^⊥(M)}，其中s^||(m)为所述重叠图像O₁和O₂的第m行图像配准后得到的沿TDI方向上的相对成像位置差，且s^⊥(m)为所述重叠图像O₁和O₂的第m行图像配准后得到的垂直于TDI方向上的相对成像位置差，且

步骤2.3、将所述沿TDI方向上的相对成像位置差集合S^||从第一个元素开始每连续的L个元素排成一行，从而形成N行×L列的二维矩阵S_∑ ^||：

其中

步骤2.4、将所述垂直于TDI方向上的相对成像位置差集合S^⊥从第一个元素开始每连续的L个元素排成一行，从而形成N行×L列的二维矩阵S_∑ ^⊥：

其中

步骤3、卫星平台颤振的估计：

步骤3.1、定义沿TDI方向上的颤振为：

其中，g^||(v·L+u)为(v·L+u)×T_r时刻沿TDI方向的颤振，g_||(v·L+u)∈Q，u∈Z，v∈Z，1≤u≤L，0≤v≤N，Z表示整数集合；

定义垂直于TDI方向上的颤振为：

其中，g^⊥(v·L+u)为(v·L+u)×T_r时刻垂直于TDI方向的颤振，g_⊥(v·L+u)∈Q，u∈Z,v∈Z，1≤u≤L，0≤v≤N；

定义颤振为G_∑＝{g_ij}_(N+1)×L，其中 为沿TDI方向的单位矢量，为垂直于TDI方向的单位矢量，1≤i≤N+1，1≤j≤L；

定义沿TDI方向上颤振G^∑||中的第一行元素[g^||(1),g^||(2),…,g^||(L)]为沿TDI方向上的颤振子块G₁ ^||；

定义垂直于TDI方向上颤振G_∑ ^⊥中的第一行元素[g^⊥(1),g^⊥(2),…,g^⊥(L)]为垂直于TDI方向上的颤振子块G₁ ^⊥；

定义沿TDI方向上颤振G_∑ ^||中的第二行元素至第N+1行元素

为沿TDI方向上的颤振子块G₂ ^||；

定义垂直于TDI方向上颤振G_∑ ^⊥中的第二行元素至第N+1行元素

为垂直于TDI方向上的颤振子块G₂ ^⊥；

步骤3.2、将所述拍摄任务中[T_r,L·T_r]时间段内的卫星在轨运行参数带入星下点对地成像模型中，分别得到沿TDI方向上颤振子块G₁ ^||在1Hz以下的低频分量和垂直于TDI方向上颤振子块G₁ ^⊥在1Hz以下的低频分量其中表示u×T_r时刻沿TDI方向上颤振在1Hz以下的低频分量，表示u×T_r时刻垂直于TDI方向上颤振在1Hz以下的低频分量，

步骤3.3、以所述沿TDI方向上颤振子块G₁ ^||与其低频分量之间的相对残差的平方和作为沿TDI方向上的目标函数H(G₁ ^||)，如式(1)所示：

式(1)中，g^||(i)表示i×T_r时刻沿TDI方向上的颤振，表示j×T_r时刻沿TDI方向上颤振在1Hz以下的低频分量；

以所述垂直于TDI方向上颤振子块G₁ ^⊥与其低频分量之间的相对残差的平方和作为垂直于TDI方向上目标函数H(G₁ ^⊥)，如式(2)所示：

式(2)中，g^⊥(i)表示i×T_r时刻垂直于TDI方向上的颤振，表示j×T_r时刻垂直于TDI方向上颤振在1Hz以下的低频分量；

步骤3.4、利用优化算法分别求解所述目标函数H(G₁ ^||)和H(G₁ ^⊥)的最小值点和并分别作为所述沿TDI方向上颤振子块G₁ ^||和垂直于TDI方向上颤振子块G₁ ^⊥的估计值，如式(3)和式(4)所示：

步骤3.5、根据空间相机的推扫成像原理和TDICCD拼接区成像特点，建立如式(5)和式(6)所示的沿TDI方向和垂直于TDI方向的颤振探测模型：

G₂ ^||＝A·G₁ ^||+B·S_∑ ^|| (5)

G₂ ^⊥＝A·G₁ ^⊥+B·S_∑ ^⊥ (6)

式(5)和式(6)中，A和B均表示系数矩阵，且

步骤3.6、根据式(5)所示沿TDI方向的颤振探测模型，利用式(7)计算沿TDI方向上颤振子块G₂ ^||的估计值

根据式(6)所示垂直于TDI方向的颤振探测模型，利用式(8)计算垂直于TDI方向上颤振子块G₂ ^⊥的估计值

步骤3.7、根据沿TDI方向上颤振的定义，利用式(9)合成沿TDI方向上的颤振G_∑ ^||的估计值

根据垂直于TDI方向上颤振的定义，利用式(10)合成垂直于TDI方向上的颤振G_∑ ^⊥的估计值

与现有技术相比，本发明的有益效果在于：

1、本发明针对现有技术颤振探测准确度低、仅能探测若干孤立频率点的问题，提出一种基于TDICCD拼接区图像的卫星平台颤振探测方法，实现了在不依赖高精度姿态传感设备和辅助成像传感器的情况下，对卫星平台颤振进行高精度、较宽频段内的全覆盖探测。

2、卫星平台颤振频率成分复杂，现有方法利用若干个三角函数的叠加拟合颤振，能够提取能量集中的颤振主分量，往往是若干个孤立的频率点，而对能量较低且占据大段连续频段的次要分量则难以探测，导致颤振探测结果与真实情况之间存在较大差距。本发明通过建立颤振探测模型，明确了颤振与相对成像位置差之间的定量关系，避免了现有技术对颤振近似拟合引入的理论误差，是提高颤振探测准确度的有力支撑。

3、现有方法通常依赖频谱分析技术得到三角函数的频率、幅值和相位，而频谱分析技术中难以避免的能量泄漏现象会引入显著的频率、幅值和相位估计误差，最终导致的颤振探测误差是不容忽视的。本发明基于颤振探测模型，采用时域分析技术，直接得到拍摄时刻卫星平台颤振的时域值，避免了现有技术中能量泄漏引入的颤振探测误差，进一步提高了颤振探测的准确度。

4、卫星上不同位置处的颤振有所差别，即使同一位置处的颤振也是随时间不断变化的，而空间相机焦面上主成像传感器处、图像拍摄时刻的颤振是人们最为关心的。对于非凝视成像且非多光谱相机而言，现有技术往往通过在焦面上主成像传感器附近增加辅助的面阵成像传感器来探测颤振，探测结果为焦面上辅助成像传感器处而非主成像传感器处的颤振，两者之间存在差异，尤其对于相对孔径较小的空间相机而言，焦平面上不同位置处的颤振差异更为显著。另外，两类成像传感器的工作方式不同无法同步曝光，上述两方面因素最终导致颤振探测结果与所需探测的位置和时刻均有偏差。本发明通过提取主成像传感器的拼接区所拍图像用于颤振探测，探测结果为主成像传感器处、图像拍摄时刻的颤振，克服了现有技术中颤振探测结果与所需探测的位置和时刻存在偏差的问题，为进一步提高颤振探测的准确度奠定了基础。

5、卫星平台颤振占据的频带较宽，现有方法可以探测到能量集中的孤立频率点处的颤振，而尚未实现较宽频段内的全覆盖探测。本发明能够通过提高图像配准的密集度，提高了颤振的采样频率，进而提高了颤振探测的带宽，颤振探测带宽最高可以提高至为相机每行图像曝光时长倒数的二分之一。

附图说明

图1为本发明方法流程框图；

图2为本发明等间距均匀取点方式示意图；

图3为本发明粗配准和精配准相结合的亚像元配准流程图；

图4为现有技术中相邻两片TDICCD拼接区成像示意图；

图5为现有技术中某高分辨率TDICCD可见光相机焦面排布示意图；

图6为现有技术中某高分辨率TDICCD可见光相机在轨拍摄的可见光全色遥感图像；

图7a为本发明垂直于TDI方向的卫星平台颤振探测结果图；

图7b为本发明沿TDI方向的卫星平台颤振探测结果图；

图7c为本发明垂直于TDI方向卫星平台颤振探测结果的功率谱密度图；

图7d为本发明沿TDI方向卫星平台颤振探测结果的功率谱密度图。

具体实施方式

本实施例中，一种基于TDICCD拼接区图像的卫星平台颤振探测方法，如图1所示，按如下步骤进行：

步骤1、重叠图像的获取：

步骤1.1、在卫星平台上搭载有高分辨率空间相机，高分辨率空间相机的焦平面上通常存在多片交错拼接排布TDICCD，任意选取两片相邻的TDICCD，记为TDICCD₁和TDICCD₂，在沿TDI方向上，TDICCD₁和TDICCD₂的间隔行数记为L，在垂直于TDI方向上，TDICCD₁和TDICCD₂的交接处存在若干列重叠像元形成拼接区，记为OA₁和OA₂，拼接区OA₁和OA₂内的重叠像元列数记为C；

高分辨率空间相机任意一次拍摄任务的起始时刻记为0时刻，拍摄持续时长记为T_w，其中T_w∈Q,T_w＞0；Q表示有理数集合；

利用拼接区OA₁和OA₂在[0,T_w]期间对其探测区域内的所有目标进行两次分时拍摄，获得图像I₁＝(I₁ ¹,I₁ ²,…,I₁ ^r,…,I₁ ^R)^T和I₂＝(I₂ ¹,I₂ ²,…,I₂ ^r,…,I₂ ^R)^T，其中为图像I₁和I₂的总行数，T_r为每行图像的曝光时长，I₁ ^r和I₂ ^r分别为图像I₁和I₂的第r行图像，T_r∈Q，T_r＞0，1≤r≤R；

步骤1.2、提取图像I₁和I₂中的重叠部分(I₁ ¹,I₁ ²,…,I₁ ^R-L)^T和(I₂ ^L+1,I₂ ^L+2,…,I₂ ^R)^T，分别记为重叠图像O₁和O₂，重叠图像O₁和O₂的总列数为C，总行数为(R-L)，记M＝R-L；

步骤2、相对成像位置差的计算：

步骤2.1、利用图像配准方法对重叠图像O₁和O₂进行逐行匹配处理。逐行匹配处理具体过程如图2所示，在重叠图像O₁的每行图像内，以等间距均匀方式选取若干个配准点，以每个配准点为中心选取一定大小的模板作为参考模板，选取重叠图像O₂内对应位置处的模板作为待配准模板，采用图像配准方法对两幅模板进行匹配处理，获得一对同名配准点。对同行图像内的其他配准点进行同样处理，并将同行内所有配准点的配准结果取算术平均值，作为整行图像的配准结果。

本发明采用粗配准和精配准相结合的亚像元配准方法，具体配准过程如图3所示，基于相位相关的粗配准方法，其基本原理为：设参考模板和待配准模板分别为u₁(m,n)和u₂(m,n)，对应的傅里叶变换分别为：

其中ξ和η为空间频率变量。模板u₁(m,n)和u₂(m,n)的互功率谱为：

其中U₁ ^*(ξ,η)为U₁(ξ,η)的复共轭。对互功率谱C_ps(ξ,η)进行傅里叶反变换，得到二维脉冲函数δ(i,j)，如式(15)所示，二维脉冲函数的峰值点(i_p,j_p)即为像元级粗配准结果。

IFFT2{C_ps(ξ,η)}＝IFFT2{e^{j2π(iξ+jη)}}＝δ(i,j) (3)

基于互相关二维曲面拟合的精配准方法，其基本原理为：根据粗配准结果，在重叠图像O₁和O₂中选取模板u₁(m,n)和u₂(m-i_p,n-j_p)，利用归一化互相关系数作为衡量两个模板相似度的评价指标，表达式如下：

其中C(u,v)为两个模板的中心点相距(u,v)时的归一化互相关系数。利用最小二乘法求解互相关系数二次曲面的表达式，曲面的峰值为归一化互相关系数C(u,v)的最大值点，此时对应的两个模板的中心点即为同名配准点。

利用上述图像配准方法对重叠图像O₁和O₂进行逐行匹配处理，获得同名配准点集合P＝{(p₁ ¹,p₂ ¹),(p₁ ²,p₂ ²),…,(p₁ ^m,p₂ ^m),,…,(p₁ ^M,p₂ ^M)}，其中(p₁ ^m,p₂ ^m)为重叠图像O₁和O₂的第m行配准后得到的一对同名配准点，1≤m≤M，并有： 和分别为点p₁ ^m在所述重叠图像O₁内的行坐标和列坐标，和分别为点p₂ ^m在所述重叠图像O₂内的行坐标和列坐标，

步骤2.2、分别计算同名配准点集合P中每对同名配准点内两点的行坐标差值和列坐标差值，获得沿TDI方向上的相对成像位置差集合S^||＝{s^||(1),s^||(2)，…,s^||(m)，…,s^||(M)}和垂直于TDI方向上的相对成像位置差集合S^⊥＝{s^⊥(1),s^⊥(2)，…,s^⊥(m)，…,s^⊥(M)}，其中s^||(m)为所述重叠图像O₁和O₂的第m行图像配准后得到的沿TDI方向上的相对成像位置差，且s^⊥(m)为所述重叠图像O₁和O₂的第m行图像配准后得到的垂直于TDI方向上的相对成像位置差，且

步骤2.3、将沿TDI方向上的相对成像位置差集合S^||从第一个元素开始每连续的L个元素排成一行，从而形成N行×L列的二维矩阵S_∑ ^||：

其中

步骤2.4、将垂直于TDI方向上的相对成像位置差集合S^⊥从第一个元素开始每连续的L个元素排成一行，从而形成N行×L列的二维矩阵S_∑ ^⊥：

其中

步骤3、卫星平台颤振的估计：

步骤3.1、定义沿TDI方向上的颤振为

其中，g^||(v·L+u)为(v·L+u)×T_r时刻沿TDI方向的颤振，g_||(v·L+u)∈Q，u∈Z,v∈Z，1≤u≤L，0≤v≤N，Z表示整数集合；

定义垂直于TDI方向上的颤振为

其中，g^⊥(v·L+u)为(v·L+u)×T_r时刻垂直于TDI方向的颤振，g_⊥(v·L+u)∈Q，u∈Z,v∈Z，1≤u≤L，0≤v≤N；

定义颤振为G_∑＝{g_ij}_(N+1)×L，其中 为沿TDI方向的单位矢量，为垂直于TDI方向的单位矢量，1≤i≤N+1，1≤j≤L；

定义沿TDI方向上颤振G^∑||中的第一行元素[g^||(1),g^||(2),…,g^||(L)]为沿TDI方向上的颤振子块G₁ ^||；

定义垂直于TDI方向上颤振G_∑ ^⊥中的第一行元素[g^⊥(1),g^⊥(2),…,g^⊥(L)]为垂直于TDI方向上的颤振子块G₁ ^⊥；

定义沿TDI方向上颤振G_∑ ^||中的第二行元素至第N+1行元素

为沿TDI方向上的颤振子块G₂ ^||；

定义垂直于TDI方向上颤振G_∑ ^⊥中的第二行元素至第N+1行元素

为垂直于TDI方向上的颤振子块G₂ ^⊥；

步骤3.2、将拍摄任务中[T_r,L·T_r]时间段内的卫星在轨运行参数带入文献“WangJiaqi,Yu Ping,Yan Changxiang,Ren Jianyue,Hebin.Space optical remote sensorimage motion velocity vector computational modeling,error budget andsynthesis[J].Chinese Optics Letters,2005,(07):414-417.”给出的星下点对地成像模型中，分别得到沿TDI方向上颤振子块G₁ ^||在1Hz以下的低频分量和垂直于TDI方向上颤振子块G₁ ^⊥在1Hz以下的低频分量其中表示u×T_r时刻沿TDI方向上颤振在1Hz以下的低频分量，表示u×T_r时刻垂直于TDI方向上颤振在1Hz以下的低频分量，

步骤3.3、以沿TDI方向上颤振子块G₁ ^||与其低频分量之间的相对残差的平方和作为沿TDI方向上的目标函数H(G₁ ^||)，如式(5)所示：

式(5)中，g^||(i)表示i×T_r时刻沿TDI方向上的颤振，表示j×T_r时刻沿TDI方向上颤振在1Hz以下的低频分量；

以垂直于TDI方向上颤振子块G₁ ^⊥与其低频分量之间的相对残差的平方和作为垂直于TDI方向上的目标函数H(G₁ ^⊥)，如式(6)所示：

式(6)中，g^⊥(i)表示i×T_r时刻垂直于TDI方向上的颤振，表示j×T_r时刻垂直于TDI方向上颤振在1Hz以下的低频分量；

步骤3.4、利用优化算法中的收敛速率和运算量均居中的共轭梯度法分别求解目标函数H(G₁ ^||)和H(G₁ ^⊥)的最小值点和并分别作为沿TDI方向上颤振子块G₁ ^||和垂直于TDI方向上颤振子块G₁ ^⊥的估计值，如式(7)和式(8)所示：

步骤3.5、根据空间相机的推扫成像原理和TDICCD拼接区成像特点推导颤振探测模型，具体推导过程如下：

根据推扫成像原理，同一目标会被相邻两片TDICCD的拼接区分时成像两次，而卫星平台颤振是随时间不断变化的，拼接区OA₁和OA₂对同一目标两次分时成像所对应的卫星平台颤振状态可能不同，进而导致该目标在拼接区OA₁和OA₂中的成像位置存在差异，如图4所示，以垂直于TDI方向为例，相对成像位置差与颤振的关系可以表示为：

s^⊥(t)＝g^⊥(t+Δt)-g^⊥(t) (9)

式(9)中，g^⊥(t)为t时刻垂直于TDI方向上的颤振，Δt为拼接区OA₁和OA₂对同一目标两次成像的时间间隔，且Δt＝L×T_r，L为相邻两片TDICCD之间的间隔行数，T_r为每行图像的曝光时长，s^⊥(t)为t时刻垂直于TDI方向上的相对成像位置差。

对式(9)以每行图像的曝光时长T_r为采样周期进行离散化并整理后得到：

g^⊥(m+L)＝g^⊥(m)+s^⊥(m) (10)

式(10)中，g^⊥(m)为m×T_r时刻垂直于TDI方向上的颤振，g^⊥(m+L)为(m+L)×T_r时刻垂直于TDI方向上的颤振，s^⊥(m)为m×T_r时刻垂直于TDI方向上的相对成像位置差，1≤m≤M。

对垂直于TDI方向上的颤振g^⊥(m+L),m＝1,2,…,M从第一个元素开始每连续的L个元素分为一组，表示为：

式(11)中

根据垂直于TDI方向上的颤振子块G₁ ^⊥、垂直于TDI方向上的颤振子块G₂ ^⊥和垂直于TDI方向上的相对成像位置差S_∑ ^⊥的定义，式(11)可以表示为：

令系数矩阵系数矩阵则式(12)可以表示为：

G₂ ^⊥＝A·G₁ ^⊥+B·S_∑ ^⊥ (13)

式(13)即为垂直于TDI方向上的颤振探测模型，同理，推导沿TDI方向上的颤振探测模型如式(14)所示：

G₂ ^||＝A·G₁ ^||+B·S_∑ ^|| (14)

式(14)中G₁ ^||、G₂ ^||和S_∑ ^||的含义同说明书中的定义，系数矩阵A和B的取值同上；

步骤3.6、根据式(14)所示沿TDI方向的颤振探测模型，利用沿TDI方向上颤振子块G₁ ^||的估计值和沿TDI方向上的相对成像位置差矩阵S_∑ ^||计算沿TDI方向上颤振子块G₂ ^||的估计值如式(15)所示：

根据式(13)所示垂直于TDI方向的颤振探测模型，利用垂直于TDI方向上颤振子块G₁ ^⊥的估计值和垂直于TDI方向上的相对成像位置差矩阵S_∑ ^⊥计算垂直于TDI方向上颤振子块G₂ ^⊥的估计值如式(16)所示：

步骤3.7、根据沿TDI方向上颤振的定义，将沿TDI方向上颤振子块G₁ ^||的估计值和沿TDI方向上颤振子块G₂ ^||的估计值合成沿TDI方向上的颤振G_∑ ^||的估计值如式(17)所示：

根据垂直于TDI方向上颤振的定义，将垂直于TDI方向上颤振子块G₁ ^⊥的估计值和垂直于TDI方向上颤振子块G₂ ^⊥的估计值合成垂直于TDI方向上的颤振G_∑ ^⊥的估计值如式(18)所示：

下面以我国某高分辨率TDICCD相机为例对本发明提出方法进行说明。该相机的焦面排布方式如图5所示，五片TDICCD交错拼接排布，在沿TDI方向上，任意相邻两片TDICCD的间隔行数L＝3028，在垂直于TDI积分方向上，任意相邻两片TDICCD拼接区内的像元列数C＝40，一次拍摄任务持续时长T_w＝30s，每行图像的曝光时长T_r＝73μs。本实验选取该相机于北纬38.8436°～39.4515°升轨阶段TDICCD₃和TDICCD₄拍摄的遥感图像，利用本发明方法探测图像拍摄期间的卫星平台颤振，TDICCD₃和TDICCD₄分别拍摄411520行×4096列图像，如图6所示，左侧为TDICCD₃拍摄的部分遥感图像，右侧为TDICCD₄拍摄的部分遥感图像，虚线内为部分拼接区重叠图像。

步骤1、重叠图像的获取：

利用TDICCD₃和TDICCD₄的拼接区OA₁和OA₂获得图像I₁＝(I₁ ¹,I₁ ²,…,I₁ ⁴¹¹⁵²⁰)^T和I₂＝(I₂ ¹,I₂ ²,…,I₂ ⁴¹¹⁵²⁰)^T，提取图像I₁和I₂中的重叠部分(I₁ ¹,I₁ ²,…,I₁ ⁴⁰⁸⁴⁹²)^T和(I₂ ³⁰²⁹,I₂ ³⁰³⁰,…,I₂ ⁴¹¹⁵²⁰)^T，分别记为重叠图像O₁和O₂；

步骤2、相对成像位置差的计算：

利用图像配准方法对重叠图像O₁和O₂进行逐行匹配处理。逐行匹配处理具体为：在重叠图像O₁中，每行图像共有40个像素点，分别选取其中的第10、20、30个像素点作为配准点，以每个配准点为中心的8行*8列邻域为参考模板，选取重叠图像O₂内对应位置处的模块作为待配准模板，采用粗配准和精配准相结合的亚像元配准方法对两幅模板进行匹配处理，获得一对同名配准点，对同行图像内三个配准点的配准结果取算术平均值作为整行图像的配准结果，对重叠图像O₁的第11行至第408482行图像按照上述步骤进行逐行匹配处理，获得同名配准点集合。分别计算同名配准点集合中每对同名配准点内两点的行坐标差值和列坐标差值，获得沿TDI方向上的相对成像位置差集合和垂直于TDI方向上的相对成像位置差集合；

步骤3、卫星平台颤振的估计：

将拍摄任务中[73μs,3028×73μs]时间段内卫星在轨运行参数带入星下点对地成像模型中，分别得到沿TDI方向上颤振子块G₁ ^||在1Hz以下的低频分量和垂直于TDI方向上颤振子块G₁ ^⊥在1Hz以下的低频分量以沿TDI方向上颤振子块G₁ ^||与其低频分量之间的相对残差的平方和作为沿TDI方向上的目标函数H(G₁ ^||)，以垂直于TDI方向上颤振子块G₁ ^⊥与其低频分量之间的相对残差的平方和作为垂直于TDI方向上的目标函数H(G₁ ^⊥)，利用共轭梯度法分别求解目标函数H(G₁ ^||)和H(G₁ ^⊥)的最小值点和并分别作为沿TDI方向上颤振子块G₁ ^||和垂直于TDI方向上颤振子块G₁ ^⊥的估计值，根据沿TDI方向的颤振探测模型，利用沿TDI方向上颤振子块G₁ ^||的估计值和沿TDI方向上的相对成像位置差矩阵S_∑ ^||计算沿TDI方向上颤振子块G₂ ^||的估计值根据垂直于TDI方向的颤振探测模型，利用垂直于TDI方向上颤振子块G₁ ^⊥的估计值和垂直于TDI方向上的相对成像位置差矩阵S_∑ ^⊥计算垂直于TDI方向上颤振子块G₂ ^⊥的估计值将沿TDI方向上颤振子块G₁ ^||的估计值和沿TDI方向上颤振子块G₂ ^||的估计值合成沿TDI方向上的颤振G_∑ ^||的估计值将垂直于TDI方向上颤振子块G₁ ^⊥的估计值和垂直于TDI方向上颤振子块G₂ ^⊥的估计值合成垂直于TDI方向上的颤振G_∑ ^⊥的估计值

颤振探测结果如图7a、图7b、图7c、图7d所示，垂直于TDI方向上，在0.24Hz、4.54Hz、9.08Hz、13.53Hz、18Hz和22.48Hz处存在振动峰值点；沿TDI方向上，除了上述峰值点，在0.13Hz处存在振动峰值。表明本发明方法能够对卫星平台颤振进行有效的探测。

具体实现可参见上述方法的相应说明。

Claims (1)

1.一种基于TDICCD拼接区图像的卫星平台颤振探测方法，其特征按如下步骤进行：

步骤1、重叠图像的获取：

步骤1.1、在卫星平台上搭载有高分辨率空间相机，所述高分辨率空间相机的焦平面上存在多片交错拼接排布的TDICCD，任意选取两片相邻的TDICCD，记为TDICCD₁和TDICCD₂，在沿TDI方向上，TDICCD₁和TDICCD₂的间隔行数记为L，在垂直于TDI方向上，TDICCD₁和TDICCD₂的交接处存在若干列重叠像元形成拼接区，记为OA₁和OA₂，所述拼接区OA₁和OA₂内的重叠像元列数记为C；

所述高分辨率空间相机任意一次拍摄任务的起始时刻记为0时刻，拍摄持续时长记为T_w，其中T_w∈Q，T_w＞0，Q表示有理数集合；

利用所述拼接区OA₁和OA₂在[0,T_w]期间对探测区域内的所有目标进行两次分时拍摄，获得图像I₁＝(I₁ ¹,I₁ ²,…,I₁ ^r,…,I₁ ^R)^T和I₂＝(I₂ ¹,I₂ ²,…,I₂ ^r,…,I₂ ^R)^T，其中为所述图像I₁和I₂的总行数，T_r为每行图像的曝光时长，I₁ ^r和I₂ ^r分别为所述图像I₁和I₂的第r行图像，T_r∈Q，T_r＞0，1≤r≤R；

步骤1.2、提取所述图像I₁和I₂中的重叠部分(I₁ ¹,I₁ ²,…,I₁ ^R-L)^T和(I₂ ^L+1,I₂ ^L+2,…,I₂ ^R)^T，分别记为重叠图像O₁和O₂，所述重叠图像O₁和O₂的总列数为C，总行数为记为M＝R-L；

步骤2、相对成像位置差的计算：

步骤2.1、利用图像配准方法对所述重叠图像O₁和O₂进行逐行匹配处理，获得同名配准点集合P＝{(p₁ ¹,p₂ ¹),(p₁ ²,p₂ ²),…,(p₁ ^m,p₂ ^m),…,(p₁ ^M,p₂ ^M)}，其中(p₁ ^m,p₂ ^m)为所述重叠图像O₁和O₂的第m行图像配准后得到的一对同名配准点，1≤m≤M，并有： 和分别为点p₁ ^m在所述重叠图像O₁内的行坐标和列坐标，和分别为点p₂ ^m在所述重叠图像O₂内的行坐标和列坐标，

步骤2.2、分别计算所述同名配准点集合P中每对同名配准点内两点的行坐标差值和列坐标差值，获得沿TDI方向上的相对成像位置差集合S^||＝{s^||(1),s^||(2)，…,s^||(m)，…,s^||(M)}和垂直于TDI方向上的相对成像位置差集合S^⊥＝{s^⊥(1),s^⊥(2)，…,s^⊥(m)，…,s^⊥(M)}，其中s^||(m)为所述重叠图像O₁和O₂的第m行图像配准后得到的沿TDI方向上的相对成像位置差，且s^⊥(m)为所述重叠图像O₁和O₂的第m行图像配准后得到的垂直于TDI方向上的相对成像位置差，且

步骤2.3、将所述沿TDI方向上的相对成像位置差集合S^||从第一个元素开始每连续的L个元素排成一行，从而形成N行×L列的二维矩阵S_∑ ^||：

其中

步骤2.4、将所述垂直于TDI方向上的相对成像位置差集合S^⊥从第一个元素开始每连续的L个元素排成一行，从而形成N行×L列的二维矩阵S_∑ ^⊥：

其中

步骤3、卫星平台颤振的估计：

步骤3.1、定义沿TDI方向上的颤振为：

其中，g^||(v·L+u)为(v·L+u)×T_r时刻沿TDI方向的颤振，g_||(v·L+u)∈Q，u∈Z，v∈Z，1≤u≤L，0≤v≤N，Z表示整数集合；

定义垂直于TDI方向上的颤振为：

其中，g^⊥(v·L+u)为(v·L+u)×T_r时刻垂直于TDI方向的颤振，g_⊥(v·L+u)∈Q，u∈Z,v∈Z，1≤u≤L，0≤v≤N；

定义颤振为G_∑＝{g_ij}_(N+1)×L，其中 为沿TDI方向的单位矢量，为垂直于TDI方向的单位矢量，1≤i≤N+1，1≤j≤L；

定义沿TDI方向上颤振G_∑ ^||中的第一行元素[g^||(1),g^||(2),…,g^||(L)]为沿TDI方向上的颤振子块G₁ ^||；

定义垂直于TDI方向上颤振G_∑ ^⊥中的第一行元素[g^⊥(1),g^⊥(2),…,g^⊥(L)]为垂直于TDI方向上的颤振子块G₁ ^⊥；

定义沿TDI方向上颤振G_∑ ^||中的第二行元素至第N+1行元素为沿TDI方向上的颤振子块G₂ ^||；

定义垂直于TDI方向上颤振G_∑ ^⊥中的第二行元素至第N+1行元素为垂直于TDI方向上的颤振子块G₂ ^⊥；

步骤3.2、将所述拍摄任务中[T_r,L·T_r]时间段内的卫星在轨运行参数带入星下点对地成像模型中，分别得到沿TDI方向上颤振子块G₁ ^||在1Hz以下的低频分量和垂直于TDI方向上颤振子块G₁ ^⊥在1Hz以下的低频分量其中表示u×T_r时刻沿TDI方向上颤振在1Hz以下的低频分量，表示u×T_r时刻垂直于TDI方向上颤振在1Hz以下的低频分量，

步骤3.3、以所述沿TDI方向上颤振子块G₁ ^||与其低频分量之间的相对残差的平方和作为沿TDI方向上的目标函数H(G₁ ^||)，如式(1)所示：

式(1)中，g^||(i)表示i×T_r时刻沿TDI方向上的颤振，表示j×T_r时刻沿TDI方向上颤振在1Hz以下的低频分量；

以所述垂直于TDI方向上颤振子块G₁ ^⊥与其低频分量之间的相对残差的平方和作为垂直于TDI方向上目标函数H(G₁ ^⊥)，如式(2)所示：

式(2)中，g^⊥(i)表示i×T_r时刻垂直于TDI方向上的颤振，表示j×T_r时刻垂直于TDI方向上颤振在1Hz以下的低频分量；

步骤3.4、利用优化算法分别求解所述目标函数H(G₁ ^||)和H(G₁ ^⊥)的最小值点和并分别作为所述沿TDI方向上颤振子块G₁ ^||和垂直于TDI方向上颤振子块G₁ ^⊥的估计值，如式(3)和式(4)所示：

步骤3.5、根据空间相机的推扫成像原理和TDICCD拼接区成像特点，建立如式(5)和式(6)所示的沿TDI方向和垂直于TDI方向的颤振探测模型：

G₂ ^||＝A·G₁ ^||+B·S_∑ ^|| (5)

G₂ ^⊥＝A·G₁ ^⊥+B·S_∑ ^⊥ (6)

式(5)和式(6)中，A和B均表示系数矩阵，且

步骤3.6、根据式(5)所示沿TDI方向的颤振探测模型，利用式(7)计算沿TDI方向上颤振子块G₂ ^||的估计值

根据式(6)所示垂直于TDI方向的颤振探测模型，利用式(8)计算垂直于TDI方向上颤振子块G₂ ^⊥的估计值

步骤3.7、根据沿TDI方向上颤振的定义，利用式(9)合成沿TDI方向上的颤振G_∑ ^||的估计值

根据垂直于TDI方向上颤振的定义，利用式(10)合成垂直于TDI方向上的颤振G_∑ ^⊥的估计值