CN114897676A - 一种无人机遥感多光谱图像拼接方法、设备及介质 - Google Patents

一种无人机遥感多光谱图像拼接方法、设备及介质 Download PDF

Info

Publication number
CN114897676A
CN114897676A CN202210299506.1A CN202210299506A CN114897676A CN 114897676 A CN114897676 A CN 114897676A CN 202210299506 A CN202210299506 A CN 202210299506A CN 114897676 A CN114897676 A CN 114897676A
Authority
CN
China
Prior art keywords
image
images
splicing
point
waveband
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202210299506.1A
Other languages
English (en)
Inventor
刘博�
孟庆民
张雅南
王铁
高有
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Vastriver Technology Co ltd
Original Assignee
Beijing Vastriver Technology Co ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beijing Vastriver Technology Co ltd filed Critical Beijing Vastriver Technology Co ltd
Priority to CN202210299506.1A priority Critical patent/CN114897676A/zh
Publication of CN114897676A publication Critical patent/CN114897676A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformation in the plane of the image
    • G06T3/40—Scaling the whole image or part thereof
    • G06T3/4038—Scaling the whole image or part thereof for image mosaicing, i.e. plane images composed of plane sub-images
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/50—Image enhancement or restoration by the use of more than one image, e.g. averaging, subtraction
    • G06T5/80—
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/37—Determination of transform parameters for the alignment of images, i.e. image registration using transform domain methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/40—Extraction of image or video features
    • G06V10/46—Descriptors for shape, contour or point-related descriptors, e.g. scale invariant feature transform [SIFT] or bags of words [BoW]; Salient regional features
    • G06V10/462—Salient features, e.g. scale invariant feature transforms [SIFT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
    • G06V10/74—Image or video pattern matching; Proximity measures in feature spaces
    • G06V10/75—Organisation of the matching processes, e.g. simultaneous or sequential comparisons of image or video features; Coarse-fine approaches, e.g. multi-scale approaches; using context analysis; Selection of dictionaries
    • G06V10/757—Matching configurations of points or features
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • G06T2207/10036—Multispectral image; Hyperspectral image
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20212—Image combination
    • G06T2207/20221—Image fusion; Image merging

Abstract

本发明实施例公开了一种无人机遥感多光谱图像拼接方法、设备及介质，所述方法包括：获取N帧单波段图像并进行波段配准校正以得到N帧多光谱图像，其中，每帧单波段图像包含M个波段的源图像；分别对所述N张多光谱图像进行波段融合得到N张灰度图像，基于所述N张灰度图像进行拼接以获得N个局部单应性矩阵集合；通过第n个局部单应性矩阵集合Dn分别对第n张单波段图像的第m波段的灰度图像进行图像配准，将配准后的相同波段图像进行拼接以获得M个单波段拼接图像，对所述M个单波段拼接图像进行融合处理，得到初步拼接图像；对所述初步拼接图像中存在的线段结构进行检验与矫正得到最终拼接结果。

Description

一种无人机遥感多光谱图像拼接方法、设备及介质
技术领域
本发明涉及图像处理领域。更具体地，涉及一种无人机遥感多光谱图像拼接方法、设备及介质。
背景技术
随着光谱成像技术的不断发展，配备更加轻便的光谱成像仪的无人机已经广泛地应用于地理遥感、农业遥感等方面。由于无人机飞行高度较低，每次采集的图像覆盖面积有限，为了获取整个场景更全面的信息，需要将多个图像拼接在一起，获得更大视角的全景图像。但是现有的拼接方法更多地针对全色图像或RGB图像，针对多光谱图像的拼接方法较少；而且在拼接过程中会存在扭曲和变形的情况，由于拼接图像数量较多，图像质量参差不齐使得目前大多数算法难以得到高效正确的拼接结果。
发明内容
本发明的目的在于提供一种无人机遥感多光谱图像拼接方法、设备及介质，以解决现有技术存在的问题中的至少一个。
为达到上述目的，本发明采用下述技术方案：
本发明第一方面提供了一种无人机遥感多光谱图像拼接方法，包括：
S10：获取N帧单波段图像并进行波段配准校正以得到N帧多光谱图像，其中，每帧单波段图像包含M个波段的源图像，其中，
N和M为正整数；
S20：分别对所述N张多光谱图像进行波段融合得到N张灰度图像，基于所述N张灰度图像进行拼接以获得N个局部单应性矩阵集合；
S30：通过第n个局部单应性矩阵集合Dn分别对第n张单波段图像的第m波段的灰度图像进行图像配准，将配准后的相同波段图像进行拼接以获得M个单波段拼接图像，对所述M个单波段拼接图像进行融合处理，得到初步拼接图像，其中，
n＝1、2、……、N；m＝1、2、……、M；
S40：对所述初步拼接图像中存在的线段结构进行检验与矫正得到最终拼接结果。
进一步地，所述S10包括：
S101：通过无人机多光谱照相机沿特定路线获取视频，所述视频包括N帧单波段图像，每帧单波段图像包含M个波段的源图像；
S102：利用SIFT算法分别对所述N帧单波段图像进行波段配准校正：通过不同尺寸的高斯核函数与当前单波段图像进行卷积，在不同的层级中对卷积之后的结果进行图像缩小以得到尺度空间，并在尺度空间中对源矩阵和网矩阵进行极值检测以得到特征点，同时移除特征点集中的边缘响应点，为关键点分配方向值；
根据所述方向值生成特征点描述子，对特征点进行匹配得到匹配点；
利用匹配点计算变换矩阵对当前单波段图像的M个波段的源图像进行图像变换和差值，得到波段配准后的多光谱图像并对M个波段的源图像进行更新。
进一步地，所述特征描述子包括特征点所在的尺度、位置以及方向。
进一步地，所述S20包括：
S201：分别对所述N张多光谱图像进行波段融合得到N张灰度图像；
S202：分别对所述N张灰度图像进行单波段配准；
S203：采用逼近投影变换算法进行所述N张灰度图像的拼接以得出对应的N个局部变换矩阵集合。
进一步地，所述S203包括：
S2031：分别对所述N张灰度图像进行特征提取和匹配得到预变换矩阵；
S2032：对灰度图像进行分块将图像均匀地划分成一个个小的网格，移动直线线性变换算法通过计算当前网格内的匹配点与网格中心的距离加权平方和误差来得到所述局部单应性变换矩阵集合Dn，其中，
Figure BDA0003564768620000021
K为当前灰度图像中的网格总数；
S2033：利用所述局部单应性变换矩阵集合对所述灰度图像进行变换，将变换后的图像与参考图像利用加权平均融合法进行图像融合与拼接。
进一步地，第n张灰度图像的第k个网格的局部单应性矩阵
Figure BDA0003564768620000022
的计算公式如式(1)所示，通过求解最小特征值对应的特征向量来估计局部变换矩阵
Figure BDA0003564768620000023
的值:
Figure BDA0003564768620000025
其中，ai为匹配点之间的线性映射函数，l为匹配点对数；d为匹配点与网格中心的距离，
Figure BDA0003564768620000031
为特征点所对应的权重，Wk为权重对应的对角矩阵，A为映射函数对应的映射矩阵形式，其中，
Figure BDA0003564768620000032
Figure BDA0003564768620000033
其中，xk和yk是第k个网格的中心点坐标，
Figure BDA0003564768620000034
和
Figure BDA0003564768620000035
是第k个网格内第i个匹配点的坐标，σ为高斯尺度因子，其范围在(0,1)之间。
进一步地，所述S40包括：
S401：进行线段结构的检测与提取，提取所述初步拼接图像中的线支撑区域，并通过所述线支撑区域得到线段特征，其中，
所述线支撑区域由有线段存在的像素区域组成；
S402：进行线段校正：利用最小二乘法在矩形区域拟合出一条线段，将所述线段的交点作为特征点，将其与矩形区域中心点的位置关系进行排序，根据变换后的特征点与原来的特征点之间的对应关系计算透视变换矩阵，将待校正的区域内所有的点进行映射，完成线段校正。
进一步地，所述S401包括：
S4011：利用每个像素点右下方的四个像素来计算该像素点对应的梯度以得到线段边缘的位置，像素点沿着x方向和y方向的梯度hx(x,y)和hy(x,y)计算公式如式(4)和式(5)所示：
Figure BDA0003564768620000036
Figure BDA0003564768620000037
其中h(x，y)表示坐标为(x，y)的像素点对应的灰度值；
梯度大小H(x,y)和方向θ(x,y)的计算公式如式(6)所示：
Figure BDA0003564768620000038
其中，hx(x,y)和hy(x,y)分别是像素点沿着x方向和y方向的梯度；
S4012：进行梯度伪排序：将梯度值从大到小进行1024等分，在每个分段中取出一个种子像素用来进行排序；
计算出由像素点构成的候选矩形区域的方向，通过比较候选像素点和候选矩形区域的方向判断该像素点是否属于矩形区域中，完成矩形区域划分；
S4013：进行矩形区域的近似计算：计算线段对应的近似矩形区域，所述矩形区域的尺寸覆盖整个线支撑区域，每个矩形区域代表一条候选线段；
矩形区域的中心点cx计算公式如式(7)所示：
Figure BDA0003564768620000041
其中(cx,cy)表示中心点的坐标，Region表示当前区域，G(j)表示区域内的像素点j的梯度值，下标j用于遍历区域内所有的像素点，x(j)和y(j)分别表示像素点j的横坐标和纵坐标；
S4014：对矩形区域进行改进，筛选出正确的线支撑区域：计算矩形区域内的误报数NFA值，若计算结果小于预设阈值，则认为当前矩形区域为正确的线支撑区域，所述NFA值的计算公式如式(8)所示：
NFA(r,t)＝(PQ)2.5γB(f(r),k(r,t),p) (8)
其中，P和Q为图像分辨率大小，B()表示二项分布，f(r)为r矩形区域内像素点的个数，k(r,t)表示r矩形区域中类内点的个数，p为梯度角对齐的概率，γ为p的不同值的可能的数量。
本发明第二方面提供了一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现本发明第一方面提供的无人机遥感多光谱图像拼接方法。
本发明第三方面提供了一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现本发明第一方面提供的无人机遥感多光谱图像拼接方法。
本发明的有益效果如下：
本申请所提供的方案，首先对获取的各个波段灰度图像进行波段配准校正来保证波段之间的一致性，再通过对融合后的单波段图像进行配准获得局部单应性矩阵集合来提升算法的效率，最后对于图片中线段结构特点，对拼接后的图像进行线段检测和校正，减少拼接可能产生的扭曲和变形现象，提高拼接结果的准确性和精度。
附图说明
下面结合附图对本发明的具体实施方式作进一步详细的说明。
图1示出本发明的一个实施例提供的一种无人机遥感多光谱图像拼接方法的流程图。
图2为本发明的另一个实施例提供的一种无人机遥感多光谱图像拼接方法的示意图。
图3为本发明的一个实施例提供的波段配准前的各个波段的图像。
图4为本发明的一个实施例提供的波段配准后的各个波段的图像。
图5示出本发明的一个实施例提供的参考图像划分的示意图。
图6示出实现本发明实施例提供的方法的计算机设备的结构示意图。
具体实施方式
为了更清楚地说明本发明，下面结合实施例和附图对本发明做进一步的说明。附图中相似的部件以相同的附图标记进行表示。本领域技术人员应当理解，下面所具体描述的内容是说明性的而非限制性的，不应以此限制本发明的保护范围。
发明人基于上述技术问题，提出了一种无人机遥感多光谱图像拼接方法、设备和介质。
如图1所示，本发明的一个实施例提出了一种无人机遥感多光谱图像拼接方法，包括：
S10：获取N帧单波段图像并进行波段配准校正以得到N帧多光谱图像，其中，每帧单波段图像包含M个波段的源图像，其中，
N和M为正整数；
S20：分别对所述N帧多光谱图像进行波段融合得到N帧灰度图像，基于所述N帧灰度图像进行拼接以获得N个局部单应性矩阵集合；
S30：通过第n个局部单应性矩阵集合Dn分别对第n张单波段图像的第m波段的灰度图像进行图像配准，将配准后的相同波段图像进行拼接以获得M个单波段拼接结果，对所述M个单波段拼接结果进行融合处理，得到初步拼接图像，其中，
n＝1、2、……、N；m＝1、2、……、M；
S40：对所述初步拼接图像中存在的线段结构进行检验与矫正得到最终拼接结果。
本方案首先对获取的各个波段灰度图像进行波段配准校正来保证波段之间的一致性，再通过对融合后的单波段图像进行配准获得局部单应性矩阵集合来提升算法的效率，最后对于农业场景中特有的农田、道路线段结构特点，对拼接后的图像进行线段检测和校正，减少拼接可能产生的扭曲和变形现象，提高拼接结果的准确性和精度。
无人机多光谱照相机是由几台单波段相机组合而成的，各台相机分别具有独立的光路，各自记录对应成像波段的数据。由于镜头间的光路是相互独立的，获取到的各个波段的灰度图像之间存在细微的差异，这些差异虽然比较小，但是也会影响到后续的拼接，因此需要针对同一场景拍摄的各个波段的灰度图像进行波段配准。
如图2所示，以具有6台单波段光谱成像仪的无人机为例，首先，无人机按照预设路线采集N帧视频，如图3所示，为波段配准前各个波段图像，对于各个波段的单波段灰度图像数据源来说，同一场景拍摄的不同波段的小图之间会存在一些差异，因此需要对这些图像进行波段配准校正，如图4所示，为波段配准后的各个波段图像，校正后的图像保证了波段之间的空间一致性，以便于后续的图像拼接，考虑到SIFT算法具有尺度不变性，并且能够快速准确地完成灰度图像配准，因此采用SIFT算法进行波段配准校正，
具体的，
进行SIFT特征提取得到图像的特征描述子，对于提取出的特征描述子进行欧式距离算法匹配，再利用匹配点计算变换矩阵，对待配准图像进行图像变换与插值，得到波段配准后的图像。
在一种可能的实现方式中，所述S10包括：
S101：通过无人机多光谱照相机沿特定路线获取视频，所述视频包括N帧单波段图像，每帧单波段图像包含M个波段的源图像；
S102：利用SIFT算法分别对所述N帧单波段图像进行波段配准校正：
通过不同尺寸的高斯核函数与当前单波段图像进行卷积，在不同的层级中对卷积之后的结果进行图像缩小以得到尺度空间，并在尺度空间中对源矩阵和网矩阵进行极值检测以得到特征点，同时移除特征点集中的边缘响应点，为关键点分配方向值；
根据所述方向值生成特征点描述子，对特征点进行匹配得到匹配点；
利用匹配点计算变换矩阵对当前单波段图像的M个波段的源图像进行图像变换和差值，得到波段配准后的多光谱图像并对M个波段的源图像进行更新。
应该理解的是，利用SIFT算法进行特征提取时，为了保证尺度恒定性，首先要创造尺度(DOG)空间，并在不同层级上提取兴趣点。为了使得尺度恒定，SIFT算法通过不同尺寸的高斯核函数与图像进行卷积，在不同的层级中对卷积之后的结果进行图像缩小生成DOG空间，
其中，二维高斯核函数G(x,y,σ)的表达式如式(1-1)所示：
Figure BDA0003564768620000071
其中，σ是二维高斯核函数的标准差，也叫尺度因子。I(x,y)是原图像，尺度空间变换公式如式(1-2)所示：
L(x,y,σ)＝I(x,y)*G(x,y,σ) (1-2)
SIFT算法会通过创建高斯金字塔来获得图像的缩略图。在下采样的基础上使用不同模型对每一级图像使用高斯平滑处理，就可以得到高斯金字塔。在高斯金字塔中，图像的组数(Octave)和层数(Interval)一样多，多层图像构成了一组。当前图像组的第一层通过前一组图像的的倒数第三张图被下采样2次获得。
SIFT算法使用构造DOG尺度空间的方式来定位鲁棒性的特征点。DOG尺度空间是通过对同一层中紧挨的两幅图像相减得到的。那么计算该尺度空间的D(x,y,σ)函数表达式如式(1-3)所示：
D(x,y,σ)＝(G(x,y,kσ)-G(x,y,σ))*I(x,y) (1-3)
然后，在DOG尺度空间中进行检测来寻找特征点：确定特征点在DOG尺度空间中位置的具体方法是在每组的每一层图像上，将每一个点和它周围相邻所有的点进行比较，包括同一尺度平面上的点以及上下平面的点。如果该像素点的值比其余的点都大或者都小，则认为这个像素点就是一个特征点。在构建的尺度空间中搜索寻找出了很多特征点，但是其中包括了一些具有边缘效应的特征点需要被筛选掉。由于特征点弯曲程度在边界处梯度方向比较大，而与之正交的方向较小。选择出的特征点的弯曲程度与2×2的Hessian矩阵特征值成正比，H矩阵如式(1-4)所示：
Figure BDA0003564768620000072
其中Dxx和Dyy分别表示沿着x和y方向的二阶偏导数，Dxy表示先沿x再沿y方向的二阶偏导数。设α和β分别是矩阵H的特征值，且α>β，则矩阵H的迹和特征值如式(1-5)和式(1-6)所示：
Tr(H)＝Dxx+Dyy＝α+β (1-5)
Det(H)＝DxxDyy-(Dxy)2＝αβ (1-6)
主曲率的变化是通过计算迹和特征值的比率来表示的。若要移除边缘响应点的话，就需要判断该点在水平和垂直边缘方向上的差别是否很大，如果有差异很小则说明是这并不是需要移除的点，反之则需要移除。因此r不能超过一定的阈值，只要控制r小于该阈值，那么迹和特征值的比率也将会小于指定的阈值，这样就使两个方向上的主曲率的差异在一定范围内，也就删除了边缘响应点。
最后，生成特征点描述子：特征描述子包含了该特征点所在的尺度、位置以及方向这三个维度的信息。因此在生成特征点描述子时，首先计算的就是关键点和它周围的点的梯度的大小和偏向角度；再利用这些信息画出梯度统计直方图，直方图中最高的对应了关键特征点的偏向角度。计算特征点的方向时，将以特征点周围区域划分成4×4的网格，那么在这16个小格中，每个小格内会包含8个方向的信息，因此总共有128维向量，也就是说特征描述子有128维。然后利用正确的匹配点计算变换矩阵，对待配准图像进行图像变换与插值，得到波段配准后的图像。
在一个具体的实施例中，所述特征描述子包括特征点所在的尺度、位置以及方向。
在一个具体的实施例中，
所述S20包括：
S201：分别对所述N张多光谱图像进行波段融合得到N张灰度图像；
S202：分别对所述N张灰度图像进行单波段配准；
S203：采用逼近投影变换算法进行所述N张灰度图像的拼接以得出对应的N个局部变换矩阵集合。
在一个具体的实施例中，将波段配准后的图像进行波段融合，对融合后的灰度图像进行灰度图像配准以获得适用于全波段的局部单应性矩阵集合，再将其应用于各个波段的图像进行图像配准与融合，得到拼接后的图像，
在具体实现过程中，首先利用SIFT算法特征提取和匹配找到正确的匹配点对，接下来进行分块计算移动直线线性变换权重，将源图像均匀地划分成一个个小的网格，利用加权方式为图像中每个小的网格计算局部的投影变换矩阵进行映射，来和目标图像对齐。得到拼接结果后，再对图像中的农田和道路等线段结构进行检测与校正，得到更为准确的多光谱图像拼接结果。
在一种可能的实现方式中，
所述S203包括：
S2031：分别对所述N张灰度图像进行特征提取和匹配得到预变换矩阵；
S2032：进行分块计算移动直线线性变换权重，对灰度图像进行分块将图像均匀地划分成一个个小的网格，划分后的图像如图5所示，利用加权方式为图像中每个小的网格计算局部的投影变换矩阵进行映射；移动直线线性变换算法通过计算当前网格内的匹配点与网格中心的距离加权平方和误差来得到所述局部单应性变换矩阵集合Dn，其中，
Figure BDA0003564768620000091
K为当前灰度图像中的网格总数；
S2033：利用所述局部单应性变换矩阵集合对所述灰度图像进行变换，将变换后的图像与参考图像利用加权平均融合法进行图像融合与拼接，
其中，
第n张灰度图像的第k个网格的局部单应性矩阵
Figure BDA0003564768620000092
的计算公式如式(1)所示，通过求解最小特征值对应的特征向量来估计局部变换矩阵
Figure BDA0003564768620000093
的值：
Figure BDA0003564768620000094
其中，ai为第n张灰度图像的第k个网格内匹配点之间的线性映射函数，l为匹配点对数；d为匹配点与网格中心的距离，
Figure BDA0003564768620000095
为特征点所对应的权重，Wk为权重对应的对角矩阵，A为映射函数对应的映射矩阵形式，其中，
Figure BDA0003564768620000096
Figure BDA0003564768620000097
其中，xk和yk是第k个网格的中心点坐标，
Figure BDA0003564768620000098
和
Figure BDA0003564768620000099
是第k个网格内第i个匹配点的坐标，σ为高斯尺度因子，其范围在(0,1)之间。
在一个具体的实施例，如图2所示，通过第n个局部单应性矩阵集合Dn分别对第n张单波段图像的第m波段的灰度图像进行图像配准，将配准后的相同波段图像进行拼接以获得M个单波段拼接图像，对所述M个单波段拼接图像进行融合处理，得到初步拼接图像，具体的，
通过对应单波段图像1的第1个局部单应性矩阵集合D1对第1张单波段图像的波段1～6的灰度图像进行波段配准，通过对应单波段图像2的第2个局部单应性矩阵集合D2对第2张单波段图像的波段1～6的灰度图像进行波段配准；依此类推；
将配准后的相同波段的灰度图像进行拼接，得到6个单波段拼接结果，对其进行图像融合，得到初步拼接图像。
在一个具体的实施例中，所述S40包括：
S401：进行线段结构的检测与提取，提取所述初步拼接图像中的线支撑区域，并通过所述线支撑区域得到线段特征，其中，
所述线支撑区域由有线段存在的像素区域组成；
S402：进行线段校正：得到线段特征后，利用最小二乘法在矩形区域拟合出一条线段，将所述线段的交点作为特征点，将其与矩形区域中心点的位置关系进行排序，根据变换后的特征点与原来的特征点之间的对应关系计算透视变换矩阵，将待校正的区域内所有的点进行映射，完成线段校正。
在一个具体的实施例中，所述S401包括：
S4011：利用每个像素点右下方的四个像素来计算该像素点对应的梯度以得到线段边缘的位置，像素点沿着x方向和y方向的梯度hx(x,y)和hy(x,y)计算公式如式(4)和式(5)所示：
Figure BDA0003564768620000101
Figure BDA0003564768620000102
其中h(x，y)表示坐标为(x，y)的像素点对应的灰度值；
梯度大小H(x,y)和方向θ(x,y)的计算公式如式(6)所示：
Figure BDA0003564768620000103
其中，hx(x,y)和hy(x,y)分别是像素点沿着x方向和y方向的梯度；
S4012：进行梯度伪排序：将梯度值从大到小进行1024等分，在每个分段中取出一个种子像素用来进行排序；
计算出由像素点构成的候选矩形区域的方向，通过比较候选像素点和候选矩形区域的方向判断该像素点是否属于矩形区域中，完成矩形区域划分；
S4013：进行矩形区域的近似计算：计算线段对应的近似矩形区域，所述矩形区域的尺寸覆盖整个线支撑区域，每个矩形区域代表一条候选线段；
矩形区域的中心点cx计算公式如式(7)所示：
Figure BDA0003564768620000111
其中(cx,cy)表示中心点的坐标，Region表示当前区域，G(j)表示区域内的像素点j的梯度值，下标j用于遍历区域内所有的像素点，x(j)和y(j)分别表示像素点j的横坐标和纵坐标；
S4014：对矩形区域进行改进，筛选出正确的线支撑区域：计算矩形区域内的误报数NFA值，若计算结果小于预设阈值，所述预设阈值例如为1，则认为当前矩形区域为正确的线支撑区域，所述矩形区域r内的误报数NFA值的计算公式如式(8)所示：
NFA(r,t)＝(PQ)2.5γB(f(r),k(r,t),p) (8)
其中，P和Q为图像分辨率大小，B()表示二项分布，f(r)为矩形区域r内像素点的个数，k(r,t)表示矩形区域r中类内点的个数，p为梯度角对齐的概率，γ为p的不同值的可能的数量。
应当说明的是，梯度角对齐指的是矩形区域中所有水平线方向角度与矩形方向角度偏差小于阈值，就认为是对齐的，所述阈值例如为5°、6°或4°。
如图6所示，本发明的第二个实施例提供的一种计算机设备的结构示意图。适于用来实现上述实施例提供的无人机遥感多光谱图像拼接方法的计算机设备，包括中央处理模块(CPU)，其可以根据存储在只读存储器(ROM)中的程序或者从存储部分加载到随机访问存储器(RAM)中的程序而执行各种适当的动作和处理。在RAM中，还存储有计算机设备操作所需的各种程序和数据。CPU、ROM以及RAM通过总线被此相连。输入/输入(I/O)接口也连接至总线。
以下部件连接至I/O接口:包括键盘、鼠标等的输入部分；包括诸如液晶显示器(LCD)等以及扬声器等的输出部分；包括硬盘等的存储部分；以及包括诸如LAN卡、调制解调器等的网络接口卡的通信部分。通信部分经由诸如因特网的网络执行通信处理。驱动器也根据需要连接至I/O接口。可拆卸介质，诸如磁盘、光盘、磁光盘、半导体存储器等等，根据需要安装在驱动器上，以便于从其上读出的计算机程序根据需要被安装入存储部分。
特别地，根据本实施例，上文流程图描述的过程可以被实现为计算机软件程序。例如，本实施例包括一种计算机程序产品，其包括有形地包含在计算机可读介质上的计算机程序，上述计算机程序包含用于执行流程图所示的方法的程序代码。在这样的实施例中，该计算机程序可以通过通信部分从网络上被下载和安装，和/或从可拆卸介质被安装。
附图中的流程图和示意图，图示了本实施例的系统、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或示意图中的每个方框可以代表一个模块、程序段或代码的一部分，上述模块、程序段或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意，在有些作为替换的实现中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个接连地表示的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，示意图和/或流程图中的每个方框、以及示意和/或流程图中的方框的组合，可以用执行规定的功能或操作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。
本申请的第四个实施例提供了一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现：
S10：获取N帧单波段图像并进行波段配准校正以得到N帧多光谱图像，其中，每帧单波段图像包含M个波段的源图像，其中，
N和M为正整数；
S20：分别对所述N张多光谱图像进行波段融合得到N张灰度图像，基于所述N张灰度图像进行拼接以获得N个局部单应性矩阵集合；
S30：通过第n个局部单应性矩阵集合Dn分别对第n张单波段图像的第m波段的灰度图像进行图像配准，将配准后的相同波段图像进行拼接以获得M个单波段拼接图像，对所述M个单波段拼接图像进行融合处理，得到初步拼接图像，其中，
n＝1、2、……、N；m＝1、2、……、M；
S40：对所述初步拼接图像中存在的线段结构进行检验与矫正得到最终拼接结果。
在实际应用中，所述计算机可读存储介质可以采用一个或多个计算机可读的介质的任意组合。计算机可读介质可以是计算机可读信号介质或者计算机可读存储介质。计算机可读存储介质例如可以是但不限于电、磁、光、电磁、红外线、或半导体的系统、装置或器件，或者任意以上的组合。计算机可读存储介质的更具体的例子(非穷举的列表)包括：具有一个或多个导线的电连接、便携式计算机磁盘、硬盘、随机存取存储器(RAM)、只读存储器(ROM)、可擦式可编程只读存储器(EPROM或闪存)、光纤、便携式紧凑磁盘只读存储器(CD-ROM)、光存储器件、磁存储器件、或者上述的任意合适的组合。在本实施例中，计算机可读存储介质可以是任何包含或存储程序的有形介质，该程序可以被指令执行系统、装置或者器件使用或者与其结合使用。
计算机可读的信号介质可以包括在基带中或者作为载波一部分传播的数据信号，其中承载了计算机可读的程序代码。这种传播的数据信号可以采用多种形式，包括但不限于电磁信号、光信号或上述的任意合适的组合。计算机可读的信号介质还可以是计算机可读存储介质以外的任何计算机可读介质，该计算机可读介质可以发送、传播或者传输用于由指令执行系统、装置或者器件使用或者与其结合使用的程序。
计算机可读介质上包含的程序代码可以用任何适当的介质传输，包括但不限于无线、电线、光缆、RF等等，或者上述的任意合适的组合。
可以以一种或多种程序设计语言或其组合来编写用于执行本申请操作的计算机程序代码，所述程序设计语言包括面向对象的程序设计语言—诸如Java、Smalltalk、C++，还包括常规的过程式程序设计语言—诸如“C”语言或类似的程序设计语言。程序代码可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或服务器上执行。在涉及远程计算机的情形中，远程计算机可以通过任意种类的网络——包括局域网(LAN)或广域网(WAN)—连接到用户计算机，或者，可以连接到外部计算机(例如利用因特网服务提供商来通过因特网连接)。
显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定，对于本领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。

Claims (10)

1.一种无人机遥感多光谱图像拼接方法，其特征在于，包括：
S10：获取N帧单波段图像并进行波段配准校正以得到N帧多光谱图像，其中，每帧单波段图像包含M个波段的源图像，其中，
N和M为正整数；
S20：分别对所述N张多光谱图像进行波段融合得到N张灰度图像，基于所述N张灰度图像进行拼接以获得N个局部单应性矩阵集合；
S30：通过第n个局部单应性矩阵集合Dn分别对第n张单波段图像的第m波段的灰度图像进行图像配准，将配准后的相同波段图像进行拼接以获得M个单波段拼接结果，对所述M个单波段拼接结果进行融合处理，得到初步拼接图像，其中，
n＝1、2、……、N；m＝1、2、……、M；
S40：对所述初步拼接图像中存在的线段结构进行检验与矫正得到最终拼接结果。
2.根据权利要求1所述的方法，其特征在于，
所述S10包括：
S101：通过无人机多光谱照相机沿特定路线获取视频，所述视频包括N帧单波段图像，每帧单波段图像包含M个波段的源图像；
S102：利用SIFT算法分别对所述N帧单波段图像进行波段配准校正：通过不同尺寸的高斯核函数与当前单波段图像进行卷积，在不同的层级中对卷积之后的结果进行图像缩小以得到尺度空间，并在尺度空间中对源矩阵和网矩阵进行极值检测以得到特征点，同时移除特征点集中的边缘响应点，为关键点分配方向值；
根据所述方向值生成特征点描述子，对特征点进行匹配得到匹配点；
利用匹配点计算变换矩阵对当前单波段图像的M个波段的源图像进行图像变换和差值，得到波段配准后的多光谱图像并对M个波段的源图像进行更新。
3.根据权利要求2所述的方法，其特征在于，
所述特征描述子包括特征点所在的尺度、位置以及方向。
4.根据权利要求1所述的方法，其特征在于，
所述S20包括：
S201：分别对所述N张多光谱图像进行波段融合得到N张灰度图像；
S202：分别对所述N张灰度图像进行单波段配准；
S203：采用逼近投影变换算法进行所述N张灰度图像的拼接以得出对应的N个局部变换矩阵集合。
5.根据权利要求4所述的方法，其特征在于，
所述S203包括：
S2031：分别对所述N张灰度图像进行特征提取和匹配得到预变换矩阵；
S2032：对灰度图像进行分块将图像均匀地划分成一个个小的网格，移动直线线性变换算法通过计算当前网格内的匹配点与网格中心的距离加权平方和误差来得到所述局部单应性变换矩阵集合Dn，其中，
Figure FDA0003564768610000021
K为当前灰度图像中的网格总数；
S2033：利用所述局部单应性变换矩阵集合对所述灰度图像进行变换，将变换后的图像与参考图像利用加权平均融合法进行图像融合与拼接。
6.根据权利要求5所述的方法，其特征在于，
第n张灰度图像的第k个网格的局部单应性矩阵
Figure FDA0003564768610000022
的计算公式如式(1)所示，通过求解最小特征值对应的特征向量来估计局部变换矩阵
Figure FDA0003564768610000023
的值：
Figure FDA0003564768610000024
其中，ai为匹配点之间的线性映射函数，l为匹配点对数；d为匹配点与网格中心的距离，
Figure FDA0003564768610000025
为特征点所对应的权重，Wk为权重对应的对角矩阵，A为映射函数对应的映射矩阵形式，其中，
Figure FDA0003564768610000026
Figure FDA0003564768610000027
其中，xk和yk是第k个网格的中心点坐标，
Figure FDA0003564768610000028
和
Figure FDA0003564768610000029
是第k个网格内第i个匹配点的坐标，σ为高斯尺度因子，其范围在(0,1)之间。
7.根据权利要求1所述的方法，其特征在于，
所述S40包括：
S401：进行线段结构的检测与提取，提取所述初步拼接图像中的线支撑区域，并通过所述线支撑区域得到线段特征，其中，
所述线支撑区域由有线段存在的像素区域组成；
S402：进行线段校正：利用最小二乘法在矩形区域拟合出一条线段，将所述线段的交点作为特征点，将其与矩形区域中心点的位置关系进行排序，根据变换后的特征点与原来的特征点之间的对应关系计算透视变换矩阵，将待校正的区域内所有的点进行映射，完成线段校正。
8.根据权利要求7所述的方法，其特征在于，
所述S401包括：
S4011：利用每个像素点右下方的四个像素来计算该像素点对应的梯度以得到线段边缘的位置，像素点沿着x方向和y方向的梯度hx(x,y)和hy(x,y)计算公式如式(4)和式(5)所示：
Figure FDA0003564768610000031
Figure FDA0003564768610000032
其中h(x，y)表示坐标为(x，y)的像素点对应的灰度值；
梯度大小H(x,y)和方向θ(x,y)的计算公式如式(6)所示：
Figure FDA0003564768610000033
其中，hx(x,y)和hy(x,y)分别是像素点沿着x方向和y方向的梯度；
S4012：进行梯度伪排序：将梯度值从大到小进行1024等分，在每个分段中取出一个种子像素用来进行排序；
计算出由像素点构成的候选矩形区域的方向，通过比较候选像素点和候选矩形区域的方向判断该像素点是否属于矩形区域中，完成矩形区域划分；
S4013：进行矩形区域的近似计算：计算线段对应的近似矩形区域，所述矩形区域的尺寸覆盖整个线支撑区域，每个矩形区域代表一条候选线段；
矩形区域的中心点cx计算公式如式(7)所示：
Figure FDA0003564768610000034
其中(cx,cy)表示中心点的坐标，Region表示当前区域，G(j)表示区域内的像素点j的梯度值，下标j用于遍历区域内所有的像素点，x(j)和y(j)分别表示像素点j的横坐标和纵坐标；
S4014：对矩形区域进行改进，筛选出正确的线支撑区域：计算矩形区域内的误报数NFA值，若计算结果小于预设阈值，则认为当前矩形区域为正确的线支撑区域，所述NFA值的计算公式如式(8)所示：
NFA(r,t)＝(PQ)2.5γB(f(r),k(r,t),p) (8)
其中，P和Q为图像分辨率大小，B()表示二项分布，f(r)为r矩形区域内像素点的个数，k(r,t)表示r矩形区域中类内点的个数，p为梯度角对齐的概率，γ为p的不同值的可能的数量。
9.一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现如权利要求1-8任一项所述的方法。
10.一种计算机存储介质，其上存储有计算机程序，其特征在于，该程序被处理器执行时实现如权利要求1-8中任一项所述的方法。
CN202210299506.1A 2022-03-25 2022-03-25 一种无人机遥感多光谱图像拼接方法、设备及介质 Pending CN114897676A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210299506.1A CN114897676A (zh) 2022-03-25 2022-03-25 一种无人机遥感多光谱图像拼接方法、设备及介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210299506.1A CN114897676A (zh) 2022-03-25 2022-03-25 一种无人机遥感多光谱图像拼接方法、设备及介质

Publications (1)

Publication Number Publication Date
CN114897676A true CN114897676A (zh) 2022-08-12

Family

ID=82714647

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210299506.1A Pending CN114897676A (zh) 2022-03-25 2022-03-25 一种无人机遥感多光谱图像拼接方法、设备及介质

Country Status (1)

Country Link
CN (1) CN114897676A (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116402693A (zh) * 2023-06-08 2023-07-07 青岛瑞源工程集团有限公司 一种基于遥感技术的市政工程图像处理方法及装置
CN117036666A (zh) * 2023-06-14 2023-11-10 北京自动化控制设备研究所 基于帧间图像拼接的无人机低空定位方法

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116402693A (zh) * 2023-06-08 2023-07-07 青岛瑞源工程集团有限公司 一种基于遥感技术的市政工程图像处理方法及装置
CN116402693B (zh) * 2023-06-08 2023-08-15 青岛瑞源工程集团有限公司 一种基于遥感技术的市政工程图像处理方法及装置
CN117036666A (zh) * 2023-06-14 2023-11-10 北京自动化控制设备研究所 基于帧间图像拼接的无人机低空定位方法
CN117036666B (zh) * 2023-06-14 2024-05-07 北京自动化控制设备研究所 基于帧间图像拼接的无人机低空定位方法

Similar Documents

Publication Publication Date Title
Feng et al. Robust registration for remote sensing images by combining and localizing feature-and area-based methods
CN110211043B (zh) 一种用于全景图像拼接的基于网格优化的配准方法
CN106023086B (zh) 一种基于orb特征匹配的航拍影像及地理数据拼接方法
CN110223330B (zh) 一种可见光和红外图像的配准方法及系统
US8755624B2 (en) Image registration device and method thereof
US10311595B2 (en) Image processing device and its control method, imaging apparatus, and storage medium
CN107578376B (zh) 基于特征点聚类四叉划分和局部变换矩阵的图像拼接方法
CN111507901B (zh) 基于航带gps及尺度不变约束的航拍图像拼接定位方法
CN106447601B (zh) 一种基于投影-相似变换的无人机遥感影像拼接方法
CN114897676A (zh) 一种无人机遥感多光谱图像拼接方法、设备及介质
CN105608667A (zh) 一种全景拼接的方法及装置
CN111383252B (zh) 多相机目标追踪方法、系统、装置及存储介质
CN111340701A (zh) 一种基于聚类法筛选匹配点的电路板图像拼接方法
CN113326763B (zh) 一种基于边界框一致性的遥感目标检测方法
CN110084743B (zh) 基于多航带起始航迹约束的图像拼接与定位方法
CN113160048A (zh) 一种缝合线引导的图像拼接方法
US8126275B2 (en) Interest point detection
Lu et al. Multiperspective image stitching and regularization via hybrid structure warping
CN113298187A (zh) 图像处理方法及装置、计算机可读存储介质
Rui et al. Research on fast natural aerial image mosaic
CN111062341B (zh) 视频图像区域的分类方法、装置、设备及存储介质
Zarei et al. MegaStitch: Robust Large-scale image stitching
US20070280555A1 (en) Image registration based on concentric image partitions
Qi et al. Image stitching based on improved SURF algorithm
Feng et al. Registration of multitemporal GF-1 remote sensing images with weighting perspective transformation model

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination