CN103983343B - 一种基于多光谱影像的卫星平台震颤检测方法及系统 - Google Patents

Abstract

本发明提供了一种基于多光谱影像的卫星平台震颤检测方法及系统，包括：数据准备，包括对原始的多光谱影像进行预处理，得到未配准多光谱影像；波段选择，包括选取用于检测的两个波段；密集匹配，包括选取样本点，然后利用影像相关匹配得到像素级匹配结果，最后利用最小二乘匹配得到子像素精度的同名像点坐标，得到匹配同名点；震颤分析，包括首先对匹配同名点逐点计算配准误差，然后进行逐行分析得到配准误差随成像行的变化，最后通过傅里叶变换对平台震颤时序数据进行频谱分析，获取震颤引起的频谱特征。本发明在无高精度高频率震颤测量设备实现了卫星平台震颤检测，可准确获取震颤频率。

Description

一种基于多光谱影像的卫星平台震颤检测方法及系统

技术领域

本发明属于卫星遥感地面预处理领域，特别是涉及到一种利用多光谱影像检测卫星平台震颤的技术方案。

背景技术

卫星平台震颤，是指卫星在轨运行期间，星上运动部件周期性运动或因变轨冷热交变等因素引发的扰动，使星体产生一种幅值较小、频率较高的颤振响应。平台震颤具有微小性、固有性、难控性和敏感性等特点。随着遥感卫星空间分辨率的提高，卫星平台震颤对成像质量的影响越来越受到卫星总体设计、相机设计以及地面处理与应用部门的重视。日本NASDA利用试验卫星ETS-VI进行了卫星振动测量试验，卫星平台角振动在0.39Hz到250Hz范围内的径向均方根为16.3μrad(约3.36角秒)，其中83.6％的振动能量集中在0.39～10Hz之间，99％以上的能量在102Hz以下。哈尔滨工业大学通过理论分析和仿真实验研究了平台运动对星载TDICCD成像质量的影响，研究结果表明平台振动会引起不同程度的图像模糊和几何变形，将直接影响卫星影像产品的质量。

对于遥感卫星平台震颤，一方面需要通过动力学方法采取隔振和抑振措施，避免和减小其对成像的影像，另外一方面需通过直接或者间接的方式进行平台震颤测量，对其引起的像质下降和几何变形进行补偿。然而平台震颤不可避免，难以控制，很难做到百分之百的抑振和隔振，因此对平台震颤的精确测量才能保证高分辨率遥感影像的高质量、高精度处理。卫星平台震颤直接测量是利用测量频率和测量精度均较高的角度传感器进行测量，获取平台角位移、角速度和角加速度等信息。常用的角度传感器有MHD角速度传感器和FDI角位移传感器，测量精度优于0.2角秒，测量频率高达上千赫兹，在国外在轨卫星中已有应用，但目前国内暂无在轨卫星搭载这类姿态测量载荷，因此还无法获取在轨卫星平台的震颤情况。平台震颤间接测量是利用非姿态测量装置对平台震颤进行检测与估计。基于遥感影像的卫星平台震颤检测是进行平台震颤分析的思路之一。法国国家空间研究中心(CNES)提出基于立体像对的卫星平台微振动探测方法，分别对翻滚、俯仰和偏航三个方向的平台震颤进行识别和移除，但仅做了仿真实验，并未对真实数据进行处理。日本东京大学也发表了利用多光谱影像进行卫星震颤检测和补偿的研究论文，但文中提出的方法和实验仅对垂直于推扫方向(即飞行方向)的配准误差进行了检测，且缺乏必要频谱分析。不论在平台震颤直接测量还是间接检测方面，国内还暂无通过遥感影像对卫星平台震颤进行检测分析方面的公开发表论文与专利。

发明内容

本发明针对卫星平台震颤难以测量的问题，提供了一种基于多光谱影像的卫星平台震颤检测技术方案。

本发明提供的技术方案提供一种基于多光谱影像的卫星平台震颤检测方法，包括以下步骤：

步骤1，数据准备，包括对原始的多光谱影像进行预处理，得到未配准多光谱影像；

步骤2，波段选择，包括选取用于检测的两个波段；

步骤3，密集匹配，包括选取样本点，然后利用影像相关匹配得到像素级匹配结果，最后利用最小二乘匹配得到子像素精度的同名像点坐标，得到匹配同名点；

