CN110608723A - 一种基于fpga的实时遥感影像地理参考方法 - Google Patents

Abstract

本发明公开一种基于FPGA的实时遥感影像地理参考方法，包括如下步骤：S1、选定地理参考影像坐标构建控制点对坐标关系模型；S2、从控制点对坐标关系模型中选取相应的坐标参数实现待校正影像坐标(x_g,y_g)转换；S3、判断待校正影像坐标(x_g,y_g)是否为正数，如果满足进入下一步；否则结束；S4、判断待校正影像坐标(x_g,y_g)是否在选定地理参考图像范围，如果满足进入下一步；否则结束；S5、对待校正影像坐标(x_g,y_g)通过最邻近插值算法获得灰度值重采样输出具有地理坐标的遥感影像；该方法基于FPGA作为硬件实现，设计硬件语言为Verilog HDL语言，算法代码仿真有效，并实现了地理参考影像的实时输出。

Description

一种基于FPGA的实时遥感影像地理参考方法

技术领域

本发明涉及遥感影像处理领域，特别是一种基于FPGA的实时遥感影像地理参考方法。

背景技术

随着技术的不断进步，遥感影像数据越来越广泛的应用于自然灾害和定位系统中，这些系统需要快速、实时的提供地理坐标，然而，传统的遥感影像数据不带有地理坐标，而且传统的处理系统还不能满足实时性的要求，主要表现在：(1)一般经过地面接收站处理分发的初级遥感影像由于受到各种因素的影响，都存在一定的几何失真，主要表现为移位、旋转、缩放、仿射、弯曲和更高阶的弯曲，或者表现为像元相对地面的实际位置产生挤压、伸展、扭曲或者偏移。(2)一般的遥感影像都没有地理坐标。因此，高精度的几何纠正是遥感影像的在投入使用之前必须要解决的问题之一，研究各种几何纠正理论、开发相关的几何纠正算法具有重要的实际意义。

在摄影测量、遥感及扫描地图等得到的影像数据或激光扫描仪扫描得到的点云数据在应用于成图或地理信息系统时，需要将其和真实的地理坐标进行匹配，即图像对地球的校正，使图像像地图一样平面化，具有大地坐标的过程称之为地理参考，又称之为地理坐标的校准。地理参考是几何纠正重要手段之一，这一过程应该包括选择合适的数学变型模型、坐标转换和像元素灰度值重采样三个过程。

地理参考分为间接地理参考和直接地理参考，间接地理参考是利用地面控制点实现的地理坐标的匹配；而直接地理参考是是利用卫星上的惯性导航和全球卫星导航系统直接完成地理坐标的计算。但对于一般的用户而言，卫星的相关参数是不公开的，因此本设计主要是基于控制点的间接地理参考算法的实现。

一般的遥感影像的处理过程是，通过地面接收站接收下载数据后，进行影像的预处理形成不同等级的遥感数据，然后根据不同的用户需求进行处理完后发送给用户，这个过程是非常耗时的，显然在火灾、地震等自然灾害发生时，需要实时给出灾情的具体位置精确的地理坐标，因此一般系统不能满足实时性的要求，基于FPGA的实时地理参考系统是解决该问题的一种方法。

发明内容

针对现有技术存在的问题，本申请提供一种基于FPGA的实时遥感影像的地理参考算法，利用二次多项式方程建立控制点与原始影像坐标之间的关系，再逐点进行地理坐标的匹配，最后利用最近邻近插值法进行灰度值重采样，实现具有地理坐标的遥感影像制作；该方法具有实时性和低功耗性。

为了解决现有技术中存在的问题，本发明采用如下技术方案：

本发明公开一种基于FPGA的实时遥感影像地理参考方法，包括如下步骤：

S1、选定地理参考影像坐标构建控制点对坐标关系模型；

S2、从控制点对坐标关系模型中选取相应的坐标参数实现待校正影像坐标(x_g,y_g)转换；

S3、判断待校正影像坐标(x_g,y_g)是否为正数，如果满足进入下一步；否则结束；

