CN112464433B - 面向fpga硬件的递推最小二乘求解rfm模型参数优化算法 - Google Patents

Abstract

本发明公开了一种面向FPGA硬件的RLS求解RFM模型参数优化算法，主要过程包括：步骤一、确定RFM模型的形式；步骤二、在影像空间中建立“虚拟控制点”格网；步骤三、在地面空间中建立“虚拟控制点”格网；步骤四、根据所建立的“虚拟控制点”，利用RLS解算RFM模型参数；步骤五、精度评定。本发明以FPGA作为硬件加速平台，以Verilog硬件语言作为设计语言，在有限的硬件资源条件下，实现利用RLS解算RFM模型参数优化算法，满足卫星遥感领域对卫星影像正射纠正的时效性、便携性和小型化的要求。

Description

面向FPGA硬件的递推最小二乘求解RFM模型参数优化算法

技术领域

本发明涉及卫星遥感影像处理领域，特别地，是一种面向现场可编程门阵列(Field Programmable Gate Array，FPGA)的递推最小二乘(Recursive Least-Square，RLS)求解有理函数模型(Rational Function Model，RFM)参数优化算法。

技术背景

随着航天、传感等技术的不断发展，卫星遥感影像分辨率不断提高、波段数不断增加、重访周期不断缩短的同时，所获取的影像数据量也以惊人的速度在增大。然而，普通用户想要获取感兴趣的卫星影像专题产品，需要经过以下过程：卫星传感器获取数据，卫星下传数据到数据接收中心，再由地面的专业人员根据用户要求处理数据，最后再将处理结果分发给用户。这一过程是漫长的，用户至少需要等待1个月的时间才能得到影像专题产品，以至于用户对时效性的需求得不到切实的保障，而且导致大量的影像数据得不到有效的利用。研究表明，在轨数据处理技术是满足用户对数据加工、信息提取实时需求的关键技术之一，星上数据实时处理对实现智能对地观测卫星系统起着至关重要的作用。经过星上影像实时处理后的专题产品的数据量会大为减少，这不仅减小了数传压力，而且还能够让普通用户直接使用专题产品。相对于中央处理器(Central Processing Unit，CPU)和图形处理器(Graphic Processing Unit，GPU)，FPGA拥有更小的尺寸、更轻的重量、和更低的功耗。FPGA不仅解决了定制电路的不足，还克服了原有可编辑器件门电路数不足的缺点。因此，FPGA是星上数据实时处理的优选方案之一。

普通用户能够直接使用星上实时处理后的影像专题产品的先决条件之一是：在制作专题产品之前对卫星影像进行正射纠正。这是因为卫星在获取影像时因受到地形起伏、相机倾斜、地球曲率等因素的影响而使卫星影像产生投影误差和几何畸变。经过正射纠正的卫星影像不仅附带地图的几何信息，而且相对于一般的地图而言，正射影像能够更直观地表达信息，更新更容易。正射影像因具有可读性强的丰富信息，已经广泛的应用于各个领域：例如，应急救灾、军事侦查、国民经济建设、“数字城市”的建设、以及土地确权等。

近几十年来，学者们已经提出了众多的正射纠正算法模型。其中，RFM模型由于拥有数学形式简洁、不需向用户提供传感器参数、以及能够达到与严格几何模型一致的纠正精度等优点，已被广泛应用于卫星遥感影像的正射纠正中。然而，现有的RFM模型参数求解算法(如基于常规最小二乘的RFM模型参数求解算法)大多是基于地面平台来实现。这类基于地面平台的RFM模型参数求解算法无法满足对时效性有较高要求的应用场景(例如防灾减害、动态目标跟踪等)，其主要原因包括：(1)若要利用地面平台求解卫星遥感影像的RFM模型参数，须等待卫星遥感影像下传到地面接收站后方可处理，这一过程十分耗时；(2)地面平台多以串行方式处理数据，导致处理速度缓慢；(3)基于地面平台的RFM模型参数求解算法，即基于常规最小二乘的RFM模型参数求解算法由于需要对大型矩阵进行复杂的乘法和求逆运算，不仅会消耗大量FPGA硬件资源，还会降低求解RFM模型参数的速度。