步骤4，震颤分析，包括首先对匹配同名点逐点计算配准误差，然后进行逐行分析得到配准误差随成像行的变化，最后通过傅里叶变换对平台震颤时序数据进行频谱分析，获取震颤引起的频谱特征。

而且，步骤2中选取的两个波段为相邻的可见光波段。

而且，步骤3中选取样本点的方式为逐点或者以一定间隔选取。

而且，步骤3中利用影像相关匹配得到像素级匹配结果，包括利用初始CCD设计值，计算同名像点的初始位置；在同名像点初始位置的八邻域范围内搜索同名像点的位置，以相关系数作为影像匹配的匹配测度，在搜索匹配目标中选取相关系数最大的点作为同名点。

而且，步骤4中对匹配同名点计算配准误差，实现方式如下，

设步骤2选择的两个波段为M、N，设第u扫描行中匹配出V个点，第u扫描行中第k个同名像点的配准误差计算如下，

${\mathrm{\Δx}}_{u,k}=x_{M_{u,k}}-x_{N_{u,k}}-\mathrm{\Δx}$

${\mathrm{\Δy}}_{u,k}=y_{M_{u,k}}-y_{N_{u,k}}-\mathrm{\Δy}$

其中，Δx_u,k，Δy_u,k分别为沿CCD方向和沿推扫方向的配准误差，分别为第u行的第k个同名像点在M，N两个波段影像上的坐标，Δx,Δy为两个波段CCD之间在沿CCD和沿推扫方向的平移设计值。

而且，步骤4中进行逐行分析得到配准误差随成像行的变化，实现方式如下，

设第u扫描行最终匹配出V个点，统计沿CCD方向的配准误差平均值Δx(u)和沿推扫方向的配准误差平均值Δy(u)，

$\mathrm{\Δx}\left(u\right)=\underset{k=1}{\overset{k=V}{\mathrm{\Σ}}}{\mathrm{\Δx}}_{u,k}/V$

$\mathrm{\Δy}\left(u\right)=\underset{k=1}{\overset{k=V}{\mathrm{\Σ}}}{\mathrm{\Δy}}_{u,k}/V$

逐行统计得到沿CCD方向和沿推扫方向的配准误差平均值Δx(u)、Δy(u)随成像行相应成像时间的变化规律曲线。

而且，步骤4中通过傅里叶变换对平台震颤时序数据进行频谱分析，实现方式如下，

设多光谱影像每一行的积分时间为Δt，那么对于U个成像行影像，成像总时间为T＝Δt·U，对于第u成像行对应成像时刻为t_u＝t₀+Δt·u，u＝0,1,……U-1，其中t₀为起始成像时间，

对Δx(u)的离散序列，进行傅里叶变换如下，

$a_K=\frac{2}{U}\underset{u=0}{\overset{U-1}{\mathrm{\Σ}}}\mathrm{\Δx}\left(u\right)\cos \frac{2\mathrm{\πKu}}{U}$

$b_K=\frac{2}{U}\underset{u=0}{\overset{U-1}{\mathrm{\Σ}}}\mathrm{\Δx}\left(u\right)\sin \frac{2\mathrm{\πKu}}{U}$

$c_K=\sqrt{a_K^2+b_K^2}$

其中，K为分解后序号，取值范围从0～U/2，a_K和b_K分别为K次谐波余弦和正弦函数的振幅，c_K表示K次谐波的振幅大小；K次谐波对应频率f_K为，

$f_K=\frac{K}{U\·\mathrm{\Δt}}=\frac{K}{T}$

得到平台震颤在沿CCD方向的频谱曲线；

对Δy(u)的离散序列，进行傅里叶变换如下，

$a_K=\frac{2}{U}\underset{u=0}{\overset{U-1}{\mathrm{\Σ}}}\mathrm{\Δy}\left(u\right)\cos \frac{2\mathrm{\πKu}}{U}$

$b_K=\frac{2}{U}\underset{u=0}{\overset{U-1}{\mathrm{\Σ}}}\mathrm{\Δy}\left(u\right)\sin \frac{2\mathrm{\πKu}}{U}$

$c_K=\sqrt{a_K^2+a_K^2}$

