留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

基于FPGA的P-H法星上解算卫星相对姿态

周国清 黄景金 舒磊

周国清, 黄景金, 舒磊. 基于FPGA的P-H法星上解算卫星相对姿态[J]. 武汉大学学报 ● 信息科学版, 2018, 43(12): 1838-1846. doi: 10.13203/j.whugis20180248
引用本文: 周国清, 黄景金, 舒磊. 基于FPGA的P-H法星上解算卫星相对姿态[J]. 武汉大学学报 ● 信息科学版, 2018, 43(12): 1838-1846. doi: 10.13203/j.whugis20180248
ZHOU Guoqing, HUANG Jingjin, SHU Lei. An FPGA-Based P-H Method On-Board Solution for Satellite Relative Attitude[J]. Geomatics and Information Science of Wuhan University, 2018, 43(12): 1838-1846. doi: 10.13203/j.whugis20180248
Citation: ZHOU Guoqing, HUANG Jingjin, SHU Lei. An FPGA-Based P-H Method On-Board Solution for Satellite Relative Attitude[J]. Geomatics and Information Science of Wuhan University, 2018, 43(12): 1838-1846. doi: 10.13203/j.whugis20180248

基于FPGA的P-H法星上解算卫星相对姿态

doi: 10.13203/j.whugis20180248
基金项目: 

国家自然科学基金 41431179

国家重点研发计划 2016YFB0502501

广西创新驱动发展专项(科技重大专项) GuikeAA18118038

详细信息

An FPGA-Based P-H Method On-Board Solution for Satellite Relative Attitude

Funds: 

The National Natural Science Foundation of China 41431179

the National Key Research and Development Program of China 2016YFB0502501

Guangxi Innovative Development Grand GuikeAA18118038

More Information
    Author Bio:

    ZHOU Guoqing, PhD, professor, specializes in the theories and methods of intelligent satellite observation systems and LiDAR technology. E-mail:gzhou@glut.edu.cn

    Corresponding author: HUANG Jingjin, PhD candidate. E-mail: jingjin_huang@tju.edu.cn
图(14) / 表(3)
计量
  • 文章访问数:  1511
  • HTML全文浏览量:  89
  • PDF下载量:  227
  • 被引次数: 0
出版历程
  • 收稿日期:  2018-08-02
  • 刊出日期:  2018-12-05

基于FPGA的P-H法星上解算卫星相对姿态

doi: 10.13203/j.whugis20180248
    基金项目:

    国家自然科学基金 41431179

    国家重点研发计划 2016YFB0502501

    广西创新驱动发展专项(科技重大专项) GuikeAA18118038

    作者简介:

    周国清, 博士, 教授, 主要从事智能卫星观测系统、激光雷达研究。gzhou@glut.edu.cn

    通讯作者: 黄景金, 博士生。jingjin_huang@tju.edu.cn
  • 中图分类号: P23

摘要: 针对目前星上遥感图像实时处理只能实现低级别算法的情况,提出了基于现场可编程门阵列(field-programmable gate array,FPGA)的P-H法星上相对姿态实时解算模型。该模型不仅避免了传统基于欧拉角的复杂三角函数计算与初值估算,还降低了迭代次数。试验选用FPGA(V7 xc7vx1140t)作为实时解算的硬件平台。在FPGA实现中,采用64位的浮点数据结构和串行/并行相结合策略;并采用LU(Lower-Upper)分解-分块算法实现矩阵求逆。试验结果表明,该模型的迭代次数比基于欧拉角的少了13次。该模型在FPGA和计算机的实现结果相差仅为5.0×10-14,加速度比为10。另外,该模型可广泛适用于实时性要求高的图像处理领域。

English Abstract