因此，为了能够实现RFM模型参数的实时解算，一方面亟需研究适用于星上求解RFM模型参数的并行算法，另一方面需将特定的数据处理和应用模型在FPGA硬件上进行固化。

发明内容

针对上述现有技术存在的问题，本发明公开一种面向FPGA硬件的RLS求解RFM模型参数优化算法，通过结合FPGA并行计算的特性，对递推最小二乘算法进行优化，使RFM模型参数的解算能够在FPGA上实现硬件加速，以满足时效性、便携性和小型化的要求。

为了解决上述现有技术存在的问题，本发明公开的面向FPGA硬件的RLS求解RFM模型参数优化算法，包括以下步骤：

步骤一、确定RFM模型的形式；

步骤二、在影像空间中建立“虚拟控制点”格网；

步骤三、在地面空间中建立“虚拟控制点”格网；

步骤四、根据所建立的“虚拟控制点”，利用递推最小二乘法解算RFM模型参数；

步骤五、精度评定。

进一步讲，本发明所述的面向FPGA硬件的RLS求解RFM模型参数优化算法，步骤一中，RFM模型利用比例多项式来建立地面点的大地坐标(Lon,Lat,Hei)与对应像素坐标(i,j)之间的关系，一般形式为，

其中，NumL、DenL、NumS和DenS可分别由以下公式计算得到，

NumL＝a₁+a₂L+a₃P+a₄H+a₅LP+a₆LH+a₇PH+a₈L²+a₉P²+a₁₀H²+a₁₁PLH+a₁₂L³+a₁₃LP²+a₁₄LH²+a₁₅L²P+a₁₆P³+a₁₇PH²+a₁₈L²H+a₁₉P²H+a₂₀H³

DenL＝b₁+b₂L+b₃P+b₄H+b₅LP+b₆LH+b₇PH+b₈L²+b₉P²+b₁₀H²+b₁₁PLH+b₁₂L³+b₁₃LP²+b₁₄LH²+b₁₅L²P+b₁₆P³+b₁₇PH²+b₁₈L²H+b₁₉P²H+b₂₀H³

NumS＝c₁+c₂L+c₃P+c₄H+c₅LP+c₆LH+c₇PH+c₈L²+c₉P²+c₁₀H²+c₁₁PLH+c₁₂L³+c₁₃LP²+c₁₄LH²+c₁₅L²P+c₁₆P³+c₁₇PH²+c₁₈L²H+c₁₉P²H+c₂₀H³

DenS＝d₁+d₂L+d₃P+d₄H+d₅LP+d₆LH+d₇PH+d₈L²+d₉P²+d₁₀H²+d₁₁PLH+d₁₂L³+d₁₃LP²+d₁₄LH²+d₁₅L²P+d₁₆P³+d₁₇PH²+d₁₈L²H+d₁₉P²H+d₂₀H³

其中，a₁～a₂₀、b₁～b₂₀，c₁～c₂₀，d₁～d₂₀为有理多项式参数(Rational PolynomialCoefficients，RPCs)；P、L和H分别为大地坐标Lon、Lat、和Hei的正则化坐标；I和J分别为像素坐标i和j的正则化坐标。正则化坐标可分别由以下公式计算得到，

其中，Lonoff、LonScale、Latoff、LatScale、hoff、hScale、ioff、iScale、joff和jScale为正则化参数，可分别由以下公式计算得到，

LonScale＝max{[|max(Lon_VGCPs)-Lonoff|],[|min(Lon_VGCPs)-Lonoff|]}

LatScale＝max{[|max(Lat_VGCPs)-Latoff|],[|min(Lat_VGCPs)-Latoff|]}

hScale＝max{[|max(Hei_VGCPs)-hoff|],[|min(Hei_VGCPs)-hoff|]}