其中，K为分解后序号，取值范围从0～U/2，a_K和b_K分别为K次谐波余弦和正弦函数的振幅，c_K表示K次谐波的振幅大小；K次谐波对应频率f_K为，

$f_K=\frac{K}{U\·\mathrm{\Δt}}=\frac{K}{T}$

得到平台震颤在沿推扫方向的频谱曲线。

本发明还相应提供了一种基于多光谱影像的卫星平台震颤检测系统，包括以下模块：

数据准备模块，用于对原始的多光谱影像进行预处理，得到未配准多光谱影像；

波段选择模块，用于选取用于检测的两个波段；

密集匹配模块，用于选取样本点，然后利用影像相关匹配得到像素级匹配结果，最后利用最小二乘匹配得到子像素精度的同名像点坐标，得到匹配同名点；

震颤分析模块，用于首先对匹配同名点逐点计算配准误差，然后进行逐行分析得到配准误差随成像行的变化，最后通过傅里叶变换对平台震颤时序数据进行频谱分析，获取震颤引起的频谱特征。

本发明提供了一种基于多光谱影像的卫星平台震颤检测技术方案，实现了无需高精度、高频率震颤测量设备，即可对卫星平台沿推扫方向和垂直于推扫方向进行震颤检测，通过频谱分析，准确获取震颤频率，为进一步提高影像处理质量和几何精度提供了依据，也为寻找平台震颤源和优化卫星平台设计提供了线索。该方法具有简单、有效、精度高、易于实现的特点。

附图说明

图1为本发明实施例的流程图。

图2为本发明实施例的多光谱相机焦平面设计示意图。

图3为本发明实施例的多光谱相机焦平面相应角度示意图。

图4为本发明实施例的震颤分析时序图。

图5为本发明实施例的震颤分析频谱图。

具体实施方式

本发明技术方案可采用计算机软件方式支持自动运行流程。以下结合附图和实施例详细说明本发明技术方案。

图1表示的是基于多光谱影像的卫星平台震颤检测技术流程图，以下针对实施例流程中的各步骤，对本发明方法做进一步详细描述。

(1)数据准备：对原始的多光谱影像进行辐射校正、CCD拼接等预处理，得到未配准多光谱影像。

为提高匹配精度，需对原始的多光谱影像进行相对辐射校正处理，消除条带噪声对后续处理的影响。另外对于共线TDICCD光学拼接影像采用去重叠处理即可实现CCD拼接，得到完整的全视场影像；对于非共线TDICCD视场拼接影像无需拼接，以避免拼接过程对影像几何特性的改变。另外，若考虑CCD变形和镜头畸变等问题，也可不进行CCD拼接处理，以减小相机内部畸变对检测结果的影响。为了更好的保留震颤对多光谱影像的配准误差影响，无需对多光谱影像进行配准处理。

(2)波段选择：选取一组最优组合的两个波段，减小匹配误差对检测结果的影响。

波段选择包括两项原则：1)尽可能选择光谱特性相似的波段，避免光谱特性差异降低匹配精度；2)尽可能选择相邻波段，由于震颤幅值通常较小、频率较高，采用相邻波段进行检测，结果更精细。以光学卫星搭载的多光谱相机为例，多光谱影像通常由蓝、绿、红、红外四个谱段组成，间隔相同距离依次排列在焦平面上。蓝、绿、红、红外谱段分别标示为B1、B2、B3、B4，见图2，图中o为像主点(此处主点在B1波段CCD上)，y轴表示推扫方向，x轴表示沿CCD方向，形成焦平面坐标系o-xy。其他三个波段则在成像视场中形成偏视场成像，见图3，O为成像投影中心，θ₂、θ₃、θ₄分别为B2、B3、B4的偏视场角。由于红外波段与可见光三个波段的成像特性存在较大差异，故采用可见光波段进行震颤检测，降低由于光谱特性差异过大对匹配精度的影响；为了更精细地检测震颤，采用相邻的可见光波段进行检测。因此可采用蓝、绿组合或者绿、红组合进行震颤检测。

(3)高精度密集匹配：为保证样本点分布均匀、数量充足，逐点或者以一定间隔选取样本点，利用影像相关匹配和最小二乘匹配逐步得到子像素精度的同名像点坐标。

