基于FPGA硬件的SPOT-6卫星影像RFM正射纠正方法与流程

基于fpga硬件的spot
‑
6卫星影像rfm正射纠正方法
技术领域
1.本发明涉及一种卫星遥感影像处理方法。特别是涉及一种基于fpga硬件的spot
‑
6卫星影像rfm正射纠正方法。


背景技术：

2.随着航天、传感等技术的不断发展，卫星遥感影像分辨率不断提高、波段数不断增加、重访周期不断缩短的同时，所获取的影像数据量也以惊人的速度在增大。然而，普通用户想要获取感兴趣的卫星影像专题产品，需要经过以下过程：卫星传感器获取数据，卫星下传数据到数据接收中心，再由地面的专业人员根据用户要求处理数据，最后再将处理结果分发给用户。这一过程是漫长的，用户至少需要等待1个月的时间才能得到影像专题产品，以至于用户对时效性的需求得不到切实的保障，而且导致大量的影像数据得不到有效的利用。研究表明，在轨数据处理技术是满足用户对数据加工、信息提取实时需求的关键技术之一，星上数据实时处理对实现智能对地观测卫星系统起着至关重要的作用。经过星上影像实时处理后的专题产品的数据量会大为减少，这不仅减小了数传压力，而且还能够让普通用户直接使用专题产品。相对于中央处理器(central processing unit，cpu)和图形处理器(graphic processing unit，gpu)，fpga拥有更小的尺寸、更轻的重量、和更低的功耗。fpga不仅解决了定制电路的不足，还克服了原有可编辑器件门电路数不足的缺点。因此，fpga是星上数据实时处理的优选方案之一。
3.普通用户能够直接使用星上实时处理后的影像专题产品的先决条件之一是：在制作专题产品之前对卫星影像进行正射纠正。这是因为卫星在获取影像时因受到地形起伏、相机倾斜、地球曲率等因素的影响而使卫星影像产生投影误差和几何畸变。经过正射纠正的卫星影像不仅附带地图的几何信息，而且相对于一般的地图而言，正射影像能够更直观地表达信息，更新更容易。正射影像因具有可读性强的丰富信息，已经广泛的应用于各个领域：例如，应急救灾、军事侦查、国民经济建设、“数字城市”的建设、以及土地确权等。
4.近几十年来，学者们已经提出了众多的正射纠正算法模型。其中，rfm模型由于拥有数学形式简洁、不需向用户提供传感器参数、以及能够达到与严格几何模型一致的纠正精度等优点，已被广泛应用于卫星遥感影像的正射纠正中。然而，现有的基于rfm模型的正射纠正方法大多是基于地面平台来实现。这类基于地面平台的rfm正射纠正方法无法满足对时效性有较高要求的应用场景(例如防灾减害、动态目标跟踪等)，其主要原因包括：(1)利用地面平台进行基于rfm模型的遥感影像正射纠正时，须等待卫星遥感影像下传到地面接收站后方可处理，这一过程十分耗时；(2)地面平台多以串行方式处理数据，导致处理速度缓慢。
5.因此，为了能够实时地进行卫星遥感影像的正射纠正，一方面亟需研究基于fpga硬件的卫星影像rfm正射纠正方法，另一方面需将特定的数据处理和应用模型在fpga硬件上进行固化。


技术实现要素：