iScale＝max{[|max(i_VGCPs)-ioff|],[|min(i_VGCPs)-ioff|]}

jScale＝max{[|max(j_VGCPs)-Lonoff|],[|min(j_VGCPs)-joff|]}

其中，(Lon_VGCPs,Lat_VGCPs,Hei_VGCPs)为“虚拟控制点”的大地坐标，(i_VGCPs,j_VGCPs)为“虚拟控制点的”像素坐标，N为“虚拟控制点”个数。

本发明所述的面向FPGA硬件的RLS求解RFM模型参数优化算法，步骤二中，在影像空间中建立“虚拟控制点”格网的过程为：在原始卫星影像平面上以一定的间隔建立大小为m×n的格网，并记录格网线交点在原始卫星遥感影像上的行列号。这些行列号即为“虚拟控制点”的像素坐标。

本发明所述的面向FPGA硬件的RLS求解RFM模型参数优化算法，步骤三中，在地面空间中建立“虚拟控制点”格网的过程为：根据影像空间的“虚拟控制点”像素坐标，利用严格传感器模型计算“虚拟控制点”的大地坐标。在计算“虚拟控制点”的大地坐标之前，利用严格传感器模型计算卫星影像四个顶点对应的大地坐标，以确定卫星影像所覆盖地面的最大最小高程值，例如可通过美国地质调查局(United States Geological Survey，USGS)提供的全球DEM来确定。在确定高程的最大最小值后，把高程以一定的间隔分为若干层，然后再利用严格传感器模型计算“虚拟控制点”的像素坐标在各个高程条件下对应的大地坐标。为了防止在计算RFM模型参数时出现病态问题，高程的分层数量应大于4层。

本发明所述的面向FPGA硬件的RLS求解RFM模型参数优化算法，步骤四中，根据“虚拟控制点”，利用RLS解算RFM模型参数的过程为：

(1)随机产生初始化矩阵F₀和

大小分别为39×1和39×39；

(2)根据以下公式，依次利用“虚拟控制点”递推的计算RFM模型参数a₁～a₂₀和b₂～b₂₀(或c₁～c₂₀和d₂～d₂₀)，直至最后一个“虚拟控制点”；

(或

)

其中，

F＝(a₁,a₂,a₃,...,a₂₀,b₂,b₃,...,b₂₀)^T

K＝(c₁,c₂,c₃,...,c₂₀,d₂,d₃,...,d₂₀)^T

m_ρ+1＝[1,L,P,H,LP,LH,PH,L²,P²,H²,PLH,L³,LP²,LH²,L²P,P³,PH²,L²H,P²H,H³,-JL,-JP,-JH,-JLP,-JLH,-JPH,-JL²,-JP²,-JH²,-JPLH,-JL³,-JLP²,-JLH²,-JL²P,-JP³,-JPH²,-JL²H,-JP²H,-JH³]

其中，F_ρ和F_ρ+1为39×1的矩阵；

和

为39×39的矩阵；E为39×39的单位矩阵；m_ρ+1为1×39的矩阵；ρ为“虚拟控制点”序号；J_ρ+1、I_ρ+1为第ρ+1个“虚拟控制点”的正则化像素坐标；

为纯量。

(3)输出最终的RFM模型参数a₁～a₂₀和b₂～b₂₀(或c₁～c₂₀和d₂～d₂₀)。

本发明所述的面向FPGA硬件的RLS求解RFM模型参数优化算法，步骤五中，精度评定过程为：通过步骤二和步骤三来建立加密的检查点。在建立检查点时，格网密度及高程的分层数量通常是“虚拟控制点”的2倍。在确定RFM模型的参数后，根据检查点的大地坐标，利用RFM模型计算检查点的像素坐标；然后计算检查点的已知像素坐标与由RFM模型计算得到的像素坐标的偏差，并以此评定RFM模型的纠正精度。

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