设选择的两个波段为M、N，实施例根据(1)所得未配准多光谱影像中M、N波段的影像进行匹配。首先逐点或者以一定间隔在M波段的影像上选取样本点，保证样本点分布均匀、数量充足，然后利用影像相关匹配得到像素级匹配结果，最后利用最小二乘匹配得到子像素精度的同名像点坐标。间隔可由本领域技术人员自行预先设置。

1)影像相关匹配

针对M波段的影像上每一个样本点，利用初始CCD设计值，计算样本点在N波段的影像上同名像点的初始位置。考虑到波段间CCD相对变形较小，可以在N波段的影像上同名像点初始位置的八邻域范围内搜索同名像点的位置。以相关系数作为影像匹配的匹配测度，在搜索匹配目标中选取相关系数最大的点作为同名点。相关系数是影像匹配中最常使用的匹配测度，它是标准化的协方差函数。在实际的应用中，采用如下公式进行计算：

$(c,r)=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}(g_{i,j}\·g_{i+r,j+c}^{\′})-\frac{1}{m\·n}\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}\right)\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i+r,j+c}^{\′}\right)}{\underset{}{\overset{}{\sqrt{[\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g_{i,j}}^2-\frac{1}{m\·n}{\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}\right)}^2][\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g_{i+r,j+c}^{\′}}^2-\frac{1}{m\·n}{\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i+r,j+c}^{\′}\right)}^2]}}}}---\left(1\right)$

其中，m,n表示匹配窗口的行列大小，i,j表示行数和列数，g_i,j、g_i'_+r,j+c分别表示左右影像匹配窗口内的灰度值，c,r相当于是搜索区影像(右影像，即N波段的影像)对于参考影像(左影像，即M波段的影像)的在行、列上的位移。在此，取c＝r＝0。相关系数是灰度线性变换不变量，可以克服灰度的线性变形，保持相关系数的不变性。

2)最小二乘匹配

为了得到子像素的匹配结果，本文采用最小二乘匹配算法得到高精度的检测结果。最小二乘影像匹配(LeastSquaresImageMatching,LSM)是德国Ackerman教授提出的高精度匹配算法，该方法充分利用了影像窗口内的信息进行平差，可得到子像素级的匹配结果。采用单点最小二乘匹配算法，在考虑灰度线性变形，几何仿射变形的条件下对左右匹配窗口内的每一个像素灰度值建立等式(公式2)，采用最小二乘方法解算灰度线性变形参数和几何仿射变形参数，再将参考影像中的匹配窗口灰度中心化获得子像素特征点，利用几何仿射变形参数求解同名点在匹配窗口中的坐标。

g₁(x₁,y₁)+n₁(x₁,y₁)＝h₀+h₁·g₂(x₂,y₂)+n₂(x₂,y₂)(2)

其中，x₁,y₁，x₂,y₂分别为左右匹配窗口对应的像素坐标，g₁(x₁,y₁)、g₂(x₂,y₂)分别为其灰度，n₁(x₁,y₁)、n₂(x₂,y₂)分别为两个灰度值的噪声，h₀、h₁分别为灰度线性变形系数。考虑几何仿射变形，则

x₂＝a₀+a₁x₁+a₂y₁

(3)

y₂=b₀+b₁x₁+b₂y₁其中，a₀、a₁、a₂、b₀、b₁、b₂分别为几何仿射变形系数。

最小二乘匹配的核心就是通过窗口内尽可能每一个点满足公式(2)的条件下，最优求解出仿射变形参数和灰度线性参数，即可得到匹配窗口内任一点在目标窗口内的子像素坐标。

(4)震颤分析：首先对匹配同名点逐点计算配准误差，然后进行逐行分析，得到配准误差随成像行(即成像时间)的变化，最后通过傅里叶变换对平台震颤时序数据进行频谱分析，获取震颤引起的配准误差周期(频率)、振幅等频谱特征。

通过密集匹配得到大量波段间同名像点，逐点计算配准误差(公式(4))。设第u扫描行中匹配出V个点，该行第k个同名像点的配准误差计算如下，

$\begin{array}{c}{\mathrm{\Δx}}_{u,k}=x_{M_{u,k}}-x_{N_{u,k}}-\mathrm{\Δx}\\ {\mathrm{\Δy}}_{u,k}=y_{M_{u,k}}-y_{N_{u,k}}-\mathrm{\Δy}\end{array}---\left(4\right)$

