一种基于FPGA的线阵CCD卫星影像星上几何纠正方法与流程

本发明涉及遥感图像处理领域，具体涉及一种基于fpga的线阵ccd卫星影像星上几何纠正方法。

技术背景

随着卫星遥感数据产品应用范围越来越广，对数据实时性的要求也越来越高。星上完成遥感产品的生产是未来的大趋势，如军事打击目标的实时跟踪和监测，洪水实时监测，森林火点跟踪等。在卫星图像预处理中，卫星影像几何纠正是最重要的计算步骤之一。星上几何纠正直接关系后续图像处理算法能否获得高精度计算结果。传统的几何纠正方法是对高分辨率卫星的成像系统进行实验室几何纠正、安装检校以及在轨几何纠正。实验室几何纠正和安装检校是在卫星发射之前进行参数的检校，由于卫星发射之后受到外界环境影响，原始的技术参数不再使用，人们开始研究在轨几何纠正方法，但是在轨几何纠正方法需要收集地面检校场区域的影像进行几何纠正，这类方法局限性比较大，降低了卫星数据的实时性和利用效率。

fpga是一种硬件可编程的芯片，它的出现解决了传统定制电路无法再编程的不足，又提供了数量丰富的逻辑资源，人们可以根据自己的需求在fpga上订制自己所需的处理芯片。

综上所述，有必要从硬件方面研究纠正方法，使得卫星影像在星上即可完成实时纠正。

技术实现要素：

本发明的目的在于提供一种基于fpga的线阵ccd卫星影像星上几何纠正方法。本发明可以精确、快速地完成星上几何纠正计算，解决了高分辨率线阵ccd卫星影实时星上几何纠正的问题。

为实现本发明之目的，采用以下技术方案予以实现：一种基于fpga的线阵ccd卫星影像星上几何纠正方法，其具体包括以下步骤：

一、外方位元素初值计算

初值计算部分主要包含模板影像选取模块、影像匹配模块、外方位元素初值计算模块，每个模块功能如下：

1、模板影像选取模块根据卫星实时的轨道与姿态测量数据计算卫星成像的地面区域中心近似的大地坐标，并与模板影像附带的控制点坐标进行比较，选出目标模板影像存入外部存储器；

2、影像匹配模块从随机存储器中读取模板影像和目标影像进行匹配，提取出像点坐标及其大地坐标，同时将控制点所在扫描行的轨道与姿态观测数据传输至外方位元素初值计算模块；

3、外方位元素初值计算模块利用匹配之后影像的若干定向线中三个姿态参数转换为地面坐标至参考坐标系的外方位角元素，利用其外方位直线元素和转换后的外方位角元素通过多项式插值模型得到该像点所在扫描行的外方位元素，作为平差计算的初值；

二、平差计算

平差计算部分基于单像空间后方交会算法采用一次多项式模型为外方位线元素引入常数项改正，为外方位角元素引入线性误差改正，主要包含参数计算模块、求解模块、阈值比较模块，每个模块功能如下：

1、参数计算模块主要完成外方位角元素的正、余弦值计算以及加减乘除的组合运算，通过三个cordic核实现正、余弦值计算模块，采用并行计算结构计算旋转矩阵，使用三个减法器计算控制点大地坐标与相应扫描行的外定向位置参数之间差值。

2、求解模块利用最小二乘平差原理进行外方位元素改正值的计算，采用乘累加运算单元(pe)实现矩阵的相乘；

pe实现乘累加运算的过程为：(a)将输入的a和b两个数据相乘；(b)reset信号将寄存器的值设置为0；(c)下一个时钟周期，信号en控制寄存器中存储的初值0输出；(d)0值与经过delay单元延迟1个时钟周期后到来的乘积项(a×b)相加，然后经过寄存器延迟1个时钟周期后输出；(e)第一次相加得到的结果又作为加法器下一次运算的输入数据，与此时到来的乘积项相加，然后又经过寄存器后输出；(f)如此反复，经过2n(n为控制点的数量)次的乘法和加法计算后，就完成了一个相应矩阵元素值的计算。a^t和l相乘、(a^ta)^-1和(a^tl)的相乘可以采用类似的计算结构，只需要分别使用9个pe单元来完成计算。

为了减少fpga的逻辑资源的使用，根据矩阵a^ta是对称正定矩阵，利用ldlt分解法对法方程系数矩阵b进行分解求逆，求解法方程，输出外方位元素的改正数；