本发明以FPGA作为硬件加速平台，以Verilog硬件语言作为设计语言，在有限的硬件资源条件下，实现利用RLS解算RFM模型参数优化算法，满足卫星遥感领域对卫星影像正射纠正的时效性、便携性和小型化的要求。

附图说明

图1面向FPGA硬件的递推最小二乘求解RFM模型参数优化算法的伪代码；

图2RLS求解RFM模型参数优化算法的FPGA硬件架构；

图3获取正则化坐标及m_ρ+1矩阵元素的FPGA硬件架构；

图4矩阵乘法的并行结构；

图5计算

的并行结构；

图6计算rtemp2＝m_ρ+1×temp1的并行结构；

图7计算temp＝temp1×invtemp2的并行结构；

图8计算rW＝E-temp×m_ρ+1的并行结构；

图9矩阵分块计算示意图；

图10计算

的并行结构。

具体实施方式

下面结合附图及具体实施例对本发明做进一步的说明，但下述实施例绝非对本发明有任何限制。

本发明公开的面向FPGA硬件的RLS求解RFM模型参数优化算法的伪代码如图1所示(以求解参数a₁～a₂₀，b₂～b₂₀为例)。其中，F₀和

为随机产生的初始矩阵，大小分别为39×1和39×39；F_ρ和F_ρ+1为39×1的矩阵；

和

为39×39的矩阵；E为39×39的单位矩阵；m_ρ+1为1×39的矩阵；ρ为“虚拟控制点”序号；J_ρ+1为第ρ+1个“虚拟控制点”的正则化像素坐标；temp1、temp和rS为39×1的矩阵；rtemp2、temp2、invtemp2、rJ、detaJ为纯量。

本发明采用Xilinx FPGA芯片，并利用Verilog硬件语言在Vivado2016.4软件平台上进行算法设计和程序编写。根据图1所示算法，设计了如图2所示的硬件架构，主要包括NORMALIZE模块、TEMP1模块、TEMP2模块、TEMP模块、UPDATE_W模块和UPDATE_S模块。其中，NORMALIZE模块主要用于获取正则化的大地坐标和像素坐标，并计算矩阵m_ρ+1中的元素term1～term39；TEMP1模块和TEMP模块分别用于计算矩阵temp1和矩阵temp中的元素；TEMP2模块主要用于计算标量temp2；UPDATE_W模块主要用于更新矩阵

UPDATE_S模块主要用于更新矩阵F_ρ+1。具体过程如下：

步骤1获取正则化的大地坐标和像素坐标，并计算矩阵元素。

在求解RFM模型参数前，需要利用优化后的公式1和公式2对“虚拟控制点”的大地坐标和像素坐标进行正则化，并利用优化后的公式3对正则化后的大地坐标和像素坐标进行相应地计算，以获取矩阵m_ρ+1中的元素。

其中，(Lon,Lat)和(i,j)分别为“虚拟控制点”的大地坐标中的经纬度和像素坐标；Hei为高程；rL、rP、rH和rJ为中间量；invLonScale、invLatScale、invhScale、invjScale分别为正则化参数LonScale、LatScale、hScale、jScale的倒数；Lonoff、Latoff、hoff、ioff和joff为正则化参数。

步骤2计算temp1矩阵中的元素。

为了能够快速获取矩阵相乘的结果，设计了如图4所示的矩阵乘法的并行结构，并通过FPGA对矩阵乘法进行了硬件加速。如图4所示，矩阵乘法的并行结构主要由相互独立的乘法累加器PE组成。由于PE之间不需要进行数据交换和通信，因此该矩阵乘法的并行结构具有很好的扩展性，在硬件资源充足的条件下，理论上可以实现任意维度的矩阵乘法。

