CN114936971A - 一种面向水域的无人机遥感多光谱图像拼接方法及*** - Google Patents
一种面向水域的无人机遥感多光谱图像拼接方法及*** Download PDFInfo
- Publication number
- CN114936971A CN114936971A CN202210640221.XA CN202210640221A CN114936971A CN 114936971 A CN114936971 A CN 114936971A CN 202210640221 A CN202210640221 A CN 202210640221A CN 114936971 A CN114936971 A CN 114936971A
- Authority
- CN
- China
- Prior art keywords
- image
- point
- images
- pixel
- matrix
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 77
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 title claims abstract description 27
- 238000012937 correction Methods 0.000 claims abstract description 53
- 230000004927 fusion Effects 0.000 claims abstract description 19
- 238000012216 screening Methods 0.000 claims abstract description 15
- 238000011156 evaluation Methods 0.000 claims abstract description 11
- 238000007781 pre-processing Methods 0.000 claims abstract description 6
- 239000011159 matrix material Substances 0.000 claims description 87
- 238000004364 calculation method Methods 0.000 claims description 45
- 230000009466 transformation Effects 0.000 claims description 16
- 239000013598 vector Substances 0.000 claims description 15
- 230000008859 change Effects 0.000 claims description 14
- 238000000513 principal component analysis Methods 0.000 claims description 13
- 238000013519 translation Methods 0.000 claims description 13
- 230000008569 process Effects 0.000 claims description 11
- PXFBZOLANLWPMH-UHFFFAOYSA-N 16-Epiaffinine Natural products C1C(C2=CC=CC=C2N2)=C2C(=O)CC2C(=CC)CN(C)C1C2CO PXFBZOLANLWPMH-UHFFFAOYSA-N 0.000 claims description 9
- 238000003384 imaging method Methods 0.000 claims description 9
- 230000002194 synthesizing effect Effects 0.000 claims description 8
- 230000002159 abnormal effect Effects 0.000 claims description 7
- 238000010276 construction Methods 0.000 claims description 7
- 230000006870 function Effects 0.000 claims description 7
- 238000003860 storage Methods 0.000 claims description 7
- PCTMTFRHKVHKIS-BMFZQQSSSA-N (1s,3r,4e,6e,8e,10e,12e,14e,16e,18s,19r,20r,21s,25r,27r,30r,31r,33s,35r,37s,38r)-3-[(2r,3s,4s,5s,6r)-4-amino-3,5-dihydroxy-6-methyloxan-2-yl]oxy-19,25,27,30,31,33,35,37-octahydroxy-18,20,21-trimethyl-23-oxo-22,39-dioxabicyclo[33.3.1]nonatriaconta-4,6,8,10 Chemical compound C1C=C2C[C@@H](OS(O)(=O)=O)CC[C@]2(C)[C@@H]2[C@@H]1[C@@H]1CC[C@H]([C@H](C)CCCC(C)C)[C@@]1(C)CC2.O[C@H]1[C@@H](N)[C@H](O)[C@@H](C)O[C@H]1O[C@H]1/C=C/C=C/C=C/C=C/C=C/C=C/C=C/[C@H](C)[C@@H](O)[C@@H](C)[C@H](C)OC(=O)C[C@H](O)C[C@H](O)CC[C@@H](O)[C@H](O)C[C@H](O)C[C@](O)(C[C@H](O)[C@H]2C(O)=O)O[C@H]2C1 PCTMTFRHKVHKIS-BMFZQQSSSA-N 0.000 claims description 6
- 238000004422 calculation algorithm Methods 0.000 claims description 6
- 238000013178 mathematical model Methods 0.000 claims description 5
- 238000005070 sampling Methods 0.000 claims description 5
- 230000000007 visual effect Effects 0.000 claims description 5
- 230000003287 optical effect Effects 0.000 claims description 4
- 240000004282 Grewia occidentalis Species 0.000 claims description 3
- 238000012952 Resampling Methods 0.000 claims description 3
- 238000009826 distribution Methods 0.000 claims description 3
- 238000012887 quadratic function Methods 0.000 claims description 3
- 230000004044 response Effects 0.000 claims description 3
- 238000012545 processing Methods 0.000 abstract description 8
- 238000010586 diagram Methods 0.000 description 15
- 238000004590 computer program Methods 0.000 description 8
- 238000009360 aquaculture Methods 0.000 description 4
- 244000144974 aquaculture Species 0.000 description 4
- 238000009825 accumulation Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 230000008030 elimination Effects 0.000 description 2
- 238000003379 elimination reaction Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012544 monitoring process Methods 0.000 description 2
- 238000010606 normalization Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 241000282320 Panthera leo Species 0.000 description 1
- 206010034719 Personality change Diseases 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000000452 restraining effect Effects 0.000 description 1
- 238000006748 scratching Methods 0.000 description 1
- 230000002393 scratching effect Effects 0.000 description 1
- NTHWMYGWWRZVTN-UHFFFAOYSA-N sodium silicate Chemical compound [Na+].[Na+].[O-][Si]([O-])=O NTHWMYGWWRZVTN-UHFFFAOYSA-N 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T3/00—Geometric image transformations in the plane of the image
- G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
- G06T3/4038—Image mosaicing, e.g. composing plane images from plane sub-images
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T3/00—Geometric image transformations in the plane of the image
- G06T3/02—Affine transformations
-
- 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 using two or more images, e.g. averaging or subtraction
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/80—Geometric correction
-
- 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
-
- 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/761—Proximity, similarity or dissimilarity measures
-
- 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/77—Processing image or video features in feature spaces; using data integration or data reduction, e.g. principal component analysis [PCA] or independent component analysis [ICA] or self-organising maps [SOM]; Blind source separation
- G06V10/7715—Feature extraction, e.g. by transforming the feature space, e.g. multi-dimensional scaling [MDS]; Mappings, e.g. subspace methods
-
- 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/20016—Hierarchical, coarse-to-fine, multiscale or multiresolution image processing; Pyramid transform
-
- 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
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/30—Assessment of water resources
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Multimedia (AREA)
- Computing Systems (AREA)
- Artificial Intelligence (AREA)
- Health & Medical Sciences (AREA)
- Databases & Information Systems (AREA)
- Evolutionary Computation (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Software Systems (AREA)
- Image Analysis (AREA)
- Image Processing (AREA)
- Processing Or Creating Images (AREA)
Abstract
本发明公开了一种面向水域的无人机遥感多光谱图像拼接方法及***,属于水域图像处理技术领域。本发明是一种面向水域的无人机遥感多光谱图像拼接方法,根据经纬度POS信息进行图像的几何校正预处理,得到校正图像,并对校正图像添加地理坐标投影信息;然后根据地理坐标投影信息,进行图像的粗配准,得到配准图像;再根据配准图像,利用计算机图形学知识,计算出重叠区域位置,得到待拼接图像;进而构建基于主成分图像的尺度不变特征变换SIFT模型,提取待拼接图像的特征点;并以空间距离作为特征点的相似性评估指标,筛选出符合要求的匹配点;最后构建多分辨率融合模型,对待拼接图形进行融合,实现无人机遥感多光谱图像的拼接。
Description
技术领域
本发明涉及一种面向水域的无人机遥感多光谱图像拼接方法及***,属于水域图像处理技术领域。
背景技术
中国专利(公开号CN 110503604 A)公开了一种基于高精度POS的航空面阵影像实时正射拼接方法,包括飞行前,进行相机检校、相机与POS***的检校;获取每张原始影像的高精度POS数据,利用相机检校结果去除原始影像的图像畸变;将地表平均高程面作为物方平面;计算每张影像的投影矩阵,再通过地表平均高程,计算物方平面到原始像方平面的唯一确定的正定矩阵;确定统一分辨率,得到纠正影像平面到原始像方平面的变换矩阵;更新总纠正影像的范围,重新计算总纠正影像的宽和高;获得当前纠正影像在总纠正影像中的位置,进行影像拼接,对于重叠区部分直接使用影像DOM上已包含的部分;迭代执行直到得到最终的总纠正影像。本发明提高了正射纠正影像拼接的精度,满足影像实时拼接的需求。
上述方案,运用地理坐标信息进行图像拼接,但由于需要提供高精度的GPS信息,否则遥感图像获取的地理信息存在较大的误差,导致拼接图像之间位置错位,而无人机的稳定性和抗风能力以及GPS测量设备性能,使得每张图像对应的经纬度POS信息精度难以保证,进而导致拼接精度将受到严重影响。
发明内容
针对现有技术的缺陷,本发明的目的一在于提供一种根据经纬度POS信息进行图像的几何校正预处理,得到校正图像,并对校正图像添加地理坐标投影信息;然后根据地理坐标投影信息,进行图像的粗配准,得到配准图像;再根据配准图像,构建计算机图形学知识模型,计算出重叠区域位置,并抠取重叠区域图像,得到待拼接图像;进而构建基于主成分图像的尺度不变特征变换SIFT模型,提取待拼接图像的特征点;并以空间距离作为特征点的相似性评估指标,筛选出符合要求的匹配点;最后构建基于拉普拉斯Laplacian金字塔的多分辨率融合模型,对待拼接图形进行融合,实现无人机遥感多光谱图像的拼接,方案科学合理,切实可行的面向水域的无人机遥感多光谱图像拼接方法。
针对现有技术的缺陷,本发明的目的二在于提供一种通过处理器实现上述的一种面向水域的无人机遥感多光谱图像拼接方法,能够有效提高图像拼接精度,方案简单,实用,便于推广使用的计算机设备。
为实现上述目的之一,本发明的第一种技术方案为:
一种面向水域的无人机遥感多光谱图像拼接方法,
包括以下步骤:
第一步,接收图像数据,并合成多波段图像;
同时获取多波段图像对应的经纬度POS信息,并剔除经纬高数据异常值;
第二步,根据第一步中对应的经纬度POS信息进行图像的几何校正预处理,得到校正图像;
第三步,基于图像对应的经纬和偏航角信息,计算地理坐标投影信息,并对第二步中的校正图像添加地理坐标投影信息;
第四步,按照拼接策略选取无人机图像序列,根据第三步中的地理坐标投影信息,进行图像的粗配准,得到平移变换矩阵和两幅配准图像;
第五步,根据第四步中的配准图像,利用计算机图形学知识,计算出重叠区域位置,并抠取重叠区域图像,得到待拼接图像;
第六步,利用主成分分析法PCA对第五步中的待拼接图像提取第一主成分图像,并对第一主成分图像构建尺度不变特征模型SIFT,以获取第一主成分图像的特征点;
第七步,判断第六步中的特征点数量,若数量少于N个,将第四步中的平移变换矩阵作为后期的单应性变化矩阵;
否则采用查询索引KD树方式,建立一对多假设匹配集,并以空间距离作为第六步中特征点的相似性评估指标,筛选出符合要求的匹配点;并使用距离直方图约束全局相似属性来剔除误匹配点对,计算单应性变换矩阵;
第八步,对第四步中的两幅配准图像构建基于拉普拉斯Laplacian金字塔的多分辨率融合模型,并结合第七步中的单应性变换矩阵,对第四步中的两幅配准图像进行融合,实现无人机遥感多光谱图像的拼接;
之后返回第四步,直到完成所有序列图像的拼接后退出。
本发明经过不断探索以及试验,根据经纬度POS信息进行图像的几何校正预处理,得到校正图像,并对校正图像添加地理坐标投影信息;然后根据地理坐标投影信息,进行图像的粗配准,得到配准图像;再根据配准图像,构建计算机图形学知识模型,计算出重叠区域位置,并抠取重叠区域图像,得到待拼接图像;进而构建基于主成分图像的尺度不变特征变换SIFT模型,提取待拼接图像的特征点;并以空间距离作为特征点的相似性评估指标,筛选出符合要求的匹配点;最后构建基于拉普拉斯Laplacian金字塔的多分辨率融合模型,对待拼接图形进行融合,实现无人机遥感多光谱图像的拼接。
进而,本发明特征点融合POS信息的方案,能够有效融合了特征点方法配准精准度高的特点,以及位姿拼接法配准速度快的特点,拼接方案可作为实时配准的基础,并降低了对POS信息的精度要求;同时特征点结合经纬度POS信息可以降低大规模图像拼接中出现的误差累积,从而使得拼接图配准精度明显提高,最终得到的拼接融合图像不仅含有地理坐标信息,而且具有良好视觉效果,满足一般的无人机水质遥感监测时的图像拼接要求,方案科学合理,切实可行。
进一步,本发明适用于渔业养殖水域环境,在特征匹配算法中,针对养殖水域图像检测的特征点幅值偏小,方向不明确,以及在图像中存在许多具有局部相似属性的特征点,提供了更好的剔除方法,使得获得的单应性变化矩阵更准确。
更进一步,本发明能够处理长条状,规划航线等图像序列的拼接,且最终的拼接图具有地理坐标。
优选的,N=4。
作为优选技术措施:
所述第一步中,多波段图像的合成方法如下:
接收无人机遥感的多光谱图像,并将同时刻拍摄的同一位置的单波段图像进行合成,得到多波段图像;
并根据获取单波段图像对应的经纬度POS信息剔除起飞和返航时的经纬高数据异常值;
所述经纬度POS信息包括经度lon、纬度lat、航高H、航向角γ、俯仰角α、滚转角β、地面分辨率∈以及测量***的地理坐标系和投影坐标系;
地面分辨率∈的计算公式如下:
∈=pixelsize*H/f (1)。
作为优选技术措施:
所述第二步中,几何校正的方法包括以下内容:
步骤21,建立图像坐标系O-xy,相机坐标系S-XsYsZs,机体坐标系P-XpYpZp,地理坐标系E-XEYEZE;
步骤22,对步骤21中的图像坐标系选择以像主点为原点,沿飞行方向为t轴的正方向,垂直于飞行方向为x轴的正方向时,根据相机与机体的相对位置关系及无人机的飞行姿态,进行图像坐标系下像点g(x,y)与地理坐标系下对应物点G(XE,YE,ZE)的变化关系转换;
变化关系转换的计算公式如下:
[XE,YE,ZE]E T=λREPRPSRSO[x,y]O T=λREPRPS[x,y,-f]O T (2)
其中,f为焦距,用于表示图像坐标系与相机坐标系在Z轴方向的平移量;
步骤23,由于相机的镜头光心与无人机的质心重合,且相机坐标轴系经过平移后与机体坐标系完全重合,即RPS为单位矩阵I;
根据步骤22中的平移量,构建相机成像模型,其计算公式如下:
[XE,YE,ZE]E T=λREP[x,y,-f]O T (3)
其中,REP=R(H)R(γ)R(β)R(α),λ为比例系数,即λ=H/f;
在正直摄影条件下,遥感图像中像点g(x,y)与校正后的像点g′(x′,y′)的满足关系式如下:
(x′,y′,-f)O T=REP(x,y,-f)O T (4)
其中,R(α)、R(β)、R(γ)、R(H)分别为基于俯仰角的校正旋转矩阵、基于滚动角的校正旋转矩阵、基于偏航角角的校正旋转矩阵、基于航高的校正矩阵;
步骤24,对步骤23中的R(α)、R(β)、R(γ)、R(H)分别建立俯仰角校正矩阵、滚动角校正矩阵、偏航角校正矩阵、高度校正矩阵,得到校正数学模型;
其中,θ表示任一像素点和焦点的连线与中心视轴之间的夹角;
步骤25,根据步骤24中的校正数学模型,计算图像四角坐标校正后的新坐标,通过旧坐标与新坐标求出变换矩阵,并采用双线性插值重采样方法计算出新的校正图像,完成图像的几何校正,并更新图像的地理坐标信息。
作为优选技术措施:
所述第三步中,地理坐标投影信息的计算方法如下:
利用仿射矩阵参数对遥感图像坐标与地理坐标进行转换,其包括6个参数,分别为XE,Xpixel,Rγ,YE,Ypixel,Rγ,描述的是图像行列号和地理坐标之间的关系,
其中,XE、YE表示图像左上角像元的地理投影坐标,Xpixel、Ypixel分别表示图像像元在经度、维度方向的地面分辨率,Rγ表示图像旋转角度的正弦值;
将图像绕中心O点顺时针旋转γ度,然后根据无人机记录的经纬度经过坐标投影计算得到O点坐标(XOE,YOE),则图像左上角像点G坐标的计算公式如下:
其中w和h分别表示图像尺寸的宽和高;
地面分辨率Xpixel=-Ypixel=-∈,此时Rγ=0;
图像坐标系下任一点g(row,col)的地理坐标G(XE′,YE′)的计算公式如下:
作为优选技术措施:
所述第四步中,粗配准的过程如下:
在拼接的前几轮使用帧到帧拼接策略,后续选择拼接图到拼接图的拼接策略以完成图像的拼接;
获取同一航线上相邻的图像I1和图像I2,图像左上角的地理坐标分别为G1(XE1,YE1)和G2(XE2,YE2),以图像I1为基准,求得图像I2相对于图像I1的偏移量;
所述偏移量的计算公式如下:
根据偏移量,将几何校正后的图像I1和图像I2分别绕O1和O2旋转γ度至航线方向,图像I2中的任一点在图像I1中的图像坐标的计算公式如下:
其中,Hrigid表示平移变化矩阵;
配准图像的仿射矩阵参数的左上角像元地理坐标G(XE new,YE new)的计算公式如下:
配准图像的地面分辨率Xpixel、Ypixel和Rγ均保持不变。
作为优选技术措施:
所述第五步中,抠取重叠区域图像包括以下内容:
步骤51,将图像I1和I2分别绕O1和O2旋转γ度至航线方向;
步骤52,步骤51中的旋转完成后,根据计算机图形学知识,计算出图像I1四角构成的四边形I1AI1BI1CI1D与图像I2四角构成的四边形I2AI2BI2CI2D的多边形重叠区域;
步骤53,在步骤52中的多边形重叠区域向外添加一定偏置δ,最终获得多边形重叠区域ABCD的坐标;
步骤54,利用步骤53中的多边形重叠区域ABCD的坐标制作掩膜图像,从图像I1和I2中获得只含重叠区域的图像,记为图像I1′和I2′。
作为优选技术措施:
所述第六步中,构建第一主成分图像的尺度不变特征变换SIFT模型的方法如下:
步骤61,输入大小为w×h,波段数为c的多光谱图像矩阵I,将I重组为w×h行c列矩阵Ireshape,对Ireshape每一列进行归一化得到I′reshape;
I′reshape的计算公式如下:
步骤62,根据步骤61中的I′reshape,计算协方差矩阵Cov;
协方差矩阵Cov的计算公式如下:
步骤63,计算步骤62中的协方差矩阵Cov最大特征值对应的特征向量Vec;
步骤64,对重组图像的Ireshape提取步骤63中的特征向量Vec的主成分图像,并重组为w×h大小的图像IPCA;
图像IPCA的计算公式如下:
IPCA=Ireshape·Vec (20)
步骤65,将步骤64中的主成分图像作为尺度不变特征变换SIFT特征点提取的输入,然后进行构建DOG尺度空间、检测DOG尺度空间的极值点,删除不稳定的特征点、对特征点方向进行赋值和生成特征点描述子;
所述构建DOG尺度空间,包括以下内容:
通过对源图像进行高斯模糊和降采样得到图像金字塔,降采样的计算公式如下:
其次对图像金字塔的每层图像使用不同的σ参数进行高斯模糊,得到的多张模糊过图像构成了高斯金字塔,其计算公式如下:
然后将相邻的两个高斯空间的图像相减得到DOG图像,其计算公式如下:
式中,Down表示降采样,其中G0=I;I(x,y)表示源图像,L(x,y,σ)表示原图像卷积后的高斯尺度空间,代表卷积运算,σ表示高斯卷积核的尺度因子;G(x,y,σ)表示高斯核函数,(m,n)代表卷积核的大小;k表示相邻尺度空间的比例因子,取
检测DOG尺度空间的极值点,包括以下内容:
将每个像素点与同一尺度领域的8个像素点和上下相邻的尺度的18个像素点进行比较,只有当该像素点的DOG值全部大于或小于对比的26个像素点的DOG值时,该像素点是DOG尺度空间上一个极值点;
删除不稳定的特征点,包括以下内容:
首先,为了得到极值点的精确位置,需要对离散空间进行像素差值,然后通过类比拟合三维二次函数以得到极值点精确位置,最后剔除低对比度的极值点;
同时,对边缘响应大的极值点进行舍弃;
对特征点方向进行赋值,包括以下内容:
计算每个特征点的模m(x,y)和方向θ(x,y)的计算公式如下:
式中,l(x,y)是特征点所在的尺度空间值;以特征点为中心,计算以3σ为半径的邻域梯度幅值和方向,利用直方图统计领域像素点的梯度方向分布,将梯度方向从0~360°分成36等份,把直方图中主峰值的所在方向确定为该特征点的主方向,若存在一个等份方向的峰值大于主峰值80%时,则将该等份方向作为该特征点的辅方向;
生成特征点描述子,包括以下内容:
将该区域的图像坐标轴旋转到与特征点梯度方向重合;
然后取以特征点为中心的16×16大小的邻域窗口,计算每一个像素的梯度,距离特征点越近则权值越大,子区域的像素梯度按照σ=d/2来进行高斯加权,权重计算公式如下;
再将该区域细分为4×4个子区域,通过对梯度值进行加权的方法统计每个子区域中8个方向的梯度直方图,形成8位的向量描述,这将构成一个4×4×8=128维的描述向量,最后将这些向量归一化后,即为尺度不变特征变换SIFT特征点的描述符。
作为优选技术措施:
所述第七步中,匹配点的筛选方法如下:
采用查询索引KD树方式,通过空间距离选择与图像I1中每个特征点最接近的图像I2中的n个特征匹配点形成一对多的假设匹配集,并采取空间距离作为特征点相似性的评估指标,空间距离包括欧氏距离Dp和像素坐标距离Dd的加权和,其计算公式如下:
其中,E1和E2分别为两个特征点的特征描述子向量;
G1和G2为像素点位置坐标;
计算每个特征点与假设匹配集点的空间距离Ds,其计算公式如下:
Ds=α·Dd+(1-α)·Dp (14)
其中,α为匹配点的像素坐标距离项影响因子;
匹配筛选策略采用最近邻空间距离与次近邻空间距离之比,其计算公式如下:
r=min_fst/min_scd (15)
其中min_fst为最近空间距离,min_scd次近空间距离;当r<T,T是预先确定的比值,则两个角点互为匹配点对,实现对所有的特征点进行粗匹配;
对于局部相似属性的特征点,匹配点的筛选过程如下:
计算匹配点对间的距离,将距离值的最大最小值均匀的分为10个区间,每个区间的频率为P={p1,…,p10},则峰值区间的频率为max(P),对应的区间为第i个,在区间[i-1,i+1]中的匹配点对是正确匹配点对,该匹配点对集为寻找的精确匹配点对;
再根据随机抽样一致RANSAC算法消除错误特征点对,从而计算出单应性变化矩阵,用单应性变换矩阵乘以图像得到待融合图像。
作为优选技术措施:
所述第八步中,基于拉普拉斯Laplacian金字塔的多分辨率融合模型的构建方法如下:
每两幅图像的融合首先建立两幅图像I1和I2的高斯金字塔G1、G2,然后建立对应的4层Laplace金字塔图像Lap1、Lap2;
制作一个和图像I1同样大小的mask掩膜图像Imask,这个掩模图像代表融合的位置,然后求mask图像的高斯金字塔Gmask,其代表每个像素点的融合权重;
在每个尺度即分辨率下,根据当前尺度的Gmask将两幅图像的Laplace金字塔图像Lap1,Lap2进行相加,最终得到拼接的Laplace金字塔图像Lapfused;
以Lapfused的最低分辨率图作为起始图,重构得到最高分辨率的拼接结果。
为实现上述目的之一,本发明的第二种技术方案为:
一种计算机设备,包括:
一个或多个处理器;
存储装置,用于存储一个或多个程序;
当所述一个或多个程序被所述一个或多个处理器执行时,使得所述一个或多个处理器实现上述的一种面向水域的无人机遥感多光谱图像拼接方法。
本发明通过处理器实现上述的一种面向水域的无人机遥感多光谱图像拼接方法,能够有效提高图像拼接精度,方案简单,实用,便于推广使用。
与现有技术相比,本发明具有以下有益效果:
本发明经过不断探索以及试验,根据经纬度POS信息进行图像的几何校正预处理,得到校正图像,并对校正图像添加地理坐标投影信息;然后根据地理坐标投影信息,进行图像的粗配准,得到配准图像;再根据配准图像,构建计算机图形学知识模型,计算出重叠区域位置,并抠取重叠区域图像,得到待拼接图像;进而构建基于主成分图像的尺度不变特征变换SIFT模型,提取待拼接图像的特征点;并以空间距离作为特征点的相似性评估指标,筛选出符合要求的匹配点;最后构建基于拉普拉斯Laplacian金字塔的多分辨率融合模型,对待拼接图形进行融合,实现无人机遥感多光谱图像的拼接。
进而,本发明特征点融合POS信息的方案,能够有效融合了特征点方法配准精准度高的特点,以及位姿拼接法配准速度快的特点,拼接方案可作为实时配准的基础,并降低了对POS信息的精度要求;同时特征点结合经纬度POS信息可以降低大规模图像拼接中出现的误差累积,从而使得拼接图配准精度明显提高,最终得到的拼接融合图像不仅含有地理坐标信息,而且具有良好视觉效果,满足一般的无人机水质遥感监测时的图像拼接要求,方案科学合理,切实可行。
进一步,本发明适用于渔业养殖水域环境,在特征匹配算法中,针对养殖水域图像检测的特征点幅值偏小,方向不明确,以及在图像中存在许多具有局部相似属性的特征点,提供了更好的剔除方法,使得获得的单应性变化矩阵更准确。
更进一步,本发明能够处理长条状,规划航线等图像序列的拼接,且最终的拼接图具有地理坐标。
附图说明
图1是本发明实施例方法的流程图;
图2是本发明各种坐标系关系示意图;
图3是本发明无人机俯仰角变化时沿无人机纵轴剖面的成像几何示意图;
图4是本发明无人机俯仰角变化时三维的成像几何示意图;
图5是本发明无人机俯仰角和滚动角变化时像平面位置几何关系示意图;
图6是本发明无人机偏航角变化时三维的成像几何示意图;
图7是本发明图像仿射矩阵参数计算示意图;
图8是本发明图像根据地理坐标粗配准计算示意图;
图9是本发明DOG尺度空间的构建过程示意图;
图10是本发明旋转关键点邻域坐标轴示意图;
图11是本发明SIFT特征描述子构造示意图
图12是本发明图像特征点匹配连线图;
图13是本发明图像帧到帧的拼接策略示意图;
图14是本发明图像拼接图到拼接图的拼接策略示意图;
图15是本发明建立高斯金字塔和拉普拉斯金字塔Laplace示意图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白,以下结合附图及实施例,对本发明进行进一步详细说明。应当理解,此处所描述的具体实施例仅仅用以解释本发明,并不用于限定本发明。
相反,本发明涵盖任何由权利要求定义的在本发明的精髓和范围上做的替代、修改、等效方法以及方案。进一步,为了使公众对本发明有更好的了解,在下文对本发明的细节描述中,详尽描述了一些特定的细节部分。对本领域技术人员来说没有这些细节部分的描述也可以完全理解本发明。
本发明面向水域的无人机遥感多光谱图像拼接方法的第一种具体实施例:
一种面向水域的无人机遥感多光谱图像拼接方法,
包括以下步骤:
第一步,接收图像数据,并合成多波段图像;
同时获取多波段图像对应的经纬度POS信息,并剔除经纬高数据异常值;
第二步,根据第一步中对应的经纬度POS信息进行图像的几何校正预处理,得到校正图像;
第三步,基于图像对应的经纬和偏航角信息,计算地理坐标投影信息,并对第二步中的校正图像添加地理坐标投影信息;
第四步,根据第三步中的地理坐标投影信息,进行图像的粗配准,得到配准图像;
第五步,根据第四步中的配准图像,构建计算机图形学知识模型,计算出重叠区域位置,并抠取重叠区域图像,得到待拼接图像;
第六步,构建基于主成分图像的尺度不变特征变换SIFT模型,提取第五步中待拼接图像的特征点;
第七步,采用查询索引KD树方式,建立一对多假设匹配集,并以空间距离作为第六步中特征点的相似性评估指标,筛选出符合要求的匹配点;并使用距离直方图约束全局相似属性来剔除误匹配点对,计算单应性变换矩阵;
第八步,根据第七步中的匹配点,按照拼接策略选取无人机图像序列,构建基于拉普拉斯Laplacian金字塔的多分辨率融合模型,对待拼接图形进行融合,实现无人机遥感多光谱图像的拼接。
如图1、图2所示,本发明面向水域的无人机遥感多光谱图像拼接方法的第二种具体实施例:
一种面向水域的无人机遥感多光谱图像拼接方法,
包括以下步骤:
步骤(1).接收数据合成多波段图像,获取图像对应的POS信息,并剔除经纬高数据异常值;
步骤(2).根据对应的POS信息进行图像的几何校正预处理;
步骤(3).基于图像对应的经纬和偏航角信息,计算地理坐标投影并给图像添加地理坐标投影信息;
步骤(4).根据相互配准的两幅图像的地理坐标信息,进行图像的粗配准;
步骤(5).根据粗配准结果,利用计算机图形学知识计算出重叠区域位置,并抠取重叠区域图像;
步骤(6).提取待拼接的两幅图像的基于第一主成分图像的SIFT特征点;
步骤(7).采用KD树索引方式,建立一对多假设匹配集,并以空间距离作为特征点相似性的评估指标,筛选出最佳匹配点,进一步使用距离直方图约束全局相似属性来剔除误匹配点对,之后计算单应性变换矩阵;
步骤(8).按照拼接策略选取无人机图像序列,每两幅图像之间采用4层Laplacian金字塔的多分辨率融合算法。
步骤(1)的具体方式为:
接收无人机遥感的多光谱图像,并将同时刻拍摄的同一位置的单波段图像进行合成,根据获取图像对应的经纬度POS信息剔除起飞和返航时的经纬高数据异常值,所属经纬度POS信息包括经(lon)维(lay)高(H)、航向角(γ)、俯仰角(α)和滚转角(β),以及测量***的地理坐标系和投影坐标系。为了得到地面分辨率(∈),还需获取相机的焦距(f)和传感器的像元尺寸(pixelsize),其计算公式如下。为了方便后续的几何校正,还需知道图像的长(w)和宽(h);
∈=pixelsize*H/f (1)
步骤(2)的具体方式为:
在几何校正的过程中,建立图像坐标系O-xy,相机坐标系S-XsYsZs,机体坐标系P-XpYpZp,地理坐标系E-XEYEZE。当图像坐标系选择以像主点为原点,沿飞行方向为y轴的正方向,垂直于飞行方向为x轴的正方向时,根据相机与机体的相对位置关系及无人机的飞行姿态,像点g(x,y)与对应物点G(XE,YE,ZE)的变化关系可通过式(2)转换,而图像坐标系与相机坐标系存在Z轴方向的平移,平移量为焦距f。
[XE,YE,ZE]E T=λREPRPSRSO[x,y]O T=λREPRPS[x,y,-f]O T (2)
假设相机的镜头光心与无人机的质心重合,且相机坐标轴系经过平移后与机体坐标系完全重合,即RPS为单位矩阵I。因此相机成像模型可简化为:
[XE,YE,ZE]E T=λREP[x,y,-f]O T (3)
其中REP=R(H)R(γ)R(β)R(α),λ为比例系数,即λ=H/f。
因此,在正直摄影条件下,遥感图像中像点g(x,y)与校正后的像点g′(x′,y′)满足:
(x′,y′,-f)O T=REP(x,y,-f)O T (4)
其中,R(α)、R(β)、R(γ)、R(H)分别为基于俯仰角的校正旋转矩阵、基于滚动角的校正旋转矩阵、基于偏航角角的校正旋转矩阵、基于航高的校正矩阵。在无人机其他飞行姿态参数不变的情况下,仅考虑某一姿态变化时的情形,可对R(α)、R(β)、R(γ)、R(H)分别建立如图3、图4、图5和图6所示的几何关系,然后计算得到如下的校正数学模型。
最后,计算图像四角坐标校正后的新坐标,通过旧坐标与新坐标可求出变换矩阵,然后采用双线性插值重采样方法计算出新的校正图像,完成几何校正。
步骤(3)的具体方式为:
仿射矩阵参数可以用来对遥感图像坐标与地理坐标进行转换,其包括6个参数(XE,Xpixel,Rγ,YE,Ypixel,Rγ),描述的是图像行列号和地理坐标之间的关系,其中XE、YE表示图像左上角像元的地理投影坐标,Xpixel、Ypixel分别表示图像像元在经度、维度方向的地面分辨率,Rγ表示图像旋转角度的正弦值。如图7所示,将图像绕O点顺时针旋转γ度,然后根据无人机记录的经纬度经过坐标投影计算可得到O点坐标(XOE,YOE),则图像左上角像点G坐标可通过式(13)获得,而地面分辨率Xpixel=-Ypixel=-∈,此时Rγ=0。其中w和h分别表示图像尺寸的宽和高。
而图像坐标系下任一点g(row,col)的地理坐标G(XE′,YE′)可以通过式(14)计算得到。
步骤(4)的具体方式为:
如图8所示,假设同一航线上相邻的图像I1和图像I2,图像左上角的地理坐标分别为G1(XE1,YE1)和G2(XE2,YE2),以图像I1为基准,可求得图像I2相对于图像I1的偏移量,公式如下:
将几何校正后的图像I1和图像I2分别绕O1和O2旋转γ度至航线方向,其仍是标准的正直摄影,视轴相互平行,且相机坐标系中心到大地平面的距离相等,因此图像I1和图像I2的拼接可认为是刚体平移变化,图像I2中的任一点可通过式(16)计算其在图像I1中的图像坐标。
其配准图像的仿射矩阵参数的左上角像元地理坐标G(XE new,YE new)可根据式(17)求得,配准结果图像的地面分辨率Xpixel、Ypixel和Rγ均保持不变。
步骤(5)的具体方式为:
将图像I1和I2分别绕O1和O2旋转γ度至航线方向后,然后可根据计算机图形学知识,计算出图像I1四角构成的四边形I1AI1BI1CI1D与图像I2四角构成的四边形I2AI2BI2CI2D的多边形重叠区域,在重叠区域向外添加一定偏置δ,最终获得多边形重叠区域ABCD的坐标,如图8所示。之后,利用多边形重叠区域ABCD的坐标制作掩膜图像,就可以从图像I1和I2中获得只含重叠区域的图像,记为图像I1′和I2′。
步骤(6)的具体方式为:
1).提取多光谱图像第一主成分图像的具体步骤如下:
a)输入大小为w×h,波段数为c的多光谱图像矩阵I,将I重组为w×h行c列矩阵Ireshape,对Ireshape每一列进行归一化得到I′reshape:
b)求I′reshape的协方差矩阵Cov:
c)求协方差矩阵Cov最大特征值对应的特征向量Vec;
d)对重组图像Ireshape提取特征向量Vec的主成分图像,并重组为w×h大小的图像IPCA。
IPCA=Ireshape·Vec (20)
2).将PCA提取的第一主成分图像作为SIFT特征点提取的输入,然后进行构建DOG尺度空间、检测DOG尺度空间的极值点,删除不稳定的特征点、对特征点方向进行赋值和生成特征点描述子五个小步骤,具体实施如下:
构建DOG尺度空间:
通过对源图像进行高斯模糊和降采样得到图像金字塔,降采样过程可以用式(21)表示;其次对图像金字塔的每层图像使用不同的σ参数进行高斯模糊,得到的多张模糊过图像构成了高斯金字塔,可用式(22)和式(23)表示;然后将相邻的两个高斯空间的图像相减就得到了DOG图像,可用式(24)表示,整个构建过程可以用图9表示。
式中,Down表示降采样,其中G0=I;I(x,y)表示源图像,L(x,y,σ)表示原图像卷积后的高斯尺度空间,代表卷积运算,σ表示高斯卷积核的尺度因子;G(x,y,σ)表示高斯核函数,具体如式(24)所示,(m,n)代表卷积核的大小;k表示相邻尺度空间的比例因子,一般取
检测DOG尺度空间的极值点:
将每个像素点与同一尺度领域的8个像素点和上下相邻的尺度的18个像素点进行比较,只有当该像素点的DOG值全部大于或小于对比的26个像素点的DOG值时,即认为该点是DOG尺度空间上一个极值点。
删除不稳定的特征点:
首先,为了得到极值点的精确位置,需要对离散空间进行像素差值,然后通过类比拟合三维二次函数以得到极值点精确位置,最后剔除低对比度的极值点;另外,对边缘响应比较大导致不太稳定的极值点也进行舍弃。
对特征点方向进行赋值:
为了使特征点具有旋转不变性,需要确定特征点的最大梯度方向,因此在高斯图像中,定义每个特征点的模m(x,y)和方向θ(x,y)如下:
式中,L(x,y)是特征点所在的尺度空间值。以特征点为中心,计算以3σ为半径的邻域梯度幅值和方向,利用直方图统计领域像素点的梯度方向分布,将梯度方向从0~360°分成36等份,把直方图中主峰值的所在方向确定为该特征点的主方向,若存在一个等份方向的峰值大于主峰值80%时,则将该等份方向作为该特征点的辅方向。
生成特征点描述子:
为了确保描述子具有旋转不变性,需要将该区域的图像坐标轴旋转到与特征点梯度方向重合,如图10所示,红色箭头表示特征点梯度方向。
然后取以特征点为中心的16×16大小的邻域窗口,计算每一个像素的梯度,距离特征点越近则权值越大,根据Lowe的研究,子区域的像素梯度按照σ=d/2来进行高斯加权,权重计算公式如下。
再将该区域细分为4×4个子区域,通过对梯度值进行加权的方法统计每个子区域中8个方向的梯度直方图,形成8位的向量描述,这将构成一个4×4×8=128维的描述向量,最后将这些向量归一化后,即为SIFT特征点的描述符,如图11所示。
步骤(7)的具体方式为:
采用KD树索引方式,通过欧式距离选择与图像I1中每个特征点最接近的图像I2中的n个特征匹配点形成一对多的假设匹配集,采取空间距离作为特征点相似性的评估指标,假设两个特征点的特征描述向量分别为E1和E2,位置坐标为G1和G2,那么欧氏距离Dp和像素坐标距离Dd的计算公式如下:
为了限制匹配点不在全局匹配中搜索匹配点,在计算空间距离时引入了匹配点间的像素坐标距离项,影响因子为α,然后计算每个特征点与假设匹配集点的空间距离Ds,其计算公式如下:
Ds=α·Dd+(1-α)·Dp (14)
匹配筛选策略采用最近邻空间距离与次近邻空间距离之比,即:
r=min_fst/min_scd (15)
其中min_fst为最近空间距离,min_scd次近空间距离。当r<T(T是预先确定的比值),则两个角点互为匹配点对,实现对所有的特征点进行粗匹配。
上面的粗匹配仅约束了全局匹配,但是对于局部相似属性的特征点,可以假设在粗匹配集中有许多正确的匹配对,错误匹配是随机分布的,那么匹配点对的距离直方图中将会出现一个峰值,则正确匹配将分布在峰值的周围。具体步骤可表示为:计算匹配点对间的距离,将距离值的最大最小值均匀的分为10个区间,每个区间的频率为P={p1,…,p10},则峰值区间的频率为max(P),对应的区间为第i个,在区间[i-1,i+1]中的匹配点对认为是正确匹配点对,该匹配点对集为寻找的精确匹配点对。最后再根据RANSAC算法进一步消除错误特征点对,从而计算出单应性变化矩阵,用单应性变换矩阵乘以图像得到待融合图像。根据上述方法,生成的特征匹配点对的连线如图12所示。
步骤(8)的具体方式为:
在拼接的前4轮使用帧到帧拼接策略,后续选择拼接图到拼接图的拼接策略以完成图像的拼接,帧到帧的拼接策略如图13所示,拼接图到拼接图的拼接策略如图14所示。每两幅图像的融合首先建立两幅图像I1和I2的高斯金字塔G1、G2,然后建立对应的4层Laplace金字塔图像Lap1、Lap2。制作一个和图像I1同样大小的mask掩膜图像Imask,这个掩模图像代表了融合的位置,然后求mask图像的高斯金字塔Gmask,其代表每个像素点的融合权重。在每个尺度(分辨率)下,根据当前尺度的Gmask将两幅图像的Laplace金字塔图像Lap1,Lap2进行相加,最终得到拼接的Laplace金字塔图像Lapfused。以Lapfused的最低分辨率图作为起始图,重构得到最高分辨率的拼接结果,整个构建构成如图15所示。
应用本发明方法的一种装置实施例:
一种计算机设备,其包括:
一个或多个处理器;
存储装置,用于存储一个或多个程序;
当所述一个或多个程序被所述一个或多个处理器执行时,使得所述一个或多个处理器实现上述的一种面向水域的无人机遥感多光谱图像拼接方法。
应用本发明方法的一种计算机介质实施例:
一种计算机可读存储介质,其上存储有计算机程序,该程序被处理器执行时实现上述的一种面向水域的无人机遥感多光谱图像拼接方法。
本领域内的技术人员应明白,本申请的实施例可提供为方法、***、或计算机程序产品。因此,本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且,本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本申请是参照根据本申请实施例的方法、设备(***)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器,使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中,使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品,该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上,使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理,从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
最后应当说明的是:以上实施例仅用以说明本发明的技术方案而非对其限制,尽管参照上述实施例对本发明进行了详细的说明,所属领域的普通技术人员应当理解:依然可以对本发明的具体实施方式进行修改或者等同替换,而未脱离本发明精神和范围的任何修改或者等同替换,其均应涵盖在本发明的权利要求保护范围之内。
Claims (10)
1.一种面向水域的无人机遥感多光谱图像拼接方法,其特征在于,
包括以下步骤:
第一步,接收图像数据,并合成多波段图像;
同时获取多波段图像对应的经纬度POS信息,并剔除经纬高数据异常值;
第二步,根据第一步中对应的经纬度POS信息进行图像的几何校正预处理,得到校正图像;
第三步,基于图像对应的经纬和偏航角信息,计算地理坐标投影信息,并对第二步中的校正图像添加地理坐标投影信息;
第四步,按照拼接策略选取无人机图像序列,根据第三步中的地理坐标投影信息,进行图像的粗配准,得到平移变换矩阵和两幅配准图像;
第五步,根据第四步中的配准图像,利用计算机图形学知识,计算出重叠区域位置,并抠取重叠区域图像,得到待拼接图像;
第六步,利用主成分分析法PCA对第五步中的待拼接图像提取第一主成分图像,并对第一主成分图像构建尺度不变特征模型SIFT,以获取第一主成分图像的特征点;
第七步,判断第六步中的特征点数量,若数量少于N个,将第四步中的平移变换矩阵作为后期的单应性变化矩阵;
否则采用查询索引KD树方式,建立一对多假设匹配集,并以空间距离作为第六步中特征点的相似性评估指标,筛选出符合要求的匹配点;并使用距离直方图约束全局相似属性来剔除误匹配点对,计算单应性变换矩阵;
第八步,对第四步中的两幅配准图像构建基于拉普拉斯Laplacian金字塔的多分辨率融合模型,并结合第七步中的单应性变换矩阵,对第四步中的两幅配准图像进行融合,实现无人机遥感多光谱图像的拼接;
之后返回第四步,直到完成所有序列图像的拼接后退出。
2.如权利要求1所述的一种面向水域的无人机遥感多光谱图像拼接方法,其特征在于,
所述第一步中,多波段图像的合成方法如下:
接收无人机遥感的多光谱图像,并将同时刻拍摄的同一位置的单波段图像进行合成,得到多波段图像;
并根据获取单波段图像对应的经纬度POS信息剔除起飞和返航时的经纬高数据异常值;
所述经纬度POS信息包括经度lon、纬度lat、航高H、航向角γ、俯仰角α、滚转角β、地面分辨率∈以及测量***的地理坐标系和投影坐标系;
地面分辨率∈的计算公式如下:
∈=pixelsize*H/f (1)。
3.如权利要求1所述的一种面向水域的无人机遥感多光谱图像拼接方法,其特征在于,
所述第二步中,几何校正的方法包括以下内容:
步骤21,建立图像坐标系O-xy,相机坐标系S-XsYsZs,机体坐标系P-XpYpZp,地理坐标系E-XEYEZE;
步骤22,对步骤21中的图像坐标系选择以像主点为原点,沿飞行方向为y轴的正方向,垂直于飞行方向为x轴的正方向时,根据相机与机体的相对位置关系及无人机的飞行姿态,进行图像坐标系下像点g(x,y)与地理坐标系下对应物点G(XE,YE,ZE)的变化关系转换;
变化关系转换的计算公式如下:
[XE,YE,ZE]E T=λREPRPSRSO[x,y]O T=λREPRPS[x,y,-f]o T (2)
其中,f为焦距,用于表示图像坐标系与相机坐标系在Z轴方向的平移量;
步骤23,由于相机的镜头光心与无人机的质心重合,且相机坐标轴系经过平移后与机体坐标系完全重合,即RPS为单位矩阵I;
根据步骤22中的平移量,构建相机成像模型,其计算公式如下:
[XE,YE,ZE]E T=λREP[x,y,-f]O T (3)
其中,REP=R(H)R(γ)R(β)R(α),λ为比例系数,即λ=H/f;
在正直摄影条件下,遥感图像中像点g(x,y)与校正后的像点g′(x′,y′)的满足关系式如下:
(x′,y′,-f)O T=REP(x,y,-f)O T (4)
其中,R(α)、R(β)、R(γ)、R(H)分别为基于俯仰角的校正旋转矩阵、基于滚动角的校正旋转矩阵、基于偏航角角的校正旋转矩阵、基于航高的校正矩阵;
步骤24,对步骤23中的R(α)、R(β)、R(γ)、R(H)分别建立俯仰角校正矩阵、滚动角校正矩阵、偏航角校正矩阵、高度校正矩阵,得到校正数学模型;
其中,θ表示任一像素点和焦点的连线与中心视轴之间的夹角;
步骤25,根据步骤24中的校正数学模型,计算图像四角坐标校正后的新坐标,通过旧坐标与新坐标求出变换矩阵,并采用双线性插值重采样方法计算出新的校正图像,完成图像的几何校正,并更新图像的地理坐标信息。
4.如权利要求1所述的一种面向水域的无人机遥感多光谱图像拼接方法,其特征在于,所述第三步中,地理坐标投影信息的计算方法如下:
利用仿射矩阵参数对遥感图像坐标与地理坐标进行转换,其包括6个参数,分别为XE,Xpixel,Rγ,YE,Ypixel,Rγ,描述的是图像行列号和地理坐标之间的关系,
其中,XE、YE表示图像左上角像元的地理投影坐标,Xpixel、Ypixel分别表示图像像元在经度、维度方向的地面分辨率,Rγ表示图像旋转角度的正弦值;
将图像绕中心O点顺时针旋转γ度,然后根据无人机记录的经纬度经过坐标投影计算得到O点坐标(XOE,YOE),则图像左上角像点G坐标的计算公式如下:
其中w和h分别表示图像尺寸的宽和高;
地面分辨率Xpixel=-Ypixel=-∈,此时Rγ=0;
图像坐标系下任一点g(row,col)的地理坐标G(XE′,YE′)的计算公式如下:
5.如权利要求1所述的一种面向水域的无人机遥感多光谱图像拼接方法,其特征在于,
所述第四步中,粗配准的过程如下:
在拼接的前几轮使用帧到帧拼接策略,后续选择拼接图到拼接图的拼接策略以完成图像的拼接;
获取同一航线上相邻的图像I1和图像I2,图像左上角的地理坐标分别为G1(XE1,YE1)和G2(XE2,YE2),以图像I1为基准,求得图像I2相对于图像I1的偏移量;
所述偏移量的计算公式如下:
根据偏移量,将几何校正后的图像I1和图像I2分别绕O1和O2旋转γ度至航线方向,图像I2中的任一点在图像I1中的图像坐标的计算公式如下:
其中,Hrigid表示平移变化矩阵;
配准图像的仿射矩阵参数的左上角像元地理坐标G(XE new,YE new)的计算公式如下:
配准图像的地面分辨率Xpixel、Ypixel和Rγ均保持不变。
6.如权利要求1所述的一种面向水域的无人机遥感多光谱图像拼接方法,其特征在于,
所述第五步中,抠取重叠区域图像,包括以下内容:
步骤51,将图像I1和I2分别绕O1和O2旋转γ度至航线方向;
步骤52,步骤51中的旋转完成后,根据计算机图形学知识,计算出图像I1四角构成的四边形I1AI1BI1CI1D与图像I2四角构成的四边形I2AI2BI2CI2D的多边形重叠区域;
步骤53,在步骤52中的多边形重叠区域向外添加一定偏置δ,最终获得多边形重叠区域ABCD的坐标;
步骤54,利用步骤53中的多边形重叠区域ABCD的坐标制作掩膜图像,从图像I1和I2中获得只含重叠区域的图像,记为图像I1′和I2′。
7.如权利要求1所述的一种面向水域的无人机遥感多光谱图像拼接方法,其特征在于,
所述第六步中,构建第一主成分图像的尺度不变特征变换SIFT模型的方法如下:
步骤61,输入大小为w×h,波段数为c的多光谱图像矩阵I,将I重组为w×h行c列矩阵Ireshape,对Ireshape每一列进行归一化得到I′reshape;
I′reshape的计算公式如下:
步骤62,根据步骤61中的I′reshape,计算协方差矩阵Cov;
协方差矩阵Cov的计算公式如下:
步骤63,计算步骤62中的协方差矩阵Cov最大特征值对应的特征向量Vec;
步骤64,对重组图像的Ireshape提取步骤63中的特征向量Vec的主成分图像,并重组为w×h大小的图像IPCA;
图像IPCA的计算公式如下:
IPCA=Ireshape·Vec (20)
步骤65,将步骤64中的主成分图像作为尺度不变特征变换SIFT特征点提取的输入,然后进行构建DOG尺度空间、检测DOG尺度空间的极值点,删除不稳定的特征点、对特征点方向进行赋值和生成特征点描述子;
所述构建DOG尺度空间,包括以下内容:
通过对源图像进行高斯模糊和降采样得到图像金字塔,降采样的计算公式如下:
其次对图像金字塔的每层图像使用不同的σ参数进行高斯模糊,得到的多张模糊过图像构成了高斯金字塔,其计算公式如下:
然后将相邻的两个高斯空间的图像相减得到DOG图像,其计算公式如下:
式中,Down表示降采样,其中G0=I;I(x,y)表示源图像,L(x,y,σ)表示原图像卷积后的高斯尺度空间,代表卷积运算,σ表示高斯卷积核的尺度因子;G(x,y,σ)表示高斯核函数,(m,n)代表卷积核的大小;k表示相邻尺度空间的比例因子,取
检测DOG尺度空间的极值点,包括以下内容:
将每个像素点与同一尺度领域的8个像素点和上下相邻的尺度的18个像素点进行比较,只有当该像素点的DOG值全部大于或小于对比的26个像素点的DOG值时,该像素点是DOG尺度空间上一个极值点;
删除不稳定的特征点,包括以下内容:
首先,为了得到极值点的精确位置,需要对离散空间进行像素差值,然后通过类比拟合三维二次函数以得到极值点精确位置,最后剔除低对比度的极值点;
同时,对边缘响应大的极值点进行舍弃;
对特征点方向进行赋值,包括以下内容:
计算每个特征点的模m(x,y)和方向θ(x,y)的计算公式如下:
式中,L(x,y)是特征点所在的尺度空间值;以特征点为中心,计算以3σ为半径的邻域梯度幅值和方向,利用直方图统计领域像素点的梯度方向分布,将梯度方向从0~360°分成36等份,把直方图中主峰值的所在方向确定为该特征点的主方向,若存在一个等份方向的峰值大于主峰值80%时,则将该等份方向作为该特征点的辅方向;
生成特征点描述子,包括以下内容:
将该区域的图像坐标轴旋转到与特征点梯度方向重合;
然后取以特征点为中心的16×16大小的邻域窗口,计算每一个像素的梯度,距离特征点越近则权值越大,子区域的像素梯度按照σ=d/2来进行高斯加权,权重计算公式如下;
再将该区域细分为4×4个子区域,通过对梯度值进行加权的方法统计每个子区域中8个方向的梯度直方图,形成8位的向量描述,这将构成一个4×4×8=128维的描述向量,最后将这些向量归一化后,即为尺度不变特征变换SIFT特征点的描述符。
8.如权利要求1所述的一种面向水域的无人机遥感多光谱图像拼接方法,其特征在于,
所述第七步中,匹配点的筛选方法如下:
采用查询索引KD树方式,通过空间距离选择与图像I1中每个特征点最接近的图像I2中的n个特征匹配点形成一对多的假设匹配集,并采取空间距离作为特征点相似性的评估指标,空间距离包括欧氏距离Dp和像素坐标距离Dd的加权和,其计算公式如下:
其中,E1和E2分别为两个特征点的特征描述子向量;
G1和G2为像素点位置坐标;
计算每个特征点与假设匹配集点的空间距离Ds,其计算公式如下:
Ds=α·Dd+(1-α)·Dp (14)
其中,α为匹配点的像素坐标距离项影响因子;
匹配筛选策略采用最近邻空间距离与次近邻空间距离之比,其计算公式如下:
r=min_fst/min_scd (15)
其中min_fst为最近空间距离,min_scd次近空间距离;当r<T,T是预先确定的比值,则两个角点互为匹配点对,实现对所有的特征点进行粗匹配;
对于局部相似属性的特征点,匹配点的筛选过程如下:
计算匹配点对间的距离,将距离值的最大最小值均匀的分为10个区间,每个区间的频率为P={p1,…,p10},则峰值区间的频率为max(P),对应的区间为第i个,在区间[i-1,i+1]中的匹配点对是正确匹配点对,该匹配点对集为寻找的精确匹配点对;
再根据随机抽样一致RANSAC算法消除错误特征点对,从而计算出单应性变化矩阵,用单应性变换矩阵乘以图像得到待融合图像。
9.如权利要求1所述的一种面向水域的无人机遥感多光谱图像拼接方法,其特征在于,
所述第八步中,基于拉普拉斯Laplacian金字塔的多分辨率融合模型的构建方法如下:
每两幅图像的融合首先建立两幅图像I1和I2的高斯金字塔G1、G2,然后建立对应的4层Laplace金字塔图像Lap1、Lap2;
制作一个和图像I1同样大小的mask掩膜图像Imask,这个掩模图像代表融合的位置,然后求mask图像的高斯金字塔Gmask,其代表每个像素点的融合权重;
在每个尺度即分辨率下,根据当前尺度的Gmask将两幅图像的Laplace金字塔图像Lap1,Lap2进行相加,最终得到拼接的Laplace金字塔图像Lapfused;
以Lapfused的最低分辨率图作为起始图,重构得到最高分辨率的拼接结果。
10.一种计算机设备,其特征在于,包括:
一个或多个处理器;
存储装置,用于存储一个或多个程序;
当所述一个或多个程序被所述一个或多个处理器执行时,使得所述一个或多个处理器实现如权利要求1-9任一所述的一种面向水域的无人机遥感多光谱图像拼接方法。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210640221.XA CN114936971A (zh) | 2022-06-08 | 2022-06-08 | 一种面向水域的无人机遥感多光谱图像拼接方法及*** |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210640221.XA CN114936971A (zh) | 2022-06-08 | 2022-06-08 | 一种面向水域的无人机遥感多光谱图像拼接方法及*** |
Publications (1)
Publication Number | Publication Date |
---|---|
CN114936971A true CN114936971A (zh) | 2022-08-23 |
Family
ID=82867391
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202210640221.XA Pending CN114936971A (zh) | 2022-06-08 | 2022-06-08 | 一种面向水域的无人机遥感多光谱图像拼接方法及*** |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN114936971A (zh) |
Cited By (12)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN115393196A (zh) * | 2022-10-25 | 2022-11-25 | 之江实验室 | 一种无人机面阵摆扫的红外多序列影像无缝拼接方法 |
CN115861136A (zh) * | 2022-11-17 | 2023-03-28 | 中国科学院空天信息创新研究院 | 一种基于航空遥感***的影像分辨率重建方法 |
CN115909113A (zh) * | 2023-01-09 | 2023-04-04 | 广东博幻生态科技有限公司 | 一种无人机遥感监测调查林业有害生物的方法 |
CN116228539A (zh) * | 2023-03-10 | 2023-06-06 | 贵州师范大学 | 一种无人机遥感图像拼接的方法 |
CN116359836A (zh) * | 2023-05-31 | 2023-06-30 | 成都金支点科技有限公司 | 一种基于超分辨率测向的无人机目标跟踪方法和*** |
CN116543309A (zh) * | 2023-06-28 | 2023-08-04 | 华南农业大学 | 一种作物异常信息获取方法、***、电子设备及介质 |
CN117036666A (zh) * | 2023-06-14 | 2023-11-10 | 北京自动化控制设备研究所 | 基于帧间图像拼接的无人机低空定位方法 |
CN117079166A (zh) * | 2023-10-12 | 2023-11-17 | 江苏智绘空天技术研究院有限公司 | 一种基于高空间分辨率遥感图像的边缘提取方法 |
CN117132913A (zh) * | 2023-10-26 | 2023-11-28 | 山东科技大学 | 基于无人机遥感与特征识别匹配的地表水平位移计算方法 |
CN117173580A (zh) * | 2023-11-03 | 2023-12-05 | 芯视界(北京)科技有限公司 | 水质参数的获取方法及装置、图像处理方法、介质 |
CN117291980A (zh) * | 2023-10-09 | 2023-12-26 | 宁波博登智能科技有限公司 | 一种基于深度学习的单张无人机图像像素定位方法 |
CN117876222A (zh) * | 2024-03-12 | 2024-04-12 | 昆明理工大学 | 一种弱纹理湖泊水面场景下的无人机影像拼接方法 |
-
2022
- 2022-06-08 CN CN202210640221.XA patent/CN114936971A/zh active Pending
Cited By (22)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN115393196B (zh) * | 2022-10-25 | 2023-03-24 | 之江实验室 | 一种无人机面阵摆扫的红外多序列影像无缝拼接方法 |
CN115393196A (zh) * | 2022-10-25 | 2022-11-25 | 之江实验室 | 一种无人机面阵摆扫的红外多序列影像无缝拼接方法 |
CN115861136B (zh) * | 2022-11-17 | 2023-09-19 | 中国科学院空天信息创新研究院 | 一种基于航空遥感***的影像分辨率重建方法 |
CN115861136A (zh) * | 2022-11-17 | 2023-03-28 | 中国科学院空天信息创新研究院 | 一种基于航空遥感***的影像分辨率重建方法 |
CN115909113A (zh) * | 2023-01-09 | 2023-04-04 | 广东博幻生态科技有限公司 | 一种无人机遥感监测调查林业有害生物的方法 |
CN116228539A (zh) * | 2023-03-10 | 2023-06-06 | 贵州师范大学 | 一种无人机遥感图像拼接的方法 |
CN116359836A (zh) * | 2023-05-31 | 2023-06-30 | 成都金支点科技有限公司 | 一种基于超分辨率测向的无人机目标跟踪方法和*** |
CN116359836B (zh) * | 2023-05-31 | 2023-08-15 | 成都金支点科技有限公司 | 一种基于超分辨率测向的无人机目标跟踪方法和*** |
CN117036666A (zh) * | 2023-06-14 | 2023-11-10 | 北京自动化控制设备研究所 | 基于帧间图像拼接的无人机低空定位方法 |
CN117036666B (zh) * | 2023-06-14 | 2024-05-07 | 北京自动化控制设备研究所 | 基于帧间图像拼接的无人机低空定位方法 |
CN116543309A (zh) * | 2023-06-28 | 2023-08-04 | 华南农业大学 | 一种作物异常信息获取方法、***、电子设备及介质 |
CN116543309B (zh) * | 2023-06-28 | 2023-10-27 | 华南农业大学 | 一种作物异常信息获取方法、***、电子设备及介质 |
CN117291980A (zh) * | 2023-10-09 | 2023-12-26 | 宁波博登智能科技有限公司 | 一种基于深度学习的单张无人机图像像素定位方法 |
CN117291980B (zh) * | 2023-10-09 | 2024-03-15 | 宁波博登智能科技有限公司 | 一种基于深度学习的单张无人机图像像素定位方法 |
CN117079166B (zh) * | 2023-10-12 | 2024-02-02 | 江苏智绘空天技术研究院有限公司 | 一种基于高空间分辨率遥感图像的边缘提取方法 |
CN117079166A (zh) * | 2023-10-12 | 2023-11-17 | 江苏智绘空天技术研究院有限公司 | 一种基于高空间分辨率遥感图像的边缘提取方法 |
CN117132913A (zh) * | 2023-10-26 | 2023-11-28 | 山东科技大学 | 基于无人机遥感与特征识别匹配的地表水平位移计算方法 |
CN117132913B (zh) * | 2023-10-26 | 2024-01-26 | 山东科技大学 | 基于无人机遥感与特征识别匹配的地表水平位移计算方法 |
CN117173580B (zh) * | 2023-11-03 | 2024-01-30 | 芯视界(北京)科技有限公司 | 水质参数的获取方法及装置、图像处理方法、介质 |
CN117173580A (zh) * | 2023-11-03 | 2023-12-05 | 芯视界(北京)科技有限公司 | 水质参数的获取方法及装置、图像处理方法、介质 |
CN117876222A (zh) * | 2024-03-12 | 2024-04-12 | 昆明理工大学 | 一种弱纹理湖泊水面场景下的无人机影像拼接方法 |
CN117876222B (zh) * | 2024-03-12 | 2024-06-11 | 昆明理工大学 | 一种弱纹理湖泊水面场景下的无人机影像拼接方法 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN114936971A (zh) | 一种面向水域的无人机遥感多光谱图像拼接方法及*** | |
CN108648240B (zh) | 基于点云特征地图配准的无重叠视场相机姿态标定方法 | |
CN111583110B (zh) | 一种航拍图像的拼接方法 | |
CN108765298A (zh) | 基于三维重建的无人机图像拼接方法和*** | |
CN103337052B (zh) | 面向宽幅遥感影像的自动几何纠正方法 | |
Zhang et al. | Vision-based pose estimation for textureless space objects by contour points matching | |
CN111507901B (zh) | 基于航带gps及尺度不变约束的航拍图像拼接定位方法 | |
CN107578376B (zh) | 基于特征点聚类四叉划分和局部变换矩阵的图像拼接方法 | |
CN108759788B (zh) | 无人机影像定位定姿方法及无人机 | |
CN113313659B (zh) | 一种多机协同约束下高精度图像拼接方法 | |
CN108765489A (zh) | 一种基于组合靶标的位姿计算方法、***、介质及设备 | |
CN115187798A (zh) | 一种多无人机高精度匹配定位方法 | |
CN109871739B (zh) | 基于yolo-sioctl的机动站自动目标检测与空间定位方法 | |
CN110084743B (zh) | 基于多航带起始航迹约束的图像拼接与定位方法 | |
JP2023530449A (ja) | 空中と地上の位置合わせのためのシステムおよび方法 | |
CN117218201A (zh) | Gnss拒止条件下无人机影像定位精度提升方法及*** | |
CN112132754A (zh) | 一种车辆移动轨迹修正方法及相关装置 | |
CN113963067B (zh) | 一种采用小靶标对大视场视觉传感器进行标定的标定方法 | |
Jhan et al. | A generalized tool for accurate and efficient image registration of UAV multi-lens multispectral cameras by N-SURF matching | |
CN112750075A (zh) | 一种低空遥感影像拼接方法以及装置 | |
CN114897676A (zh) | 一种无人机遥感多光谱图像拼接方法、设备及介质 | |
Zhao et al. | Digital Elevation Model‐Assisted Aerial Triangulation Method On An Unmanned Aerial Vehicle Sweeping Camera System | |
CN117073669A (zh) | 一种飞行器定位方法 | |
Hu et al. | Planetary3D: A photogrammetric tool for 3D topographic mapping of planetary bodies | |
Tjahjadi et al. | Single image orientation of UAV's imagery using orthogonal projection 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 |