6.本发明所要解决的技术问题是，提供一种在有限的硬件资源条件下，能够实现遥感影像的rfm正射纠正的基于fpga硬件的spot
‑
6卫星影像rfm正射纠正方法。
7.本发明所采用的技术方案是：一种基于fpga硬件的spot
‑
6卫星影像rfm正射纠正方法，包括如下步骤：
8.1)获取正射影像像素点的大地坐标；
9.2)利用rfm模型进行坐标转换，即将地面点的大地坐标转换为原始影像的像素行坐标和列坐标；
10.3)进行正射影像灰度插值；
11.4)重复步骤1)～步骤3)，直至得到正射影像上所有像素点的灰度值，从而获得经过正射纠正的影像。
12.本发明的基于fpga硬件的spot
‑
6卫星影像rfm正射纠正方法，以fpga作为硬件加速平台，以verilog硬件语言作为设计语言，在有限的硬件资源条件下，实现遥感影像的rfm正射纠正，满足卫星遥感领域对卫星影像正射纠正的时效性、便携性和小型化的要求。
附图说明
13.图1是基于fpga硬件的spot
‑
6卫星影像rfm正射纠正方法的fpga硬件架构；
14.图2是图1中getcord模块的fpga硬件架构；
15.图3是图2中getlonlat模块的fpga硬件架构；
16.图4是图2中interpolatehei模块的fpga硬件架构；
17.图5是图1中orthom模块的fpga硬件架构；
18.图6是图5中getterm模块的fpga硬件架构；
19.图7是图5中coeffmtermpe模块的fpga硬件架构；
20.图8是图5中getrowclm模块的fpga硬件架构；
21.图9是图5中interpolategrey模块的fpga硬件架构；
22.图10是由pc和fpga得到的spot
‑
6影像rfm正射纠正结果对比图；
23.图11是pc得到的纠正影像与fpga得到的纠正影像的差值图像。
具体实施方式
24.下面结合实施例和附图对本发明的基于丝杠和直线导轨的竖直方向波浪浮标检测装置及方法做出详细说明。
25.本发明的基于fpga硬件的spot
‑
6卫星影像rfm正射纠正方法，包括如下步骤：
26.1)获取正射影像像素点的大地坐标；包括：
27.设原始影像的4个顶点分别为p1、p2、p3和p4，它们对应的地面点p1、p2、p3和p4的大地坐标分别为(lon1,lat1)、(lon2,lat2)、(lon3,lat3)和(lon4,lat4)，那么正射影像覆盖的地面范围根据以下公式来确定，即：
[0028][0029]
在确定了正射影像的范围后，根据给定的经、纬度地面采样间距steplon和steplat来确定正射影像的行、列数量，如下公式所示：
[0030][0031]
式中，n
jorth
和n
iorth
分别为正射影像的行、列数量；
[0032]
设正射影像上任意一个像素点p
orth
的像素坐标为(i
orth
,j
orth
)，那么像素点p
orth
对应的地面点p
i
，i＝1,2,3,
…
,n
jorth
×
n
iorth
，地面点p
i
的大地坐标(lon,lat)通过如下公式得到：
[0033][0034]
式中，(i
lb
，j
lb
)为正射影像左下角像素点的像素坐标，即i
lb
＝0，j
lb
＝n
jorth
；lon为经度，lat为纬度；
[0035]
在得到地面点p
i
的大地坐标(lon,lat)后，根据以下公式计算地面点p
i
在dem中的行列坐标(i
d
,j
d
)，即：
[0036][0037]
其中，lond
min
和latd
min
分别为dem中经度和纬度的最小值；stlond和stlatd为dem的采样间距。在获取地面点p
i
在dem中的行列坐标(i
d
,j
d
)后，采用以下并行双线性插值方法计算地面点p
i
对应的高程值hei，
[0038][0039]
式中，h11和h12为中间变量；clmfdem、rowfdem为地面点p在dem(数字高程模型(digital elevation model)简称dem)中的行列坐标(i
d
,j
d
)的小数部分；h1、h2、h3、h4为四邻域像素对应的高程值；hei为最终的高程插值结果。
[0040]
2)利用rfm模型进行坐标转换，即将地面点的大地坐标转换为原始影像的像素行坐标和列坐标；
[0041]
是利用如下rfm模型把地面点的大地坐标坐标lon,lat,hei转换为原始影像的像素列坐标和行坐标i
orig
,j
orig
，
[0042][0043]
其中，ioff、iscale、joff和jscale为像素坐标正则化参数；多项式numl、denl、
nums和dens分别由以下公式计算得到，
[0044]
numl＝a1 a2l a3p a4h a5lp a6lh a7ph
[0045]
a8l2 a9p2 a
10
h2 a
11
plh a
12
l3 a
13
lp2
[0046]
a
14
lh2 a
15
l2p a
16
p3 a
17
ph2 a
18
l2h a
19
p2h a
20
h3[0047]
denl＝b1 b2l b3p b4h b5lp b6lh b7ph
[0048]
b8l2 b9p2 b
10
h2 b
11
plh b
12
l3 b
13
lp2
[0049]
b
14
lh2 b
15
l2p b
16
p3 b
17
ph2 b
18
l2h b
19
p2h b
20
h3[0050]
nums＝c1 c2l c3p c4h c5lp c6lh c7ph
[0051]
c8l2 c9p2 c
10
h2 c
11
plh c
12
l3 c
13
lp2
[0052]
c
14
lh2 c
15
l2p c
16
p3 c
17
ph2 c
18
l2h c
19
p2h c
20
h3[0053]
dens＝d1 d2l d3p d4h d5lp d6lh d7ph
[0054]
d8l2 d9p2 d
10
h2 d
11
plh d
12
l3 d
13
lp2
[0055]
d
14
lh2 d
15
l2p d
16
p3 d
17
ph2 d
18
l2h d
19
p2h d
20
h3[0056]
其中，a1～a
20
、b1～b
20
，c1～c
20
，d1～d
20
为rfm模型参数；p、l和h分别为大地坐标lon、lat、和hei的正则化坐标；正则化坐标分别由以下公式计算得到，
[0057][0058]
其中，lonoff、lonscale、latoff、latscale、hoff和hscale为大地坐标正则化参数。
[0059]
3)进行正射影像灰度插值；
[0060]
是采用如下并行双线性插值方法计算原始影像的像素坐标(i
orig
，j
orig
)的灰度值greyfpga，即：
[0061][0062]
式中，clmfgrey、rowfgrey分别为原始影像的像素坐标(i
orig
,j
orig
)的小数部分；gy1、gy2、gy3、gy4为原始影像的像素坐标(i
orig
,j
orig
)四邻域像素所对应的灰度值，gy11和gy12为中间变量；最后把灰度值greyfpga赋给像素坐标为(clmfpga,rowfpga)的正射影像像素点。
[0063]
4)重复步骤1)～步骤3)，直至得到正射影像上所有像素点的灰度值，从而获得经过正射纠正的影像。
[0064]
下面给出具体的实施例。
[0065]
本发明的实施例采用xilinx fpga芯片，并利用verilog硬件语言在vivado2016.4软件平台上进行算法设计和程序编写，设计了如图1所示的硬件架构，主要包括getcord模块和orthom模块。当信号isstart＝1'b1时，getcord模块将开始通过计算为orthom模块提供输入数据lonorth[63:0]、latorth[63:0]和heiorth[63:0]，它们分别为正射影像的经度、纬度和高程，其中lonorth[63:0]和latorth[63:0]可根据采样距离计算得到，而heiorth[63:0]可根据lonorth[63:0]和latorth[63:0]从dem内插得到。
[0066]
orthom模块主要负责影像的正射纠正，即进行坐标变换和灰度插值。当getcord模块开始产生大地坐标时，即信号isgetcord＝1'b1时，orthom模块将开始进行坐标变换的计算。在完成坐标变换后，orthom模块将通过行坐标rowgrey[63:0]和列坐标clmgrey[63:0]从rom中获取4个邻域的灰度值gy1[63:0]～gy4[63:0]，并进行双线性插值。
[0067]
本发明实施例的具体过程如下：
[0068]
步骤1获取正射影像像素的地面坐标。
[0069]
原始影像的4个顶点分别为p1、p2、p3和p4，它们对应的地面点p1、p2、p3和p4的大地坐标分别为(lon1,lat1)、(lon2,lat2)、(lon3,lat3)和(lon4,lat4)，那么正射影像覆盖的地面范围可根据公式(1)来确定，即：
[0070][0071]
在确定了正射影像的范围后，即可根据给定的地面采样间距steplon和steplat来确定正射影像的行、列数量，如公式(2)所示：
[0072][0073]
式中，n
jorth
和n
iorth
分别为正射影像的行、列数量。
[0074]
设正射影像上任意一个像点p
orth
的像素坐标为(i
orth
,j
orth
)，那么像点p
orth
对应的地面点p的大地坐标(lon,lat)可通过公式(3)得到，即：
[0075][0076]
式中，(i
lb
，j
lb
)为正射影像左下角顶点的像素坐标，即i
lb
＝0，j
lb
＝n
jorth
。
[0077]
在得到地面点p的大地坐标(lon,lat)后，根据公式(4)计算地面点p在dem中的行列坐标(i
d
,j
d
)，即：
[0078][0079]
其中，lond
min
和latd
min
分别为dem中经度和纬度的最小值；stlond和stlatd为dem的采样间距。在获取地面点p在dem中的行列坐标(i
d
,j
d
)后，采用公式(5)所示的并行双线性插值方法计算地面点p对应的高程值hei，
[0080][0081]
式中，h11和h12为中间量；clmfdem、rowfdem为地面点p在dem中的行列坐标(i
d
,j
d
)的小数部分；h1、h2、h3、h4为四邻域像素对应的高程值；hei为最终的高程插值结果。
[0082]
根据上述分析，为getcord模块设计了如图2所示的fpga硬件架构。该硬件架构主
要包含了用于获取大地坐标lonorth[63:0]和latorth[63]的getlonlat模块和用于内插高程heiorth[63:0]的interpolatehei模块。为了保证lonorth[63:0]和latorth[63:0]能够与heiorth[63:0]同步地输入到orthom模块，利用移位寄存器对lonorth[63:0]和latorth[63:0]进行了延迟。
[0083]
根据公式(3)，为getlonlat模块设计了如图3所示的fpga硬件架构。当开始信号isstart＝1'b1时，getlonlat模块将开始执行有效操作。经过一定的潜伏时间后，“fixed2float”将陆续有lonorth[63:0]和latorth[63:0]坐标数据输出，而且信号isgetlonlat将被拉高为高电平。当计数器countlat与计数器countlon分别等于正射影像的最大行列数orthnumlat和orthnumlon，即countlat＝orthnumlat
‑
1且countlon＝orthnumlon
‑
1时，说明正射影像上所有像素点的大地坐标均已获取完毕，因此，完成信号isdoneorth将被拉高为高电平。
[0084]
当interpolatehei模块检测到getlonlat模块发送而来的信号isgetlonlat为高电平时，interpolatehei模块将根据公式(4)和公式(5)对dem进行双线性插值。interpolatehei模块的fpga硬件架构如图4所示。在该硬件架构中，rowdem[63:0]和clmdem[63:0]均为64bits浮点数，为了便于分离它们的整数部分和小数部分，需要把64bits浮点数rowdem[63:0]和clmdem[63:0]分别转化为q40格式的64bits的定点数rowfx[63:0]和clmfx[63:0]。由此可得到rowfx[63:0]和clmfx[63:0]的整数部分分别为rowdem＝{rowfx[63:40],40'd0}和clmdem＝{clmfx[63:40],40'd0}，小数部分分别为rowfdem＝{24'd,rowfx[39:0]}和clmfdem＝{24'd,clmfx[39:0]}。一方面，利用rowdem[63:0]、clmdem[63:0]以及信号isgetneighborh，可从存储器中读取4个邻域地面点的高程h1[63:0]、h1[63:0]、h1[63:0]和h4[63:0]。另一方面，定点数rowfdem[63:0]和clmfdem[63:0]则重新转换为64bits的浮点数rowfdem[63:0]和clmfdem[63:0]，以便后续的计算。
[0085]
步骤2利用rfm模型进行坐标转换。
[0086]
orthom模块的fpga硬件架构如图5所示。orthom模块主要是利用rfm模型把正射影像像素点的大地坐标lonorth[63:0]、latorth[63:0]和heiorth[63:0]转换为原始影像的像素坐标rowfpga[63:0]和clmfpga[63:0]，然后再根据rowfpga[63:0]和clmfpga[63:0]内插出该像素点的灰度值greyfpga[63:0]。因此，该硬件架构主要由getterm模块、coeffmtermpe1～coeffmtermpe4模块、getrowclm模块和interpolategrey模块。各模块的功能如下：
[0087]
·
getterm模块
[0088]
getterm模块主要是通过公式(6)～公式(9)所示的4级运算来获取orthterm1[63:0]～orthterm20[63:0]。getterm模块采用流水线结构来获取orthterm1[63:0]～orthterm20[63:0]。根据公式(6)～公式(9)，为getterm模块设计了如图6所示的fpga硬件架构。由于getterm模块采用了流水线结构，为了保证能够同时输出orthterm1[63:0]～orthterm20[63:0]，因此需要利用延迟单元对orthterm1[63:0]～orthterm10[63:0]进行延迟。当信号isdoneorthnor为高电平时，说明getterm模块已开始向coeffmtermpe1等模块发送orthterm1[63:0]～orthterm20[63:0]的数据流。
[0089]
第1级：
[0090]
rl＝lonorth
‑
lonoff，rp＝latorth
‑
latoff，rh＝heiorth
‑
hoff
ꢀꢀꢀꢀ
(6)
[0091]
第2级：
[0092]
l＝rl
×
invlonscale，p＝rp
×
invlatscale，h＝rh
×
invhscale
ꢀꢀꢀꢀꢀ
(7)
[0093]
第3级：
[0094][0095]
第4级：
[0096][0097]
·
coeffmtermpe模块
[0098]
coeffmtermpe1～coeffmtermpe4模块具有相同的硬件结构。它们主要根据公式(10)～公式(15)来获取多项式numl、denl、nums和dens的值，即orthnuml[63:0]、orthdenl[63:0]、orthnums[63:0]和orthdens[63:0]。coeffmtermpe模块采用的是流水线结构。
[0099]
第1级：
[0100][0101]
第2级：
[0102][0103]
第3级：
[0104][0105]
第4级：
[0106][0107]
第5级：
[0108]
cmtl51＝cmtl41 cmtl42，cmtl35_delay2＝cmtl35_delay1，
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(14)
[0109]
第6级：
[0110]
resultscmt＝cmtl51 cmtl35_delay2
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(15)
[0111]
其中，cmtl11～cmtl42为中间变量；resultscmt为多项式的值。
[0112]
根据式公式(10)～公式(15)，设计了一个可复用的coeffmtermpe模块，其fpga硬件架构如图7所示。coeffmtermpe模块的fpga硬件架构为流水线结构，总共使用了20个浮点数乘法“mul”和19个浮点数加法“add”。当getterm模块发送到coeffmtermpe1～coeffmtermpe4模块的信号isdoneorthnor为高电平时，coeffmtermpe1～coeffmtermpe4模块开始执行计算，经过一定的潜伏期后，coeffmtermpe1～coeffmtermpe4模块将同时分别输出数据流orthnuml[63:0]、orthdenl[63:0]、orthnums[63:0]和orthdens[63:0]，而且信号isdonecmt(即图5中的信号isgetnuml)将被拉高为高电平。
[0113]
·
getrowclm模块
[0114]
在获得rfm模型多项式orthnuml[63:0]、orthdenl[63:0]、orthnums[63:0]和orthdens[63:0]的值后，根据公式(16)～公式(18)，即可获得原始影像像素点的像素坐标clmfpga[63:0]和rowfpga[63:0]。
[0115]
第1级：
[0116][0117]
第2级：
[0118]
tempclm＝norclm
×
iscale，temprow＝norrow
×
jscale(17)
[0119]
第3级：
[0120]
clmfpga＝tempclm ioff，rowfpga＝temprow joff(18)
[0121]
式中，iscale、jscale、ioff和joff为像素坐标的正则化参数；norclm和norrow为正则化的像素坐标；tempclm和temprow为中间量。
[0122]
根据公式(16)～公式(18)，采用流水线结构，为getrowclm模块设计了如图8所示的fpga硬件架构。在该硬件架构中，首先并行地计算norrow[63:0]和norclm[63:0]，然后再利用正则化参数iscale[63:0]、jscale[63:0]、ioff[63:0]和joff[63:0]反求原始影像的像素坐标rowfpga[63:0]和clmfpga[63:0]。getrowclm模块得到的rowfpga[63:0]和clmfpga[63:0]均为64bits浮点数。
[0123]
·
interpolategrey模块
[0124]
interpolategrey模块采用如公式(19)所示的并行双线性插值算法对像素点的灰度进行插值，因此需要对像素坐标rowfpga和clmfpga的整数部分和小数部分进行分离。根据整数部分rowgrey和clmgrey，可确定4个邻域像素点在影像上的位置分别为(clmgrey,rowgrey)、(clmgrey 1,rowgrey)、(clmgrey,rowgrey 1)和(clmgrey 1,rowgrey 1)，由此可得到它们的灰度值gy1、gy2、gy3和gy4。像素坐标rowfpga和clmfpga的小数部分分别为rowfgrey和clmfgrey。像素点(clmfpga,rowfpga)的灰度greyfpga可通过公式19计算得到，
即：
[0125][0126]
式中，gy11和gy12为中间量。
[0127]
根据公式(19)，为interpolategrey模块设计了如图9所示的fpga硬件架构。当getrowclm模块发送到interpolategrey模块的信号isdonegetrowclm为高电平时，说明getrowclm模块开始陆续地输出数据rowfpga[63:0]和clmfpga[63:0]，interpolategrey模块将开始执行灰度插值：
[0128]
首先，利用“float2fixed”ip核把像素坐标rowfpga[63:0]和clmfpga[63:0]分别转换为q40格式的64bits定点数rowfxg[63:0]和clmfxg[63:0]。
[0129]
然后，采用位拼接的方式分离定点数rowfxg[63:0]和clmfxg[63:0]的整数部分(如公式(20)所示)和小数部分(如公式(21)所示)，即，
[0130][0131][0132]
式中，rowgrey、clmgrey、rowfgreyq40、clmfgreyq40均为q40格式的64bits定点数。
[0133]
当信号isgetneighborg为高电平时，根据rowgrey[63:0]和clmgrey[63:0]的值即可从存储器中获取4个邻域像素点的灰度值，即gy1[63:0]、gy2[63:0]、gy3[63:0]和gy4[63:0]。
[0134]
rowfgreyq40[63:0]和clmfgreyq40[63:0]则将重新转换为64bits的浮点数rowfgrey[63:0]和clmfgrey[63:0]。当利用“fixed2float”ip核完成定点数到浮点数的转换后，mtvalid_clmfgrey信号将被拉高为高电平，以驱动公式20的计算，即进行灰度插值。当信号isgetgrey由低电平变为高电平时，说明将有像素点的灰度值被陆续地输出。
[0135]
实验：
[0136]
在windows10系统下，利用verilog硬件语言在vivado2016.4软件平台上完成了算法设计和程序编写，并利用spot
‑
6卫星影像数据进行了无控制点的正射纠正仿真实验。在进行基于fpga硬件的spot
‑
6卫星影像rfm正射纠正后，分别随机选取10个独立检查点(independent check points，icps)对纠正精度进行验证。如表1所示，对于spot
‑
6影像，icps的平面rmse分别为5.2399像素。此外，对fpga与pc的纠正影像进行了对比(如图10所示)，并分别对随机选取的50个检查点进行了像素坐标偏差统计(如表2所示)。
[0137]
由表1、表2和图10、图11可知，基于fpga硬件的spot
‑
6卫星影像rfm正射纠正方法有较大的纠正精度潜力。
[0138]
表1面向fpga硬件的rfm正射纠正算法的纠正精度(单位：像素)
[0139][0140]
表2检查点的像素坐标偏差统计(单位：像素)
[0141][0142]
尽管上面结合附图对本发明进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨的情况下，还可以做出很多变形，这些均属于本发明的保护之内。