下面以A×B＝C为例对矩阵乘法的并行计算原理进行了描述，其中矩阵A、B、C的大小均为39×39。如图4所示，A0101～A3939和B0101～B3939分别为矩阵A和矩阵B中的元素；C0101～C3939为矩阵C中的元素。为了能够获得正确的计算结果，需要在同一时刻并行的把矩阵A和矩阵B中的元素发送到对应的PE处理单元中，其中矩阵A和矩阵B分别以列优先和行优先的方式向对应的PE处理单元发送数据。例如，如图4所示，在t时刻，矩阵A的第1列数据，即A0101被广播到第1行的PE处理单元中(即PE0101～PE0139)，A0201被广播第2行的PE处理单元中(即PE0201～PE0239)，以此类推，A3901被广播第39行的PE处理单元中(即PE3901～PE3939)；而B矩阵的第1行数据，即B0101被广播到第1列的PE处理单元中(即PE0101～PE3901)，以此类推，B0139被广播到第39列的PE处理单元中(即PE0139～PE3939)。当上述数据完成乘法累加运算后，矩阵A的第2列数据和矩阵B的第2行将按照上述过程发送到对应的PE处理单元中进行乘法累加运算，矩阵C中的元素C0101～C3939将得到更新。当矩阵A最后一列的元素和矩阵B最后一行的元素按照上述过程完成乘法累加运算后即可得到最终的矩阵C。

TEMP1模块的主要功能为计算矩阵temp1中的元素。在计算矩阵

时，涉及到39×39矩阵乘以39×1矩阵。通过对图4中的矩阵乘法并行结构进行适当调整后，即可得到TEMP1模块的矩阵乘法并行结构(如图5所示)，其中W_ρ0101～W_ρ3939为矩阵

的元素、term1～term39为矩阵

的元素、temp1_0101～temp1_3901为矩阵temp1的元素。如图5所示，在计算矩阵temp1时，总共需要39个PE处理单元。

步骤3计算temp2。

TEMP2模块的主要功能为计算temp2。根据图1所示的算法可知，temp2可通过1+rtemp2计算得到，中间量rtemp2＝m_ρ+1×temp1所涉及的矩阵乘法类型为1×39矩阵乘以39×1矩阵。temp2和rtemp2皆为纯量。rtemp2可通过如图6所示的结构进行矩阵乘法运算得到。如图6所示，在计算中间变量rtemp2时只需要一个PE处理单元即可。在得到中间变量rtemp2后，只需再进行一次浮点数加法即可得到temp2。

步骤4计算temp矩阵中的元素。

TEMP模块的主要功能为计算temp矩阵中的元素。根据图1所示的算法可知，矩阵temp可由temp1×invtemp2计算得到，其中invtemp2为temp2的倒数。之所以利用invtemp2乘矩阵temp1，是为了避免对矩阵temp1中的每个元素进行除法运算。temp＝temp1×invtemp2为所涉及的矩阵乘法类型为39×1矩阵乘以纯量，其矩阵乘法结构如图7所示。

步骤5更新

矩阵中的元素。

UPDATE_W模块的主要功能为更新

矩阵中的元素。根据图1所示的算法可知，矩阵

的更新需要进行两次的矩阵乘法，即temp×m_ρ+1和

其中矩阵rW＝E-temp×m_ρ+1。

为了能够利用相互独立的乘法累加器来计算rW＝E-temp×m_ρ+1，设计了如图8所示的并行计算结构，其中PE'处理单元为针对rW＝E-temp×m_ρ+1所设计的第2类乘法累加器。如图8所示，用于计算矩阵rW的并行结构总共包含39个PE'处理单元，因此，每次运算都会同时得到矩阵rW一整行的元素。

综合考虑处理速度和FPGA硬件资源，采用图9所示的矩阵分块方法来进行

的矩阵乘法运算。如图9所示，目标矩阵

被划分为9块(每个分块大小为13×13)，并按图9中红色箭头所指方向进行计算。矩阵

的任意分块都可以通过如图10所示的并行结构计算得到。该并行结构总共包含了13×13个PE处理单元。

步骤6更新F_ρ+1矩阵中的元素。