周国清, 黄景金, 舒磊. 基于FPGA的P-H法星上解算卫星相对姿态[J]. 武汉大学学报 ● 信息科学版, 2018, 43(12): 1838-1846. doi: 10.13203/j.whugis20180248
引用本文: 周国清, 黄景金, 舒磊. 基于FPGA的P-H法星上解算卫星相对姿态[J]. 武汉大学学报 ● 信息科学版, 2018, 43(12): 1838-1846. doi: 10.13203/j.whugis20180248
ZHOU Guoqing, HUANG Jingjin, SHU Lei. An FPGA-Based P-H Method On-Board Solution for Satellite Relative Attitude[J]. Geomatics and Information Science of Wuhan University, 2018, 43(12): 1838-1846. doi: 10.13203/j.whugis20180248
Citation: ZHOU Guoqing, HUANG Jingjin, SHU Lei. An FPGA-Based P-H Method On-Board Solution for Satellite Relative Attitude[J]. Geomatics and Information Science of Wuhan University, 2018, 43(12): 1838-1846. doi: 10.13203/j.whugis20180248
  • 随着各方面对卫星产品的多样性、实时性等需求不断提高,特别是军事行动、应急救灾等方面,国内外相继提出了智能地球观测卫星系统[1-3]、对地观测脑[4]等概念,这些系统都属于从数字地球发展到智慧地球的必然过程[5]。这些系统具有基于事件驱动,遥感数据星上实时处理,遥感产品实时分发的特点。而遥感数据具有空间、时间、光谱分辨率高等优点,但遥感卫星硬件资源有限、功耗低、体积小且实时性要求高,因此,遥感数据的星上实时处理是亟需解决的核心技术之一。

    针对星上实时处理的特点,一般采用基于现场可编程门阵列(field-programmable gate array,FPGA)[3-4]或者FPGA+数字信号处理(digital signal processing,DSP)[6]作为硬件架构。基于上述架构,目前能够实现星上实时处理的算法有分类[7]、数据压缩[8]、云检测[9-10]、图像复原[11]等。这类算法均属于遥感图像的预处理阶段(低级别算法),大部分算法采用定点数据结构即可完成。针对计算量大、时间复杂度高、浮点数据结构的算法,如三维重建、目标跟踪与识别、姿态解算等,目前还无法在星上进行实时处理。

    笔者在基于FPGA的局部特征检测与匹配[12-13]和基于FPGA的几何校正[14]的基础上,提出卫星相对姿态解算的FPGA实现。主要创新有:①在算法方面,采用单位四元数的P-H法相对姿态解算模型,该模型避免了旋转矩阵中大量三角函数的计算,同时减少了迭代次数,提高了迭代稳定性[15-16];②在矩阵求逆方面,直接的矩阵求逆方法会消耗大量的硬件资源和加大时间延迟。虽然传统的LU(lower-upper)分解算法降低了矩阵求逆的时间复杂度,但由于矩阵元素之间存在严重的依赖关系,并行度不大,在FPGA上实现难度较大。为解决上述问题,本文提出基于FPGA的LU分解-分块算法。该算法不仅降低了矩阵的维度和时间复杂度,更提高了算法本身的并行特性,适用于具有并行特性的FPGA硬件平台;③在数据结构和实现策略方面,为保证计算精度,采用了64位的双精度浮点数据结构和串/并行相结合的实现策略。这3个创新点保证了迭代结果精度高,处理速度快,硬件资源消耗少。

    • 令$ q = d + {\rm{i}}a + {\rm{j}}b + {\rm{k}}c$是单位四元数(i、j、k是虚数单位,dabc为实数,且$ {d^2} + {a^2} + {b^2} + {c^2} = 1$)。旋转矩阵R的四元数表达式为:

      $$ \mathit{\boldsymbol{R}} = \left[ {\begin{array}{*{20}{c}} {{d^2} + {a^2} - {b^2} - {c^2}}&{2\left( {ab - cd} \right)}&{2\left( {ac + bd} \right)}\\ {2\left( {ab + cd} \right)}&{{d^2} - {a^2} + {b^2} - {c^2}}&{2\left( {bc - ad} \right)}\\ {2\left( {ac - bd} \right)}&{2\left( {bc + ad} \right)}&{{d^2} - {a^2} - {b^2} + {c^2}} \end{array}} \right] $$ (1)

      R具有正交性,即RRT=I,求微分可得:

      $$ {\rm{d}}\left( {\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{R}}^{\rm{T}}}} \right) = {\rm{d}}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{R}}^{\rm{T}}} + \mathit{\boldsymbol{R}}{\rm{d}}{\mathit{\boldsymbol{R}}^{\rm{T}}} = 0 $$ (2)

      变换为:

      $$ {\rm{d}}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{R}}^{\rm{T}}} = - \mathit{\boldsymbol{R}}{\rm{d}}{\mathit{\boldsymbol{R}}^{\rm{T}}} = - {\left( {{\rm{d}}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{R}}^{\rm{T}}}} \right)^{\rm{T}}} $$ (3)

      式(3)为反对称阵,令:

      $$ {\mathit{\boldsymbol{S}}_w} = {\rm{d}}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{R}}^{\rm{T}}} $$ (4)

      则:

      $$ {\mathit{\boldsymbol{S}}_w}\mathit{\boldsymbol{R}} = {\rm{d}}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{R}}^{\rm{T}}}\mathit{\boldsymbol{R}} $$ (5)

      又因R-1=RT,得:

      $$ {\mathit{\boldsymbol{S}}_w}\mathit{\boldsymbol{R}} = {\rm{d}}\mathit{\boldsymbol{R}} $$ (6)

      R进行全微分,得:

      $$ {\rm{d}}\mathit{\boldsymbol{R}} = \frac{{\partial \mathit{\boldsymbol{R}}}}{{\partial d}}\Delta d + \frac{{\partial \mathit{\boldsymbol{R}}}}{{\partial a}}\Delta a + \frac{{\partial \mathit{\boldsymbol{R}}}}{{\partial b}}\Delta b + \frac{{\partial \mathit{\boldsymbol{R}}}}{{\partial c}}\Delta c $$ (7)

      则:

      $$ {S_w} = \left( {\frac{{\partial \mathit{\boldsymbol{R}}}}{{\partial d}}\Delta d + \frac{{\partial \mathit{\boldsymbol{R}}}}{{\partial a}}\Delta a + \frac{{\partial \mathit{\boldsymbol{R}}}}{{\partial b}}\Delta b + \frac{{\partial \mathit{\boldsymbol{R}}}}{{\partial c}}\Delta c} \right){\mathit{\boldsymbol{R}}^{\rm{T}}} $$ (8)

      对${d^2} + {a^2} + {b^2} + {c^2} = 1 $两边微分得:

      $$ \Delta d = - \frac{{a\Delta a + b\Delta b + c\Delta c}}{d} $$ (9)

      将式(9)代入式(8)中,消去Δd,得:

      $$ {\mathit{\boldsymbol{S}}_w} = \left[ {\begin{array}{*{20}{c}} 0&{{w_3}}&{ - {w_2}}\\ { - {w_3}}&0&{{w_1}}\\ {{w_2}}&{ - {w_1}}&0 \end{array}} \right] $$ (10)

      式中,

      $$ \begin{array}{*{20}{c}} {{w_1} = \frac{2}{d}\left[ {\left( {{d^2} + {a^2}} \right)\Delta a + \left( {ab + dc} \right)\Delta b + } \right.}\\ {\left. {\left( {ac - db} \right)\Delta c} \right]} \end{array} $$
      $$ \begin{array}{*{20}{c}} {{w_2} = \frac{2}{d}\left[ {\left( {ab - dc} \right)\Delta a + \left( {{d^2} + {b^2}} \right)\Delta b + } \right.}\\ {\left. {\left( {da + bc} \right)\Delta c} \right]} \end{array} $$
      $$ \begin{array}{*{20}{c}} {{w_3} = \frac{2}{d}\left[ {\left( {db + ac} \right)\Delta a + \left( {bc - da} \right)\Delta b + } \right.}\\ {\left. {\left( {{d^2} + {c^2}} \right)\Delta c} \right]} \end{array} $$

      令$\mathit{\boldsymbol{W}} = {[{w_1}\;\;{w_2}\;\;{w_3}]^{\rm{T}}} $,并整理得:

      $$ \mathit{\boldsymbol{W}} = \frac{2}{d}\left[ {\begin{array}{*{20}{c}} {{d^2} + {a^2}}&{ab + dc}&{ac - db}\\ {ab - dc}&{{d^2} + {b^2}}&{da + bc}\\ {db + ac}&{bc - da}&{{d^2} + {c^2}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\Delta a}\\ {\Delta b}\\ {\Delta c} \end{array}} \right] $$ (11)

      令:

      $$ \mathit{\boldsymbol{C}} = \frac{2}{d}\left[ {\begin{array}{*{20}{c}} {{d^2} + {a^2}}&{ab + dc}&{ac - db}\\ {ab - dc}&{{d^2} + {b^2}}&{da + bc}\\ {db + ac}&{bc - da}&{{d^2} + {c^2}} \end{array}} \right] $$ (12)

      式中,C有唯一逆矩阵:

      $$ {\mathit{\boldsymbol{C}}^{ - 1}} = \frac{1}{2}\left[ {\begin{array}{*{20}{c}} d&{ - c}&b\\ c&d&{ - a}\\ { - b}&a&d \end{array}} \right] $$ (13)

      则式(11)变为:

      $$ \left[ {\begin{array}{*{20}{c}} {\Delta a}\\ {\Delta {\rm{b}}}\\ {\Delta c} \end{array}} \right] = {\mathit{\boldsymbol{C}}^{ - 1}}\mathit{\boldsymbol{W}} = \frac{1}{2}\left[ {\begin{array}{*{20}{c}} {d{w_1} - c{w_2} + b{w_3}}\\ {c{w_1} + d{w_2} - a{w_3}}\\ { - b{w_1} + a{w_2} + d{w_3}} \end{array}} \right] $$ (14)

      从式(14)中可看到,P-H法的最大特点是通过求解中间参数W反算dabc的改正数,其本质上是一组基于反对称矩阵的求解方案。

    • 本文通过相对定向模型来确定两幅遥感图像在拍摄时刻的相对姿态。两幅有重叠遥感图像中的某一个点与摄像机在两个拍摄时刻位置形成一个面。利用共面约束条件,可以得到如下方程[19]

      $$ F = \left| {\begin{array}{*{20}{c}} {{B_x}}&{{B_y}}&{{B_z}}\\ {{X_t}}&{{Y_t}}&{{Z_t}}\\ p&q&r \end{array}} \right| = 0 $$ (15)

      式中,BxByBz为摄像机拍摄时的空间位置Bxyz方向上的分量;(Xt, Yt, Zt)为地面点在t时刻的像空间辅助坐标(该坐标系以t时刻的像空间坐标系为基准);pqr为同一地面点在t+1时刻的像空间辅助坐标。定义:

      $$ \left[ {\begin{array}{*{20}{c}} p\\ q\\ r \end{array}} \right] = \mathit{\boldsymbol{R}}\left[ {\begin{array}{*{20}{c}} {{X_{t + 1}}}\\ {{Y_{t + 1}}}\\ {{Z_{t + 1}}} \end{array}} \right] $$ (16)

      式中,$({X_{t + 1}}, {\rm{ }}{Y_{t + 1}}, {\rm{ }}{Z_{t + 1}}) $为t+1时刻的像空间坐标;Rt+1时刻的图像相对t时刻的图像的旋转矩阵,通过该矩阵可获得相对姿态参数φκω

      对式(15)进行微分求导得:

      $$ \begin{array}{*{20}{c}} {{\rm{d}}F = \frac{{\partial F}}{{\partial {B_y}}}{\rm{d}}{B_y} + \frac{{\partial F}}{{\partial {B_z}}}{\rm{d}}{B_z} + \frac{{\partial F}}{{\partial p}}{\rm{d}}p + }\\ {\frac{{\partial F}}{{\partial q}}{\rm{d}}q + \frac{{\partial F}}{{\partial r}}{\rm{d}}r} \end{array} $$ (17)

      式中,$ \frac{{\partial F}}{{\partial {B_y}}} = {Z_t}p - {X_t}r;\frac{{\partial F}}{{\partial {B_z}}} = {X_t}q - {Y_t}p;\frac{{\partial F}}{{\partial p}} = {B_y}{Z_t} - {B_z}{Y_t};\frac{{\partial F}}{{\partial q}} = {B_z}{X_t} - {B_x}{Z_t};\frac{{\partial F}}{{\partial r}} = {B_x}{Y_t} - {B_y}{X_t}$。

      根据R的正交性和P-H法的推导结果,对式(16)进行微分求导得:

      $$ \left[ {\begin{array}{*{20}{c}} {{\rm{d}}p}\\ {{\rm{d}}q}\\ {{\rm{d}}r} \end{array}} \right] = {\mathit{\boldsymbol{S}}_w}\mathit{\boldsymbol{R}}\left[ {\begin{array}{*{20}{c}} {{X_{t + 1}}}\\ {{Y_{t + 1}}}\\ {{Z_{t + 1}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {q{w_3} - r{w_2}}\\ {r{w_1} - p{w_3}}\\ {p{w_2} - q{w_1}} \end{array}} \right] $$ (18)

      将式(15)线性化,并结合式(17)、(18)得:

      $$ \begin{array}{*{20}{c}} {F = {F_0} + {\rm{d}}F = {F_0} + \left( {p{Z_t} - r{X_t}} \right){\rm{d}}{B_y} + }\\ {\left( {q{X_t} - p{Y_t}} \right){\rm{d}}{B_z} + \left( {r\left( {{B_z}{X_t} - {B_x}{Z_t}} \right) - } \right.}\\ {\left. {q\left( {{B_x}{Y_t} - {B_y}{X_t}} \right)} \right){w_1} + \left( {p\left( {{B_x}{Y_t} - {B_y}{X_t}} \right) - } \right.}\\ {\left. {r\left( {{B_y}{Z_t} - {B_z}{Y_t}} \right)} \right){w_2} + \left( {q\left( {{B_y}{Z_t} - {B_z}{Y_t}} \right) - } \right.}\\ {\left. {p\left( {{B_z}{X_t} - {B_x}{Z_t}} \right)} \right){w_3}} \end{array} $$ (19)

      式中,

      $$ \mathit{\boldsymbol{A}} = {\left[ {\begin{array}{*{20}{c}} {p{Z_t} - r{X_t}}\\ {q{X_t} - p{Y_t}}\\ {r\left( {{B_z}{X_t} - {B_x}{Z_t}} \right) - q\left( {{B_x}{Y_t} - {B_y}{X_t}} \right)}\\ {p\left( {{B_x}{Y_t} - {B_y}{X_t}} \right) - r\left( {{B_y}{Z_t} - {B_z}{Y_t}} \right)}\\ {q\left( {{B_y}{Z_t} - {B_z}{Y_t}} \right) - p\left( {{B_z}{X_t} - {B_x}{Z_t}} \right)} \end{array}} \right]^{\rm{T}}}; $$
      $$ \mathit{\boldsymbol{X}} = {\left[ {\begin{array}{*{20}{c}} {{\rm{d}}{B_y}}&{{\rm{d}}{B_z}}&{{w_1}}&{{w_2}}&{{w_3}} \end{array}} \right]^{\rm{T}}};\mathit{\boldsymbol{L}} = {F_0} - F。 $$

      由于观测值存在误差,式(19)的误差方程应变为:

      $$ \mathit{\boldsymbol{V}} = \mathit{\boldsymbol{AX}} - \mathit{\boldsymbol{L}} $$ (20)

      根据最小二乘法平差原理,可列出间接平差的法方程:

      $$ {\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{PAX}} = {\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{PL}} $$ (21)

      式中,定义每个像点获取的权重是相同的,即P=I,则未知数的向量解为:

      $$ \mathit{\boldsymbol{X}} = {\left( {{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}}} \right)^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{L}} $$ (22)

      将求得的w1w2w3代入式(14),求出单位四元数的改正数,得到更新值,再重复式(16)至式(22)的步骤,直至满足迭代条件为止。

    • 根据FPGA的并行处理特性和最小二乘算法的特点,提出如图 1所示的FPGA硬件架构。图 1中,实线框代表数值正在计算(结果不可调用),虚线框代表数值计算完成(结果可调用)。不同颜色箭头代表各变量的去向。T1-T0为一次迭代所需时间,在两次迭代之间(T1~T2),复位信号有效,各模块计算结果清零,并等待下一次迭代开始,每次迭代结果作为下一次迭代的输入,直至迭代结果满足规定阈值。当若干点对输入到系数模块中,采用串行处理,通过牺牲部分处理速度来减少大量的资源消耗,如pqrALATAATLX。对于计算复杂、时序要求高的模块,则采用并行处理,通过牺牲部分硬件资源来提高处理速度,如B-1。这种策略使处理速度和资源消耗达到平衡。下面详细介绍各个子模块。

      图  1  FPGA硬件架构

      Figure 1.  FPGA Architecture

    • FPGA在计算R时(图 2(a)),采用了3层结构,第一层为9个并行乘法,T1-T0为乘法运行时间; 第二层为10个并行加/减法,T2-T1为加/减法运行时间; 第三层为9个并行加/减法,运行时间为T3-T2,经过3层运算后,直接输出R。计算pqr时(图 2(b)),同样采用了3层结构,第一层为9个并行乘法,第二、三层都是3个并行加法。在串行计算中,只需调用一次该模块即可。

      图  2  Rpqr子模块

      Figure 2.  R and p, q, r Sub-modules

      图 3为系数AL的模块,A模块有4层结构,共有16个乘法器,8个减法器。L模块共有5层结构,其中有9个乘法器,5个加/减法器。L模块比A模块多一层加法运算,因此L运行时间比A多出一个加法器的时间。

      图  3  AL子模块

      Figure 3.  A and L Sub-modules

      图 4(a)为乘加子模块(MD_i),用于计算结构如“ $ {a_1}{b_1} + {a_2}{b_2} + \ldots + {a_i}{b_i}$ ”的结果。由于该架构使用9组点对(i=9),在ATAATL的计算中,均基于MD_9子模块运算。在串行执行中,ATA使用了5个该模块,ATL则使用了1个。

      图  4  ATAATL子模块

      Figure 4.  ATA and ATL Sub-modules

    • 在FPGA的矩阵求逆中,若直接通过伴随矩阵方式,会消耗大量硬件资源,增加计算时间,不利于实时性要求高的场景,因此,需要从另外角度进行矩阵求逆。一般的矩阵求逆方法有LU分解、LDLT(lower-diagonally-lower-transpose)分解及乔里斯基分解等。这些方法虽然能减少计算量,但计算速度有待提高,并行度不高。针对FPGA并行处理的特性,提出LU分解-分块算法[20]的矩阵求逆。由于ATA的大小为5×5,首先定义:

      $$ \mathit{\boldsymbol{B}} = \left[ {\begin{array}{*{20}{c}} {{b_{11}}}&{{b_{12}}}&{{b_{13}}}&{{b_{14}}}&{{b_{15}}}\\ {{b_{21}}}&{{b_{22}}}&{{b_{23}}}&{{b_{24}}}&{{b_{25}}}\\ {{b_{31}}}&{{b_{32}}}&{{b_{33}}}&{{b_{34}}}&{{b_{35}}}\\ {{b_{41}}}&{{b_{42}}}&{{b_{43}}}&{{b_{44}}}&{{b_{45}}}\\ {{b_{51}}}&{{b_{52}}}&{{b_{53}}}&{{b_{54}}}&{{b_{55}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{B}}_{11}}}&{{\mathit{\boldsymbol{B}}_{12}}}\\ {{\mathit{\boldsymbol{B}}_{21}}}&{{\mathit{\boldsymbol{B}}_{22}}} \end{array}} \right] $$ (23)

      式中,B11大小为3×3;B12为3×2;B21为2×3;B22为2×2。

      LU分解-分块算法为:

      $$ \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{B}}_{11}}}&{{\mathit{\boldsymbol{B}}_{12}}}\\ {{\mathit{\boldsymbol{B}}_{21}}}&{{\mathit{\boldsymbol{B}}_{22}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{L}}_{11}}}&{\bf{0}}\\ {{\mathit{\boldsymbol{L}}_{21}}}&{{\mathit{\boldsymbol{L}}_{22}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{U}}_{11}}}&{{\mathit{\boldsymbol{U}}_{12}}}\\ {\bf{0}}&{{\mathit{\boldsymbol{U}}_{22}}} \end{array}} \right] $$ (24)

      式中,B11=L11U11B12=L11U12B21=L21U11B22=L21U12+L22U22

      由式(24)计算得:

      $$ {\mathit{\boldsymbol{L}}_{11}} = \left[ {\begin{array}{*{20}{c}} 1&0&0\\ {{l_{21}}}&1&0\\ {{l_{31}}}&{{l_{32}}}&1 \end{array}} \right],{\mathit{\boldsymbol{U}}_{11}} = \left[ {\begin{array}{*{20}{c}} {{u_{11}}}&{{u_{12}}}&{{u_{13}}}\\ 0&{{u_{22}}}&{{u_{23}}}\\ 0&0&{{u_{33}}} \end{array}} \right] $$ (25)

      式中$ {u_{11}} = {b_{11}};{u_{12}} = {b_{12}};{u_{13}} = {b_{13}};{l_{21}} = {b_{21}}/{b_{11}}; = {b_{31}}/{b_{11}};$ ${u_{22}} = {b_{22}} - {l_{21}}{u_{12}};{u_{23}} = {b_{23}} - {l_{21}}{u_{13}}; $ $ _{32} = \left( {{b_{32}} - {l_{31}}{u_{12}}} \right)/{u_{22}};{u_{33}} = {b_{33}} - {l_{31}}{u_{13}} - {l_{32}}{u_{23}}$。

      计算L11U11后,可以得到:

      $$ {\mathit{\boldsymbol{U}}_{12}} = \mathit{\boldsymbol{L}}_{11}^{ - 1}{\mathit{\boldsymbol{B}}_{12}},{\mathit{\boldsymbol{L}}_{21}} = {\mathit{\boldsymbol{B}}_{21}}\mathit{\boldsymbol{U}}_{11}^{ - 1} $$ (26)

      令:

      $$ {{\mathit{\boldsymbol{B'}}}_{22}} = {\mathit{\boldsymbol{B}}_{22}} - {\mathit{\boldsymbol{L}}_{21}}{\mathit{\boldsymbol{U}}_{12}} = \left[ {\begin{array}{*{20}{c}} {{{b'}_{11}}}&{{{b'}_{12}}}\\ {{{b'}_{21}}}&{{{b'}_{22}}} \end{array}} \right] $$ (27)

      则:

      $$ {\mathit{\boldsymbol{L}}_{22}} = \left[ {\begin{array}{*{20}{c}} 1&0\\ {{{l'}_{21}}}&1 \end{array}} \right],{\mathit{\boldsymbol{U}}_{22}} = \left[ {\begin{array}{*{20}{c}} {{{u'}_{11}}}&{{{u'}_{12}}}\\ 0&{{{u'}_{22}}} \end{array}} \right] $$ (28)

      式中,$ {{u'}_{11}} = {{b'}_{11}};{{u'}_{12}} = {{b'}_{12}};{{l'}_{21}} = {{b'}_{21}}/{{u'}_{11}};{{u'}_{22}} = {{b'}_{22}} - {{l'}_{21}}{{u'}_{12}}$。

      因此,有:

      $$ {\mathit{\boldsymbol{B}}^{ - 1}} = {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{U}}_{11}}}&{{\mathit{\boldsymbol{U}}_{12}}}\\ {\bf{0}}&{{\mathit{\boldsymbol{U}}_{22}}} \end{array}} \right]^{ - 1}}{\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{L}}_{11}}}&{\bf{0}}\\ {{\mathit{\boldsymbol{L}}_{21}}}&{{\mathit{\boldsymbol{L}}_{22}}} \end{array}} \right]^{ - 1}} $$ (29)

      式中,$ {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{U}}_{11}}}&{{\mathit{\boldsymbol{U}}_{12}}}\\ {\bf{0}}&{{\mathit{\boldsymbol{U}}_{22}}} \end{array}} \right]^{ - 1}} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{U}}_{11}^{ - 1}}&{\mathit{\boldsymbol{M}}_{12}^{ - 1}}\\ {\bf{0}}&{\mathit{\boldsymbol{U}}_{22}^{ - 1}} \end{array}} \right], $

      $$ \mathit{\boldsymbol{M}}_{12}^{ - 1} = - \mathit{\boldsymbol{U}}_{11}^{ - 1}{\mathit{\boldsymbol{U}}_{12}}\mathit{\boldsymbol{U}}_{22}^{ - 1}; $$
      $$ {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{L}}_{11}}}&{\bf{0}}\\ {{\mathit{\boldsymbol{L}}_{21}}}&{{\mathit{\boldsymbol{L}}_{22}}} \end{array}} \right]^{ - 1}} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{L}}_{11}^{ - 1}}&{\bf{0}}\\ {\mathit{\boldsymbol{L}}_{21}^{ - 1}}&{\mathit{\boldsymbol{L}}_{22}^{ - 1}} \end{array}} \right], $$
      $$ \mathit{\boldsymbol{N}}_{21}^{ - 1} = - \mathit{\boldsymbol{L}}_{22}^{ - 1}{\mathit{\boldsymbol{L}}_{21}}\mathit{\boldsymbol{L}}_{11}^{ - 1}。 $$

      根据式(23)到式(29),LU分解-分块算法的流程见图 5

      图  5  LU分解-分块算法流程图

      Figure 5.  Flowchart of LU Decomposition Block Algorithm

      下面介绍各变量在FPGA中的实现。图 6L11-1U11-1的FPGA实现过程,B11为输入数据,由于L11-1U11-1分别属于上三角型矩阵和下三角型矩阵,因此只需求部分元素,其余元素为常量。图中“-x”为符号位取反模块,执行不需要时间。

      图  6  L11-1U11-1子模块

      Figure 6.  L11-1 and U11-1 Sub-modules

      图 7为式(27)的FPGA实现过程,其中分别使用了4个MD_2和MD_3模块。

      图  7  L21U12子模块

      Figure 7.  L21 and U12 Sub-modules

      图 8L22-1U22-1的FPGA实现过程,其中L21U12B22为输入,直接输出L22-1U22-1中的元素。

      图  8  L22-1U22-1子模块

      Figure 8.  L22-1 and U22-1 Sub-modules

      图 9为式(30)中M12-1N21-1的FPGA实现过程,其中输入分别为U11-1U12U22-1L22-1L21L11-1

      图  9  M12-1N21-1子模块

      Figure 9.  M12-1 and N21-1 Sub-modules

      图 10U11-1L11-1M12-1N21-1的FPGA实现过程,使用了多个并行MD_3和MD_2模块,属于计算B-1的中间结果。

      图  10  U11-1L11-1M12-1N21-1子模块

      Figure 10.  U11-1L11-1 and M12-1N21-1 Sub-modules

      图 11为FPGA实现B-1的最后一步,图 11(a)11(d)分别为B-1(1:3, 1:3)、B-1(1:3, 4:5)、B-1(4:5, 1:3)、B-1(4:5, 4:5)。图 6图 11为求解B-1的全部过程。

      图  11  B-1子模块

      Figure 11.  B-1 Sub-module

    • 图 12(a)为式(9)和式(14)的实现过程;图 12(b)为式(22)中X的实现过程;图 12(c)为更新后的abcdByBz。重复图 2图 12过程,直至Δa、Δb、Δc、Δd满足给定阈值。

      图  12  改正数模块

      Figure 12.  Correction Module

    • 使用Verilog HDL硬件语言在VIVADO软件实现上述各个模块,经过综合后,即可进行Modelsim仿真。图 13为顶层模块的仿真波形,其中$ L\_{X_t}、L\_{Y_t}、L\_{Z_t}$和$ R\_{X_t}、R\_{Y_t}、R\_{Z_t}$为输入的点对,d3a3b3c3为每次迭代后的更新结果。经过7次迭代后,it0_en信号拉低,表示迭代终止。图 14为LU分解-分块算法的仿真波形,每次迭代都需要进行矩阵求逆,图 14中显示的数值是第7次的求逆结果。

      图  13  迭代过程的仿真波形

      Figure 13.  Simulation Waveform of Iteration Processing

      图  14  LU分解-分块算法的仿真波形

      Figure 14.  Simulation Waveform of LU Decomposition Block Algorithm

    • FPGA参数为Xilinx, XC7VX1140T, 其主要资源有触发器(FF)1 424 000个,查找表(LUT)712 000个,块内存为1 880 kB,3 360个DSP单元等。

      PC机为Win10系统,CPU为Intel (R) Core (TM) i7-4770 CPU@3.40 GHz,内存为8 GB,软件为Dev-C++5.11。

      FPGA端和PC端输入参数为两幅图像的9个点对(表 1)。迭代终止条件为w1w2w3都大于或等于1.0×10-7,迭代结果见表 2。由表 2可知,PC1与PC2的最大偏差在ω,约为2.4×10-5;PC2与FPGA的最大偏差值约为5.0×10-14。基于P-H法的迭代次数比欧拉角的少13次。在计算速度方面,PC1与PC2的运行时间分别为2.269 ms和2.873 ms,FPGA的计算运行时间为0.308 ms(100 MHz),加速比约为10倍。FPGA资源消耗情况见表 3表 3中资源消耗最多的是DSP48,约为80.7%。由于整个模型属于数值运算,调用了大量的浮点运算IP核。其次是LUT和FF,百分比约为50%。所选的FPGA(V7 xc7vx1140tflg1930-1)满足整个模型所需的硬件资源。

      表 1  两幅图像点对(f=0.1 m)

      Table 1.  Point Correspondences of Two Images (f=0.1 m)

      点号 左图像 右图像
      行号 列号 行号 列号
      1 -37.104 -20.849 -43.682 -21.725
      2 -34.718 -21.669 -41.243 -22.652
      3 -32.359 -22.520 -38.839 -23.610
      4 -29.959 -23.372 -36.396 -24.567
      5 -27.631 -24.187 -34.030 -25.483
      6 -36.121 -18.675 -42.539 -19.529
      7 -33.601 -19.385 -39.958 -20.350
      8 -31.343 -20.141 -37.657 -21.206
      9 -28.945 -20.954 -35.218 -22.125

      表 2  结果对比

      Table 2.  Results Comparison

      参数 PC1 (欧拉角) PC2 (P-H法) FPGA (P-H法) |PC1-PC2| |PC2-FPGA|
      φ(picth)
      航向倾角
      0.043 373 3
      6 510 974
      0.043 373 3
      6 316 752
      0.043 373 3
      6 316 747
      1.9×10-9 5×10-14
      κ(yaw)
      旁向倾角
      0.044 063 6
      4 319 247
      0.044 063 6
      5 064 706
      0.044 063 6
      5 064 708
      7.5×10-9 2×10-14
      ω(roll)
      像片旋角
      0.021 523 4
      0 090 194
      0.021 547 6
      7 515 801
      0.021 547 6
      7 515 795
      2.4×10-5 5×10-14
      迭代次数 20 7 7 13 0
      运行时间/ms 2.269
      (3.40 GHz)
      2.873
      (3.40 GHz)
      0.308
      (100 MHz)
      / /

      表 3  FPGA硬件资源消耗

      Table 3.  Utilization of FPGA Resources

      资源 消耗情况 占总资源百分比/%
      FF 656 464 46.1
      LUT 382 344 53.7
      Memory LUT 45 596 16.1
      DSP48 2 712 80.7
    • 本文针对遥感图像复杂算法星上实时处理技术难点,提出了基于FPGA的P-H法星上相对姿态解算。对比试验发现,在算法层面,旋转矩阵的计算由基于欧拉角转换为基于单位四元数,采用P-H法平差模型转换为间接平差。在迭代处理中,直接取d=1和a=b=c=0,避免了初值估算和大量的三角函数计算,使迭代次数减少了13次;在矩阵求逆方面,提出LU分解-分块算法。该算法不仅使最小单元的矩阵维度由5×5降到了3×3,更提高了算法的并行性和处理速度;在FPGA方面,通过采用64位的双精度浮点数据结构,保证了计算结果的精度。如本方法在PC与FPGA计算结果偏差仅为5.0×10-14。另外,采用串行、并行相结合的策略,在系数(pqr, A, L, ATAATL)模块和改正数(X)模块中采用的是串行方式处理,在矩阵求逆模块和更新模块则采用并行处理,这种策略使得处理速度和资源消耗达到平衡。

参考文献 (20)

目录

    /

    返回文章
    返回