S4、判断待校正影像坐标(x_g,y_g)是否在选定地理参考图像范围，如果满足进入下一步；否则结束；

S5、对待校正影像坐标(x_g,y_g)通过最邻近插值算法获得灰度值重采样输出具有地理坐标的遥感影像。

所述S1步骤中控制点对坐标关系模型建立过程：

1.1、通过如下二次多项式方程建立地理参考影像坐标(x,y)和大地坐标(X,Y)关系；

1.2、选取控制点建立与地理参考影像坐标(x,y)和大地坐标(X,Y)相应坐标矩阵模型；

(A^TPA)Δ_a＝A^TPL_x

其中A为控制点大地坐标矩阵，P为控制点的权重矩阵，一般取1，Δ_a为a_i系数矩阵，L_x为控制点待校正影像坐标。

所述S2步骤实现待校正影像坐标(x_g,y_g)转换过程：

2.1、通过待校正影像四个角点的坐标通过如下公式输出选定影像区域：

X_min＝min(X_ul,X_ur,X_lr,X_ll)

X_max＝max(X_ul,X_ur,X_lr,X_ll)

Y_min＝min(Y_ul,Y_ur,Y_lr,Y_ll)

Y_max＝max(Y_ul,Y_ur,Y_lr,Y_ll)

2.2、对选定影像区域进行扫描坐标的转换；

(1)通过如下公式计算选定影像区域的行(M)和列(N)坐标范围：

M＝(X_max-X_min)/X_GSD+1

N＝(Y_max-Y_min)/Y_GSD+1

其中,X_GSD和Y_GSD为别为地面采样间隔；

(2)根据扫描坐标按照如下公式获得地理参考影像的大地坐标：