根据图1所示的算法可知，在更新矩阵F_ρ+1时，所涉及的矩阵乘法主要包含rJ＝m_ρ+1×F_ρ和rS＝temp×detaJ，运算过程分别与步骤3和步骤4相同。因此，rJ＝m_ρ+1×F_ρ和rS＝temp×detaJ可分别通过如图6和图7所示的并行结构计算得到。在此不再赘述。

步骤7精度评定。

在建立了35×35×10的“虚拟控制点”(即在平面上建立35×35的格网点，并把高程为10层)后，分别利用常规最小二乘法和递推最小二乘法求解了三阶且分母不相同的RFM模型参数。然后建立了70×70×20的“虚拟检查点”(即在平面上建立70×70的格网点，并把高程为20层)，并随机选取其中的100个检查点，利用公式4～公式6计算均方根误差(RootMean Square Error,RMSE)，以评定RFM模型的纠正精度，结果如表1所示。

式中，N_VCP为选择的“虚拟检查点”个数，k_VCP＝1，2，3，…，N_VCP。RMSE_i、RMSE_j、RMSE_ijPlane分别为i方向、j方向和平面的RMSE。(i_LS,j_LS)和(i_RLS,j_RLS)分别为由常规最小二乘法和递推最小二乘法确定的RFM模型计算得到的检查点的像素坐标。

表1检查点像素坐标的偏差(单位：像素)

尽管上面结合附图对本发明进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨的情况下，还可以做出很多变形，这些均属于本发明的保护之内。

Claims (5)

1.面向FPGA的递推最小二乘RLS求解有理函数模型RFM参数优化算法，其特征在于，包括以下步骤：

步骤一、确定RFM模型的形式；

步骤二、在影像空间中建立“虚拟控制点”格网；

步骤三、在地面空间中建立“虚拟控制点”格网；

步骤四、根据所建立的“虚拟控制点”，利用RLS求解RFM模型参数a₁～a₂₀、b₂～b₂₀、c₁～c₂₀、和d₂～d₂₀；模型参数[c₁,c₂,…,c₂₀,d₂,d₃,…,d₂₀]的求解过程与模型参数[a₁,a₂,…,a₂₀,b₂,b₃,…,b₂₀]的过程相同；模型参数[a₁,a₂,…,a₂₀,b₂,b₃,…,b₂₀]的求解过程具体为：

分别随机初始化参数矩阵F和系数矩阵Q^-1为F₀和

其中矩阵F₀和

的大小分别为39×1和39×39；把N个“虚拟控制点”依次带入递推公式(1)和(2)，分别更新系数矩阵Q^-1和参数矩阵F为

和F_ρ+1，直至最后一个“虚拟控制点”；

其中，ρ为“虚拟控制点”序号，取值范围为[0,N-1]；参数矩阵F＝[a₁,a₂,…,a₂₀,b₂,b₃,…,b₂₀]^T；m_ρ+1为根据第ρ+1个“虚拟控制点”计算得到的系数向量，即：

L、P、H分别为第ρ+1个“虚拟控制点”的正则化大地坐标；J为第ρ+1个“虚拟控制点”的列像素坐标的正则化坐标；

为纯量；E为单位矩阵；

步骤五、精度评定。

2.根据权利要求1所述的面向FPGA的递推最小二乘RLS求解有理函数模型RFM参数优化算法，其特征在于，步骤四中，所述的利用RLS求解RFM模型参数a₁～a₂₀、b₂～b₂₀、c₁～c₂₀、和d₂～d₂₀的步骤为：

S1、分别随机初始化参数矩阵F、参数矩阵K、系数矩阵Q^-1、和系数矩阵O^-1为F₀、K₀、

和

矩阵大小分别为39×1、39×1、39×39和39×39；

S2、把N个“虚拟控制点”依次带入递推公式(3)、(4)、(5)、和(6)，分别更新系数矩阵Q^-1、系数矩阵O^-1、参数矩阵F、和参数矩阵K为