其中，Δx_u,k，Δy_u,k分别为沿CCD方向和沿推扫方向的配准误差，分别为第u行的第k个同名像点在M，N两个波段影像上的坐标，Δx,Δy为两个波段CCD之间在沿CCD和沿推扫方向的平移设计值。

考虑到多光谱相机视场角较小，假定在同一成像时间内每一行受到的平台震颤影响具有一致性，假设第u扫描行最终匹配出V个点，统计配准误差平均值，作为该行配准误差的最优估值(见公式(5))，通过逐行统计得到沿CCD方向(Δx(u))和沿推扫方向(Δy(u))的配准误差随成像行(即成像时间)的变化规律曲线。

$\mathrm{\Δx}\left(u\right)=\underset{k=1}{\overset{k=V}{\mathrm{\Σ}}}{\mathrm{\Δx}}_{u,k}/v$ (5)

$\mathrm{\Δy}\left(u\right)=\underset{k=1}{\overset{k=V}{\mathrm{\Σ}}}{\mathrm{\Δy}}_{u,k}/V$

在此基础上，对平台震颤进行频谱分析，主要是分析由平台震颤引起的配准误差的周期(频率)、振幅等频谱特征，从而可得到平台震颤的频率。此处采用离散Fourier变换对配准误差时序曲线进行频谱分析。若多光谱影像每一行的积分时间为Δt，那么对于U个成像行影像，成像总时间为T＝Δt·U，对于第u成像行对应成像时刻为t_u＝t₀+Δt·u(u＝0,1,……U-1)，其中t₀为起始成像时间。那么对Δx(u)(u＝0,1,……U-1)的离散序列，进行Fourier变换：

$a_K=\frac{2}{U}\underset{u=0}{\overset{U-1}{\mathrm{\Σ}}}\mathrm{\Δx}\left(u\right)\cos \frac{2\mathrm{\πKu}}{U}$ (6)

$b_K=\frac{2}{U}\underset{u=0}{\overset{U-1}{\mathrm{\Σ}}}\mathrm{\Δx}\left(u\right)\sin \frac{2\mathrm{\πKu}}{U}$

$c_K=\sqrt{a_K^2+b_K^2}---\left(7\right)$

其中，K为分解后序号，取值范围从0～U/2，a_K和b_K分别为K次谐波余弦和正弦函数的振幅，c_K表示K次谐波的振幅大小。K次谐波对应频率f_K为：

$f_K=\frac{K}{U\·\mathrm{\Δt}}=\frac{K}{T}---\left(8\right)$

由此得到配准误差在沿CCD方向(垂直于推扫方向)的频谱曲线。将公式(6)中的Δx(u)改成Δy(u)，并按照公式(7)、(8)计算，即可得到配准误差在推扫方向的频谱曲线。由于平台震颤表现为周期性运动，其引起的配准误差也是周期性变化，那么在配准误差频谱曲线峰值处的非零频率即为卫星平台震颤频率。

下面以资源三号一景多光谱数据为例对本发明提出方法进行说明：首先对原始数据进行辐射校正、CCD拼接等预处理，得到未配准的多光谱影像；然后选择资源三号多光谱影像的第二波段与第三波段(分别为绿波段B2和红波段B3)进行逐点密集匹配，根据两个波段的相对位置设计值Δx＝0,Δy＝128得到同名点的初始位置，采用相关系数在同名点初始位置的八邻域范围内寻找相关系数最大且优于0.9的点作为同名点，并继续进行最小二乘匹配得到子像素匹配结果；最后根据公式(4)计算同名点配准误差，利用公式(5)计算每一行的配准误差，得到误差随成像行的变化曲线(如图4所示，横坐标为行号，纵坐标为B2-B3波段配准误差)，可以明显看到配准误差随时间有周期性变化，表明卫星在成像阶段存在震颤。接着利用公式(6)～(8)进行频谱分析，所选取数据的成像积分时间Δt＝803ms。对沿CCD方向的配准误差Δx(u)进行频谱分析，Fourier变化结果显示：在0赫兹和0.66赫兹处出现了峰值，两处幅值分别为0.75和0.46像素。频率为0赫兹处的幅值表示配准误差的平移量；0.66赫兹处的幅值则表示由卫星平台震颤引起的周期性配准误差(如图5所示，横坐标为频率、单位为Hz，纵坐标为振幅、单位为像素)。由此说明资源三号在该景影像成像阶段存在约0.66赫兹的周期性平台震颤运动，即约1900成像行，这与时序图结果保持一致。