X_g＝X_min+X_GSD(x'-1)

Y_g＝Y_max-Y_GSD(y'-1)

其中，x'和y'是扫描坐标(x'＝1,2,3,…,M；y'＝1,2,3,…,N)，(X_g,Y_g)为对应的大地坐标。

有益效果

本发明基于FPGA作为硬件实现，设计硬件语言为Verilog HDL语言，算法代码仿真有效。本设计方案具有运行速率快，功耗低的特性，能够实现遥感影像的地理参考的实时输出。

附图说明

图1基于FPGA的二次多项式的遥感影像实时地理参考结构框图；

图2二次多项方程模型解算结构框图；

图3原始图像坐标和输出图像坐标结构框图；

图4坐标转换和最近邻近插值算法流程图；

图5最近邻近插值结构图；

图6最近邻近插值算法电路原理图；

图7FPGA的功耗仿真图；

图8原始SPOT影像；

图9具有地理参考的SPOT影像；

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。

第一步：二次多项式系数解算

首先，建立待校正影像与参考影像之间的坐标转换方程，目前主要有严格成像模型、多项式模和有理函数模型，本专利中采用二次多项式模型，根据已知控制点数据对构造二次多项式方程，最后，利用最小二乘法循环求解最优二次多项式系数。

基于FPGA的二阶多项式的遥感影像实时地理参考实现结构框图如图1所示，主要包括二次多项式模型的建立、坐标转换和最近邻近插值三个主要部分，下面对三个方面分别详细地介绍。

二次多项方程如1式所示：

其中，(x,y)是待地理参考影像坐标，(X,Y)是相应同名点的大地坐标，a_i和b_i(i＝0,1,…, 5)是二次多项式的系数，总计有12个未知数。理论上说，只要有6个点对就可以求解未知方程系数，但必须保证6个点应该均匀的分布在影像上，这个一般是很大达到的，因此为了保证整个地理参考过程经度很高，因此选择了10个控制点对。设10个控制点对坐标表示为： (x₁,y₁,X₁,Y₁),(x₂,y₂,X₂,Y₂),…,(x₁₀,y₁₀,X₁₀,Y₁₀)，由于a_i和b_i的系数求解过程完全一样，这里以a_i系数求解来说明。将控制点对代入式1中转换成矩阵形式如式2所示：

根据最小二乘法原理，式(2)可以写成式(3)：

(A^TPA)Δ_a＝A^TPL_x (3)

其中A为控制点大地坐标矩阵，P为控制点的权重矩阵，一般取1，Δ_a为a_i系数矩阵，L_x为控制点待校正影像坐标，分别表示如下。

Δ_a＝[a₀ a₁ a₂ a₃ a₄ a₅]^T，L_x＝[x₁ x₂ x₃ x₄ x₅ …x₁₀]^T。

则,式(3)中a_i系数可以通过式(4)求得。

Δ_a＝(A^TA)^-1(A^TL_x) (4)

同理可求得Δ_b＝(A^TA)^-1(A^TL_y)的值。

二次多项式系数解算结构框图如图2所示。当时系数a₀，…,a₅和b₁，…,b₅求得后，则方程(1)就变成了参考影像坐标(x,y)和大地坐标(X,Y)对应方程的二元二次方程。每给出一对大地坐标(X,Y)就可以求得一组惟一的影响坐标(x,y)，实现了坐标映射，具体映射过程如下。

第二步：坐标转换；

当完成二次多项式系数求解后，就可以按照多项式方程建立参考影像与待校正影像的同名点之间坐标的关系，然后逐像素点进行坐标转换，最终将实现影像坐标转换为大地最标，使得原始图像既具有图像坐标，又具有大地坐标。

(1)输出影像的区域

为了实现坐标转换，首先估算校正影像的输出范围，设待校正影像四个角点的坐标为： (x_ul,y_ul),(x_ur,y_ur),(x_lr,y_lr)和(x_ll,y_ll)(其中，第一个字母：l代表lower，下边；u,代表up,上边。第二个字母：r代表right，右边；l代表left，左边)带入式(1)中可以求得四组坐标(X_ul,Y_ul), (X_ur,Y_ur),(X_lr,Y_lr),和(X_ll,Y_ll)，分别选出坐标中最小和最大的坐标作为校正后输出影像的范围坐标，即：

各坐标的位置关系如图3所示，其中3(a)为原始图像四个顶点的坐标图，3(b)为矫正后影像的输出范围坐标。

(2)输出影像与扫描坐标的转换

校正影像的行号(M)和列(N)号坐标范围由式(6)求得：

式中，X_GSD和Y_GSD为别为地面采样间隔。

地理参考影像的大地坐标就可以通过式(7)由扫描坐标求得。

其中，x'和y'是扫描坐标(x'＝1,2,3,…,M；y'＝1,2,3,…,N)，(X_g,Y_g)为对应的大地坐标。将求得的大地坐标(X_g,Y_g)带入式代入式(1)中就可以求得相应待校正影像的坐标(x_g,y_g)，下一步就是通过插值算法将大地坐标(X_g,Y_g)和待校正影像的坐标(x_g,y_g)建立一一对应关系，并将待校正影像最标的灰度值赋予给相应的大地坐标。

基于FPGA的并行坐标转换与最近邻近插值流程图如图4所示。

第三步：最近邻近插值。

在二次多项式方程进行坐标转换中，有时计算出来的最标值不可能恰好落在待校正影像的网格上，为了保证说有计算出来的坐标落在网格线上，就必须将不在网格线上的点近似移动到与其最邻近的网格上，这一过程是通过插值运算来完成的，本文采用最近邻近插值算法。所以，可以将第二步和第三步合在一起，同时完成。

通过二次多项式方程建立影像坐标与大地坐标之间的坐标转换关系，利用已知的地面控制点数据对求解方二次多项式系数，理论上有6对分布均匀的坐标点对就可以完成二次多项式系数的解算，但这在实际上是不可能找到非常准确的6个坐标点的，因此，实际采取的坐标点对都要大于理论点对数，本专利中为了折衷FPGA有限的资源和数据解算的精度，采用 10对坐标点对来实现。

由于式(7)解算出的待校正影像坐标(x_g,y_g)可能不是整数，这时需要进行插值，本专利中采用最邻近插值法(Nearest Neighbor Interpolation，简称NNI)。NNI算法原理如图5所示。设求得待校正影像的坐标为(i+u,j+v)(其中i,j为整数部分，u,v表示小数部分)，与其最相邻的4个像素坐标分别(i,j)，(i+1,j)，(i,j+1)和(i+1,j+1)。则，NNI算法可用式(8)表示。

其中，I(i+u,j+v)，I(i,j)，I(i+1,j)，I(i,j+1)和I(i+1,j+1)表示相应坐标(i,j)， (i+1,j)，(i,j+1)和(i+1,j+1)的灰度值。

基于FPGA的最邻近插值电路图如图6所示。

实验：

在Windows系统下使用Xilinx Vivado2014环境进行硬件语言设计，通过时序仿真与 MATLAB生成校正影像数据对比，仿真数据一致，图8为原始的SPOT TM遥感图，图9为校正后的带有地理参考信息的遥感图，表1为校正模型求解和坐标转换以及最近临近插值算法的FPGA资源消耗。图7为系统的功耗仿真图。综上所述，本系统实现了实时低功耗的二阶多项式遥感影像的地理参考算法。

表1(a).坐标转换资源消耗列表

Claims (3)

1.一种基于FPGA的实时遥感影像地理参考方法，其特征在于，包括如下步骤：

S1、选定地理参考影像坐标构建控制点对坐标关系模型；即：

Δ_a＝(A^TA)^-1(A^TL_x)

Δ_b＝(A^TA)^-1(A^TL_y)；

S2、从控制点对坐标关系模型中选取相应的坐标参数实现待校正影像坐标(x_g,y_g)转换；

S3、判断待校正影像坐标(x_g,y_g)是否为正数，如果满足进入下一步；否则结束；

S4、判断待校正影像坐标(x_g,y_g)是否在选定地理参考图像范围，如果满足进入下一步；否则结束；

S5、对待校正影像坐标(x_g,y_g)通过最最近邻近插值算法获得灰度值重采样输出具有地理坐标的遥感影像。

2.根据权利要求1所述的一种基于FPGA的实时遥感影像地理参考方法，其特征在于，所述S1步骤中控制点对坐标关系模型建立过程：

1.1、通过如下二次多项式方建立地理参考影像坐标(x,y)和大地坐标(X,Y)关系；

1.2、选取控制点建立与地理参考影像坐标(x,y)和大地坐标(X,Y)相应坐标矩阵模型；

(A^TPA)Δ_a＝A^TPL_x

其中A为控制点大地坐标矩阵，P为控制点的权重矩阵，一般取1，Δ_a为a_i系数矩阵，L_x为控制点待校正影像坐标。

3.根据权利要求1所述的一种基于FPGA的实时遥感影像地理参考方法，其特征在于，所述S2步骤实现待校正影像坐标(x_g,y_g)转换过程：

2.1、通过待校正影像四个角点的坐标通过如下公式输出选定影像区域：

X_min＝min(X_ul,X_ur,X_lr,X_ll)

X_max＝max(X_ul,X_ur,X_lr,X_ll)

Y_min＝min(Y_ul,Y_ur,Y_lr,Y_ll)

Y_max＝max(Y_ul,Y_ur,Y_lr,Y_ll)

2.2、对选定影像区域进行扫描坐标的转换；

(1)通过如下公式计算选定影像区域的行(M)和列(N)坐标范围：

M＝(X_max-X_min)/X_GSD+1

N＝(Y_max-Y_min)/Y_GSD+1

其中,X_GSD和Y_GSD为别为地面采样间隔；

(2)根据扫描坐标按照如下公式获得地理参考影像的大地坐标：

X_g＝X_min+X_GSD(x'-1)

Y_g＝Y_max-Y_GSD(y'-1)

其中，x'和y'是扫描坐标(x'＝1,2,3,…,M；y'＝1,2,3,…,N)，(X_g,Y_g)为对应的大地坐标。