F_ρ+1、和K_ρ+1，直至最后一个“虚拟控制点”；

其中，ρ为“虚拟控制点”序号，取值范围为[0,N-1]；

为纯量；

为纯量；E为单位矩阵；参数矩阵F＝[a₁,a₂,…,a₂₀,b₂,b₃,…,b₂₀]^T；

参数矩阵K＝[c₁,c₂,…,c₂₀,d₂,d₃,…,d₂₀]^T；m_ρ+1和n_ρ+1为根据第ρ+1个“虚拟控制点”计算得到的系数向量，即：

其中，I、J为第ρ+1个“虚拟控制点”的正则化像素坐标；L、P、H分别为第ρ+1个“虚拟控制点”的正则化大地坐标；

S3、输出最终的RFM模型参数a₁～a₂₀、b₂～b₂₀、c₁～c₂₀、和d₂～d₂₀。

3.根据权利要求2所述的面向FPGA的递推最小二乘RLS求解有理函数模型RFM参数优化算法，其特征在于，在S2步骤中，为了便于FPGA硬件实现，公式(3)～公式(6)优化为以下计算过程：

s1、初始化：F₀、

s2、forρ＝0:N-1

s4、rtemp2＝m_ρ+1×temp1；

s5、temp2＝1+rtemp2；

s6、invtemp2＝1/temp2；

s7、temp＝temp1×invtemp2；

s8、rJ＝m_ρ+1×F_ρ；

s9、detaJ＝J_ρ+1－rJ；

s10、rS＝temp×detaJ；

s11、F_ρ+1＝F_ρ+rS；

s12、rW＝E-temp×m_ρ+1；

s14、end；

其中，F₀和

为随机产生的初始矩阵，大小分别为39×1和39×39；F_ρ和F_ρ+1为39×1的矩阵；

和

为39×39的矩阵；E为39×39的单位矩阵；m_ρ+1为1×39的矩阵；ρ为“虚拟控制点”序号；J_ρ+1为第ρ+1个“虚拟控制点”的正则化像素坐标；temp1、temp和rS为39×1的矩阵；rW为39×39的矩阵；rtemp2、temp2、invtemp2、rJ、detaJ为纯量。

4.根据权利要求3所述的面向FPGA的递推最小二乘RLS求解有理函数模型RFM参数优化算法，其特征在于，为了能够快速获取矩阵相乘的结果，设计了面向FPGA硬件的矩阵乘法并行结构；矩阵乘法的并行结构由相互独立的乘法累加器PE组成；PE之间不需要进行数据交换和通信。

5.根据权利要求4所述的面向FPGA的递推最小二乘RLS求解有理函数模型RFM参数优化算法，其特征在于，所述矩阵乘法并行结构计算过程如下：

对于大小均为39×39的矩阵A、B、C，A0101～A3939和B0101～B3939分别为矩阵A和矩阵B中的元素；C0101～C3939为矩阵C中的元素；为了能够正确计算A×B＝C，需要在同一时刻并行的把矩阵A和矩阵B中的元素发送到对应的PE处理单元中，其中矩阵A和矩阵B分别以列优先和行优先的方式向对应的PE处理单元发送数据；在t时刻，矩阵A的第1列数据，即A0101被广播到第1行的PE处理单元PE0101～PE0139中，A0201被广播第2行的PE处理单元PE0201～PE0239中，以此类推，A3901被广播第39行的PE处理单元PE3901～PE3939中；而B矩阵的第1行数据，即B0101被广播到第1列的PE处理单元PE0101～PE3901中，以此类推，B0139被广播到第39列的PE处理单元PE0139～PE3939中；当上述数据完成乘法累加运算后，矩阵A的第2列数据和矩阵B的第2行将按照上述过程发送到对应的PE处理单元中进行乘法累加运算，矩阵C中的元素C0101～C3939将得到更新；当矩阵A最后一列的元素和矩阵B最后一行的元素按照上述过程完成乘法累加运算后即可得到最终的矩阵C。

Images