具体实施时，还可采用系统方式实现，实施例提供的一种基于多光谱影像的卫星平台震颤检测系统，包括以下模块：

数据准备模块，用于对原始的多光谱影像进行预处理，得到未配准多光谱影像；

波段选择模块，用于选取用于检测的两个波段；

密集匹配模块，用于选取样本点，然后利用影像相关匹配得到像素级匹配结果，最后利用最小二乘匹配得到子像素精度的同名像点坐标，得到匹配同名点；

震颤分析模块，用于首先对匹配同名点逐点计算配准误差，然后进行逐行分析得到配准误差随成像行的变化，最后通过傅里叶变换对平台震颤时序数据进行频谱分析，获取震颤引起的频谱特征。

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

本文中所描述的具体实例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (6)

1.一种基于多光谱影像的卫星平台震颤检测方法，包括以下步骤：

步骤1，数据准备，包括对原始的多光谱影像进行预处理，得到未配准多光谱影像；

步骤2，波段选择，包括选取用于检测的两个波段；

步骤3，密集匹配，包括选取样本点，然后利用影像相关匹配得到像素级匹配结果，最后利用最小二乘匹配得到子像素精度的同名像点坐标，得到匹配同名点；

步骤4，震颤分析，包括首先对匹配同名点逐点计算配准误差，然后进行逐行分析得到配准误差随成像行的变化，最后通过傅里叶变换对平台震颤时序数据进行频谱分析，获取震颤引起的频谱特征；

对匹配同名点计算配准误差，实现方式如下，

设步骤2选择的两个波段为M、N，设第u扫描行中匹配出V个点，第u扫描行中第k个同名像点的配准误差计算如下，

${\mathrm{\Δx}}_{u,k}=x_{M_{u,k}}-x_{N_{u,k}}-\mathrm{\Δ}x$

${\mathrm{\Δy}}_{u,k}=y_{M_{u,k}}-y_{N_{u,k}}-\mathrm{\Δ}y$

其中，Δx_u,k，Δy_u,k分别为沿CCD方向和沿推扫方向的配准误差，分别为第u行的第k个同名像点在M，N两个波段影像上的坐标，Δx,Δy为两个波段CCD之间在沿CCD和沿推扫方向的平移设计值；

进行逐行分析得到配准误差随成像行的变化，实现方式如下，

设第u扫描行最终匹配出V个点，统计沿CCD方向的配准误差平均值Δx(u)和沿推扫方向的配准误差平均值Δy(u)，

$\mathrm{\Δ}x\left(u\right)=\underset{k=1}{\overset{k=V}{\Σ}}{\mathrm{\Δx}}_{u,k}/V$

$\mathrm{\Δ}y\left(u\right)=\underset{k=1}{\overset{k=V}{\Σ}}{\mathrm{\Δy}}_{u,k}/V$

逐行统计得到沿CCD方向和沿推扫方向的配准误差平均值Δx(u)、Δy(u)随成像行相应成像时间的变化规律曲线。

2.根据权利要求1所述基于多光谱影像的卫星平台震颤检测方法，其特征在于：步骤2中选取的两个波段为相邻的可见光波段。

3.根据权利要求1或2所述基于多光谱影像的卫星平台震颤检测方法，其特征在于：步骤3中选取样本点的方式为逐点或者以一定间隔选取。

4.根据权利要求1或2所述基于多光谱影像的卫星平台震颤检测方法，其特征在于：步骤3中利用影像相关匹配得到像素级匹配结果，包括利用初始CCD设计值，计算同名像点的初始位置；在同名像点初始位置的八邻域范围内搜索同名像点的位置，以相关系数作为影像匹配的匹配测度，在搜索匹配目标中选取相关系数最大的点作为同名点。

5.根据权利要求1或2所述基于多光谱影像的卫星平台震颤检测方法，其特征在于：步骤4中通过傅里叶变换对平台震颤时序数据进行频谱分析，实现方式如下，