ldlt分解法具体实现过程包括：(1)利用数据选择器将矩阵b的元素按列方向构建为向量形式；(2)进行ldlt分解运算，得到单位下三角矩阵l和对角矩阵d的元素；(3)对矩阵l和d进行求逆，并对矩阵l^-1进行转置：a)利用数据选择器将单位下三角矩阵l的元素重排序，分别得到沿主对角线方向排列的多个向量数据；b)采用脉动阵列结构对矩阵l进行求逆，得到单位下三角矩阵l^-1；c)对求逆结果进行重排序，得到矩阵(l^-1)^t沿行方向的向量数据；d)利用数据选择器将对角矩阵d的对角线上所有元素构建为一个行向量，并通过一个除法器依次求解每个元素的倒数，完成d的逆矩阵的元素的求解；(4)将矩阵(l^-1)^t和d^-1的向量数据分别相乘得到上三角矩阵(l^-1)^td^-1，再将矩阵(l^-1)^td^-1与l^-1的向量数据分别相乘完成对称矩阵b^-1的对角线及其上方的元素的求解，然后利用12个数据选择器将矩阵元素构建为行向量，从而得到矩阵b^-1按行方向的9个行向量；

3、阈值比较模块判断外方位元素的改正数是否符合设定的阈值范围，符合则在一次平差之后输出结果，不符合则进行二次平差直至符合为止；

附图说明

图1为本发明实施例fpga硬件结构示意图。

图2为本发明实施例流程示意图。

图3为本发明实施例平差计算fpga硬件结构示意图

图4为本发明实施例参数计算子模块内部架构图

图5为本发明实施例乘累加运算单元(pe)计算结构图

图6为本发明实施例矩阵a^ta和a^tl相乘结构示意图

图7为本发明实施例矩阵b的ldlt分解与求逆的fpga实现结构

具体实施方式

下面结合附图和实施例对本发明的步骤做进一步的详细描述。

实施例

如图1所示，本发明主要应用于基于fpga的线阵ccd卫星影像星上几何纠正，以三线阵ccd相机moms-2p为例，其硬件方面，主要分为外方位元素初值计算和平差计算两部分。

而基于上述硬件结构，如图2所示，本发明实现流程如下：

一、外方位元素初值计算

初值计算部分主要包含模板影像选取模块、影像匹配模块、外方位元素初值计算模块；

1、模板影像选取模块根据卫星实时的轨道与姿态测量数据计算卫星成像的地面区域中心近似的大地坐标，并与模板影像附带的控制点坐标进行比较，选出目标模板影像存入外部存储器；

2、影像匹配模块从随机存储器中读取模板影像和目标影像进行匹配，提取出像点坐标及其大地坐标，同时将控制点所在扫描行的轨道与姿态观测数据传输至外方位元素初值计算模块；

3、外方位元素初值计算模块利用匹配之后影像的若干定向线中三个姿态参数转换为地面坐标至参考坐标系的外方位角元素，利用其外方位直线元素和转换后的外方位角元素通过多项式插值模型得到该像点所在扫描行的外方位元素，作为平差计算的初值；

二、平差计算

平差计算部分针对moms-2p影像的摄影测量处理采用的坐标系统，由参考坐标系到地面坐标系所构建的共线方程的形式如下：

式中，(xr,yr,zr)为影像参考坐标系下的像点坐标，(xsi,ysi,zsi)为历元t时刻相机曝光中心在wgs-84坐标系中的坐标，(x,y,z)为像点相应的地面点坐标；rij(i,j＝1,2,3)为矩阵的元素，为采用转角系统时影像的外方位角元素构成的地面坐标系到影像参考坐标系的旋转矩阵；

为了避免引入过多的未知参数所带来的参数间的强相关性问题，采用一次多项式模型表示moms-2p影像的外方位元素误差，为外方位线元素引入常数项改正，为外方位角元素引入线性误差改正，得到如下的误差方程：

写成矩阵形式为：

v＝ax－l,e(3)

式中，v＝[vxvy]^t；x＝[δa0δb0δc0δd0δe0δf0δd1δe1δf1]^t，为待求未知数向量；a为未知数x的系数矩阵，形式为：

其中为利用影像外方位元素初值按式(1)计算得到的影像坐标近似值；

根据误差方程式设计如图3中的参数计算模块所示包含七个子模块，用于得到法方程系数矩阵的元素，其中，通过三个cordic核实现正、余弦值计算模块；控制点大地坐标与外方位线元素相减模块使用了三个减法器实现；ti计算模块的内部将扫描行的行号乘以相邻扫描行之间的成像时间间隔，从而得到相应的成像时间；r11～r13计算模块采用了图4(a)所示的并行结构，用于计算用转角系统表示影像的外方位角元素时旋转矩阵的9个元素值；根据式(3)的说明部分给出的lx和ly的计算方法，可以设计出图4(b)所示的lx、ly计算模块的并行结构；计算结构的设计遵循以下两点要求：(1)利用延迟单元实现各参数计算的同步开始和结果的同时输出，采用并行结构来提高计算速度；(2)不同参数的计算若包含相同的输入数据或计算过程时，只给出唯一的实现，从而减少硬件资源的占用；按照同样的方法，可以设计得到如图4(c)所示的计算系数矩阵a的18个元素的并行结构；

根据最小二乘法求解误差方程设计了如图3中求解模块，设计了如图5所示的乘累加运算单元(pe)计算结构，根据pe设计了如图6实现矩阵a^t分别与矩阵a和l的相乘；

pe实现乘累加运算的过程为：(a)将输入的a和b两个数据相乘；(b)reset信号将寄存器的值设置为0；(c)下一个时钟周期，信号en控制寄存器中存储的初值0输出；(d)0值与经过delay单元延迟1个时钟周期后到来的乘积项(a×b)相加，然后经过寄存器延迟1个时钟周期后输出；(e)第一次相加得到的结果又作为加法器下一次运算的输入数据，与此时到来的乘积项相加，然后又经过寄存器后输出；(f)如此反复，经过2n(n为控制点的数量)次的乘法和加法计算后，就完成了一个相应矩阵元素值的计算。a^t和l相乘、(a^ta)^-1和(a^tl)的相乘可以采用类似的计算结构，只需要分别使用9个pe单元来完成计算。

根据矩阵a^ta是对称正定矩阵为了减少fpga的逻辑资源的使用，采用如下公式所示的ldlt分解法进行矩阵的求逆：

b^-1＝(ldl^t)^-1＝(l^t)^-1(ld)^-1＝(l^-1)^td^-1l^-1(5)

根据式(5)设计了如图7矩阵b的ldlt分解与求逆的fpga实现结构，主要步骤包含：

(1)利用数据选择器将矩阵b的元素按列方向构建为向量形式；

(2)进行ldlt分解运算，得到单位下三角矩阵l和对角矩阵d的元素；

(3)对矩阵l和d进行求逆，并对矩阵l^-1进行转置：a)利用数据选择器将单位下三角矩阵l的元素重排序，分别得到沿主对角线方向排列的多个向量数据；b)采用脉动阵列结构对矩阵l进行求逆，得到单位下三角矩阵l^-1；c)对求逆结果进行重排序，得到矩阵(l^-1)^t沿行方向的向量数据；d)利用数据选择器将对角矩阵d的对角线上所有元素构建为一个行向量，并通过一个除法器依次求解每个元素的倒数，完成d的逆矩阵的元素的求解；

(4)将矩阵(l^-1)^t和d^-1的向量数据分别相乘得到上三角矩阵(l^-1)^td^-1，再将矩阵(l^-1)^td^-1与l^-1的向量数据分别相乘完成对称矩阵b^-1的对角线及其上方的元素的求解，然后利用9个数据选择器将矩阵元素构建为行向量，从而得到矩阵b^-1按行方向的9个行向量；

利用如下公式进行法方程的求解：

x＝(a^ta)^-1(a^tl)(4)

阈值比较模块判断外方位元素的改正数x是否符合设定的阈值范围，符合则在一次平差之后输出结果，不符合则进行二次平差直至符合为止；

本实例应用在星上几何纠正上，如表1、2所示，表1是姿态参数转换的fpga计算结果与软件计算结果的对比，表2是一次平差的fpga计算结果与软件计算结果的对比。表1结果显示三组数据存在不同的差值，但是差值的大小均在可以忽略不计。表2结果显示两组数据中外方位线元素在厘米级甚至毫米级，外方位角元素的差值可以忽略不计，总体的差值在合理范围之内。

以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变形，因此，所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。

表1：

表2：