设多光谱影像每一行的积分时间为Δt，那么对于U个成像行影像，成像总时间为T＝Δt·U，对于第u成像行对应成像时刻为t_u＝t₀+Δt·u，u＝0,1,……U-1，其中t₀为起始成像时间，

对Δx(u)的离散序列，进行傅里叶变换如下，

$a_K=\frac{2}{U}\underset{u=0}{\overset{U-1}{\Σ}}\mathrm{\Δ}x\left(u\right)cos\frac{2\mathrm{\π}Ku}{U}$

$b_K=\frac{2}{U}\underset{u=0}{\overset{U-1}{\Σ}}\mathrm{\Δ}x\left(u\right)sin\frac{2\mathrm{\π}Ku}{U}$

$c_K=\sqrt{a_K^2+b_K^2}$

其中，K为分解后序号，取值范围从0～U/2，a_K和b_K分别为K次谐波余弦和正弦函数的振幅，c_K表示K次谐波的振幅大小；K次谐波对应频率f_K为，

$f_K=\frac{K}{U\·\mathrm{\Δ}t}=\frac{K}{T}$

得到平台震颤在沿CCD方向的频谱曲线；

对Δy(u)的离散序列，进行傅里叶变换如下，

$a_K=\frac{2}{U}\underset{u=0}{\overset{U-1}{\Σ}}\mathrm{\Δ}y\left(u\right)cos\frac{2\mathrm{\π}Ku}{U}$

$b_K=\frac{2}{U}\underset{u=0}{\overset{U-1}{\Σ}}\mathrm{\Δ}y\left(u\right)sin\frac{2\mathrm{\π}Ku}{U}$

$c_K=\sqrt{a_K^2+b_K^2}$

其中，K为分解后序号，取值范围从0～U/2，a_K和b_K分别为K次谐波余弦和正弦函数的振幅，c_K表示K次谐波的振幅大小；K次谐波对应频率f_K为，

$f_K=\frac{K}{U\·\mathrm{\Δ}t}=\frac{K}{T}$

得到平台震颤在沿推扫方向的频谱曲线。

6.一种基于多光谱影像的卫星平台震颤检测系统，包括以下模块：

数据准备模块，用于对原始的多光谱影像进行预处理，得到未配准多光谱影像；

波段选择模块，用于选取用于检测的两个波段；

密集匹配模块，用于选取样本点，然后利用影像相关匹配得到像素级匹配结果，最后利用最小二乘匹配得到子像素精度的同名像点坐标，得到匹配同名点；

震颤分析模块，用于首先对匹配同名点逐点计算配准误差，然后进行逐行分析得到配准误差随成像行的变化，最后通过傅里叶变换对平台震颤时序数据进行频谱分析，获取震颤引起的频谱特征；

对匹配同名点计算配准误差，实现方式如下，

设步骤2选择的两个波段为M、N，设第u扫描行中匹配出V个点，第u扫描行中第k个同名像点的配准误差计算如下，

${\mathrm{\Δx}}_{u,k}=x_{M_{u,k}}-x_{N_{u,k}}-\mathrm{\Δ}x$

${\mathrm{\Δy}}_{u,k}=y_{M_{u,k}}-y_{N_{u,k}}-\mathrm{\Δ}y$

其中，Δx_u,k，Δy_u,k分别为沿CCD方向和沿推扫方向的配准误差，分别为第u行的第k个同名像点在M，N两个波段影像上的坐标，Δx,Δy为两个波段CCD之间在沿CCD和沿推扫方向的平移设计值；

进行逐行分析得到配准误差随成像行的变化，实现方式如下，

设第u扫描行最终匹配出V个点，统计沿CCD方向的配准误差平均值Δx(u)和沿推扫方向的配准误差平均值Δy(u)，

$\mathrm{\Δ}x\left(u\right)=\underset{k=1}{\overset{k=V}{\Σ}}{\mathrm{\Δx}}_{u,k}/V$

$\mathrm{\Δ}y\left(u\right)=\underset{k=1}{\overset{k=V}{\Σ}}{\mathrm{\Δy}}_{u,k}/V$

逐行统计得到沿CCD方向和沿推扫方向的配准误差平均值Δx(u)、Δy(u)随成像行相应成像时间的变化规律曲线。