CN107274380A - 一种无人机多光谱影像快速拼接方法 - Google Patents

一种无人机多光谱影像快速拼接方法 Download PDF

Info

Publication number
CN107274380A
CN107274380A CN201710552651.5A CN201710552651A CN107274380A CN 107274380 A CN107274380 A CN 107274380A CN 201710552651 A CN201710552651 A CN 201710552651A CN 107274380 A CN107274380 A CN 107274380A
Authority
CN
China
Prior art keywords
image
mrow
msubsup
msub
band
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.)
Granted
Application number
CN201710552651.5A
Other languages
English (en)
Other versions
CN107274380B (zh
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.)
Peking University
Original Assignee
Peking University
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 Peking University filed Critical Peking University
Priority to CN201710552651.5A priority Critical patent/CN107274380B/zh
Publication of CN107274380A publication Critical patent/CN107274380A/zh
Application granted granted Critical
Publication of CN107274380B publication Critical patent/CN107274380B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • 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
    • 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
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/70—Determining position or orientation of objects or cameras
    • 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Image Processing (AREA)

Abstract

本发明提供一种无人机多光谱影像快速拼接方法，包括以下步骤：确定最佳波段和最佳匹配算法；使用确定的最佳波段和最佳匹配算法，对剩余的存在潜在重叠区域的影像对进行特征匹配；以每张影像对应的相机投影中心和航向角为优化参数，对目标函数进行迭代求解，求得每张影像对应的最佳相机投影中心坐标和最佳航向角的值；进而实现多光谱影像拼接。优点为：使用最佳波段描述子选取代表性波段，用于特征点匹配，并将单波段匹配计算得到的外方位元素赋予其余波段，从而保证多波段影像在拼接后，波段之间的独立性与一致性，同时避免了每个波段单独处理产生的冗余计算，或波段分组拼接处理产生的拼接结果之间的不一致问题。

Description

一种无人机多光谱影像快速拼接方法
技术领域
本发明属于影像拼接技术领域，具体涉及一种无人机多光谱影像快速拼接方法。
背景技术
随着无人机(Unmanned aerial vehicle，UAV)技术和多光谱相机的发展，基于无人机的多光谱遥感越来越被人们所关注。以无人机为飞行平台的遥感系统具有低成本、高分辨率和高灵活性等优势，近年得到了快速的发展，可在国土测绘、环境灾害、农情监测等方面得到广泛应用。
由于无人机通常飞行高度较低，所以，单幅影像覆盖面积有限，需要经过拼接才能得到大范围的影像。影像拼接是无人机多光谱影像应用中的重要环节。
现有多光谱影像拼接方法主要包括两种：(1)分组拼接：由于多光谱影像含有的波段数较多，大于3个波段，因此，通常以3个波段为一组，划分为多组波段；然后，对每组波段独立进行拼接；其存在的问题为：由于波段分组后，各组在独立拼接过程中，所提取的特征点、计算得到的匹配结果以及影像的外方位元素，乃至后续影像融合过程中的接缝线的选取等，均会有所差异，从而使得基于波段分组拼接后的影像之间，其波段无法再随意组合，换言之，分组拼接后的结果，只能按组使用，无法随意组合使用，而这一点完全违背了遥感影像处理的初衷；此外，分组拼接时，存在大量的重复性冗余工作，因而拼接过程耗时长，严重影响影像拼接的处理速度；(2)每个波段单独处理，提取每个波段所对应的特征点；然后再对各个波段进行特征匹配，得到最终的影像拼接结构；其存在的问题为：由于某些波段所对应的地物光谱成像较弱，因此，对于一些特殊区域(如农田)，无法提取稳定的特征点，从而造成特征匹配失败；或者，由于不同波段提取的特征点分布不均匀，不同波段计算得到的影像外方位元素不一致，由此也造成拼接后的结果之间无法保持一致性。
另外，影像拼接过程中的相关研究重点通常是保证影像的几何信息准确，例如，利用光束法平差等方法恢复地物的三维信息。而对于多光谱遥感影像而言，由于需要基于地物的光谱特性对地物进行分析，所以在影像拼接过程中，需要尽可能保证光谱特性不受拼接过程影响而失真。以精准农业的多光谱遥感应用为例，多光谱影像主要用于分析作物的生长以及病虫害等相关的光谱信息，一旦光谱信息在影像拼接过程中产生变化，则病害及作物长势的分析得不到准确的结果，因此，多光谱影像拼接过程中，应该尽可能避免对原始影像的重采样或过渡融合引入的原始信息失真等问题。而传统的影像拼接过程的融合阶段，拼接方法主要考虑视觉效果，对相邻重叠区域的多张影像采用融合算法，比如采用加权平均的采样策略，强制使拼接后的影像的像素不再与任一原始影像相同，由此会造成拼接后的影像不再与原始影像一致，从而会损失原始影像的光谱信息。
由此可见，传统的多光谱影像拼接方法，具有多光谱影像拼接效率低、匹配精度有限等问题。
发明内容
针对现有技术存在的缺陷，本发明提供一种无人机多光谱影像快速拼接方法，可有效解决上述问题。
本发明采用的技术方案如下：
本发明提供一种无人机多光谱影像快速拼接方法，包括以下步骤：
步骤1，获取无人机采集的待拼接的n张多光谱遥感影像，同时获取无人机采集所述n张多光谱遥感影像时的航迹数据；
根据所述航迹数据，分析得到各张多光谱遥感影像之间的相邻位置关系，从而确定所有存在潜在重叠区域的影像对，假设存在潜在重叠区域的影像对共有P组；
步骤2，确定最佳波段和最佳匹配算法，步骤如下：
步骤2.1，从P组存在潜在重叠区域的影像对中，挑选出R组存在潜在重叠区域的影像对，分别记为：第1组影像对、第2组影像对...第R组影像对；其中，R远小于P；
步骤2.2，确定需要筛选的波段集合和匹配算法集合，假设波段集合共有a个波段，分别为：第1波段、第2波段...第a波段；假设匹配算法集合共有b个匹配算法，分别为：第1匹配算法、第2匹配算法...第b匹配算法；
步骤2.3，使用不同波段和不同匹配算法的组合进行正交试验，即：使用每种波段-匹配算法组合对挑选出的每组影像对进行特征匹配，并对特征匹配结果进行记录，得到每种波段-匹配算法组合所对应的特征匹配描述子平均值最高值的特征匹配描述子平均值所对应的波段-匹配算法组合即为最佳的波段-匹配算法组合，由此得到最佳波段和最佳匹配算法；
其中，对于任意的第f波段和第g匹配算法，其中，f∈(1，2...a)；g∈(1，2...b)，其所对应的特征匹配描述子平均值基于以下公式计算得到：
其中，Pi代表使用第f波段和第g匹配算法，对第i组影像对进行特征点提取与匹配计算时，得到的第i特征匹配描述子；i∈(1，2...R)；
其中：Pi＝Viri；
其中：Vi＝C(Ii-1,Ii-2)/t；Ii-1和Ii-2分别代表第i组影像对所包含的两个影像；C(Ii-1,Ii-2)代表第i组影像对成功匹配的点对数量；t代表匹配第i组影像对所消耗的时间；Vi代表单位数量匹配点对所耗费的时间；
ri＝C(Ii-1,Ii-2)/mean(Ii-1,Ii-2)；其中，mean(Ii-1,Ii-2)代表第i组影像对提取到的特征点数量的平均值，即：第Ii-1影像提取到的特征点数量和第Ii-2影像提取到的特征点数量的平均值；ri代表第i组影像对的匹配成功率；
步骤3，使用步骤2确定的最佳波段和最佳匹配算法，对剩余的P-R组存在潜在重叠区域的影像对进行特征匹配；
基于步骤2确定的最佳波段和最佳匹配算法，查找步骤2.3的特征匹配结果记录，得到与最佳波段和最佳匹配算法对应的R组存在潜在重叠区域的影像对的特征匹配结果；
由此得到全部P组影像对在最佳波段和最佳匹配算法下的特征匹配结果；
步骤4，对所述特征匹配结果进行分析，识别出所有的同名点；
步骤5，对于待拼接的n张多光谱遥感影像，分别记为：影像1、影像2...影像n；
步骤5.1，令i＝1；
步骤5.2，根据航迹数据，获得影像i对应的原始的相机投影中心和航向角κi的值；还获得对应的航高Hi的值，其中，航高Hi为固定值；由此得到每张影像Vi对应的λi的值；其中，λi为相机焦距fi与航高Hi的比值，为固定值；
步骤5.3，获得影像i的相邻且未处理过的相邻影像，假设相邻且未处理过的相邻影像的数量为wi，将wi个相邻影像依次记为：相邻影像1、相邻影像2...相邻影像wi；
步骤5.4，令v＝1；
步骤5.5，根据步骤4的结果，获得影像i和相邻影像v匹配得到的mv个同名点对的像点坐标，将mv个同名点对的像点坐标按顺序从1开始编号；
步骤5.6，令j＝1；
步骤5.7，按下式计算第j个同名点对所对应的两个地面三维点坐标之间的距离：
其中：代表影像i上的像点所对应的地面二维坐标，并且，该像点与相邻影像v上的某个像点为同名点，并且，该同名点为影像i和相邻影像v之间的第j个同名点；
代表相邻影像v上的像点所对应的地面二维坐标，并且，该像点为与影像i上的某个像点为同名点，并且，该同名点为影像i和相邻影像v之间的第j个同名点；
代表影像i和相邻影像v之间第j个同名点对所对应的地面点对的距离；
表示影像i的像点坐标，并且，该像点与相邻影像v上的某个像点为同名点，并且，该同名点为影像i和相邻影像v之间的第j个同名点；
表示相邻影像v上与对应的同名点的像点坐标；
ki和kv分别为影像i与相邻影像v的航向角；
分别为影像i与相邻影像v的相机投影中心坐标；
步骤5.8，判断j是否等于mv，如果等于，则执行步骤5.9；否则，令j＝j+1，返回步骤5.7；
步骤5.9，判断v是否等于wi，如果等于，则执行步骤5.10；否则，令v＝v+1，返回步骤5.5；
步骤5.10，判断i是否等于n，如果等于，则执行步骤5.11；否则，令i＝i+1，返回步骤5.2；
步骤5.11，对步骤1到步骤5.10所计算得到的地面三维点坐标之间的距离进行求和计算，由此得到如下的目标函数：
步骤6，以每张影像i对应的相机投影中心和航向角κi为优化参数，对步骤5.11的目标函数进行迭代求解，求得每张影像对应的最佳相机投影中心坐标和最佳航向角的值；
步骤7，利用步骤6得到的每张影像的最佳相机投影中心和最佳航向角，对每张影像的所有波段统一进行平移与旋转处理，使每张影像均平移旋转到最佳的位置与姿态；
步骤8，对经过步骤7处理后得到的影像进行融合处理，从而得到最终拼接形成的完整影像。
优选的，步骤2.1中，R组存在潜在重叠区域的影像对为：顺序匹配得到的存在潜在重叠区域的影像对。
优选的，步骤3中，剩余的P-R组存在潜在重叠区域的影像对是指：航带间匹配得到的存在潜在重叠区域的影像对。
优选的，步骤8中，所述融合处理具体为：
当每张影像均平移旋转到最佳的位置与姿态后，对于任意两张存在重叠区域的影像，分别记为影像V1和影像V2；影像V1和影像V2分别具有G个波段，G个波段记为：第1波段、第2波段...第G波段；因此，影像V1和影像V2采用各个对应波段逐一进行融合的方法，即：影像V1的第1波段和影像V2的第1波段进行融合处理；影像V1的第2波段和影像V2的第2波段进行融合处理...依此类推，直到影像V1的第G波段和影像V2的第G波段进行融合处理；
其中，对于任意一个第h波段，其中，h∈(1、2...G)，影像V1的第h波段和影像V2的第h波段采用以下方法进行融合：
步骤8.1，对于影像V1的第h波段和影像V2的第h波段之间的重叠区域，在该重叠区域采样到多个采样点，任意一个采样点Oc的平面坐标为(x,y)；设影像V1的第h波段的中心位置点为OV1，影像V2的第h波段的中心位置点为OV2；
步骤8.2，如果采样点Oc到OV1的直线距离小于采样点Oc到OV2的直线距离，则：Or(x,y)点的像素值为：对影像V1的第h波段在(x,y)位置的像素点进行灰度采样，从而得到Or(x,y)点的像素值；其中，Or(x,y)点为：影像V1的第h波段和影像V2的第h波段进行拼接处理后得到影像V3，影像V3在(x,y)点的表达；
反之，如果采样点Oc到OV1的直线距离大于等于采样点Oc到OV2的直线距离，则：Or(x,y)点的像素值为：对影像V2的第h波段在(x,y)位置的像素点进行灰度采样，从而得到Or(x,y)点的像素值；
步骤8.3，因此，多个采样点Oc对应得到多个Or(x,y)点的像素值，多个Or(x,y)点的像素值即形成拼接影像V3在重叠区域的像素，由此实现影像融合。
本发明提供的一种无人机多光谱影像快速拼接方法具有以下优点：
使用最佳波段描述子选取代表性波段，用于特征点匹配，并将单波段匹配计算得到的外方位元素赋予其余波段，从而保证多波段影像在拼接后，波段之间的独立性与一致性，同时避免了每个波段单独处理产生的冗余计算，或波段分组拼接处理产生的拼接结果之间的不一致问题。
附图说明
图1为本发明提供的一种无人机多光谱影像快速拼接方法的流程示意图；
图2为本发明提供的顺序拼接方法示意图；
图3为本发明提供的航带间拼接方法示意图；
图4为采用本发明提供的多光谱影像快速拼接方法对数据集1拼接得到的影像图；
图5为采用现有技术中的软件Pix4UAV对数据集1进行影像拼接得到的影像图；
图6为采用本发明提供的多光谱影像快速拼接方法对数据集2拼接得到的影像图；
图7为本发明提供的测试区域选取图。
具体实施方式
为了使本发明所解决的技术问题、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
本发明提供一种无人机多光谱影像快速拼接方法，主要特点包括：
(1)当多光谱影像波段数量较多(如大于3个波段)时，本发明使用最佳波段描述子选取代表性波段，用于特征点匹配，并将单波段匹配计算得到的外方位元素(影像的平移量与旋转量)赋予其余波段，由此使一张影像的所有波段使用一致的外方位元素进行融合处理，从而保证多波段影像在拼接后，波段之间的独立性与一致性，同时避免了每个波段单独处理产生的冗余计算，或波段分组拼接处理产生的拼接结果之间的不一致问题。
(2)为了避免多光谱影像在拼接过程中，由于过多的重采样引起的光谱失真问题，将现有光束法平差的三维投影模型简化为二维仿射模型，避免生成三维点云，以及基于三维点云进行重采样计算等复杂处理，同时也回避了由这些计算引入的不规则重采样对原始影像光谱信息所产生的损失。
以上两方面的创新加快并简化了影像拼接过程，具有处理速度快和光谱精度保真度高的优点。
下面对本发明提供的无人机多光谱影像快速拼接方法进行详细介绍：
参考图1，无人机多光谱影像快速拼接方法，包括以下步骤：
步骤1，获取无人机采集的待拼接的n张多光谱遥感影像，同时获取无人机采集所述n张多光谱遥感影像时的航迹数据；
根据所述航迹数据，分析得到各张多光谱遥感影像之间的相邻位置关系，从而确定所有存在潜在重叠区域的影像对，假设存在潜在重叠区域的影像对共有P组；
步骤2，确定最佳波段和最佳匹配算法，步骤如下：
步骤2.1，从P组存在潜在重叠区域的影像对中，挑选出R组存在潜在重叠区域的影像对，分别记为：第1组影像对、第2组影像对...第R组影像对；其中，R远小于P；
本步骤中，可以采用顺序匹配得到的存在潜在重叠区域的影像对。参考图2，即为顺序匹配的示意图。
步骤2.2，确定需要筛选的波段集合和匹配算法集合，假设波段集合共有a个波段，分别为：第1波段、第2波段...第a波段；假设匹配算法集合共有b个匹配算法，分别为：第1匹配算法、第2匹配算法...第b匹配算法；
步骤2.3，使用不同波段和不同匹配算法的组合进行正交试验，即：使用每种波段-匹配算法组合对挑选出的每组影像对进行特征匹配，并对特征匹配结果进行记录，得到每种波段-匹配算法组合所对应的特征匹配描述子平均值特征匹配描述子指标可对影像波段与匹配算法结合的效率和稳定性进行综合评估。较高的值意味着更好的匹配效果，因此，最高值的特征匹配描述子平均值所对应的波段-匹配算法组合即为最佳的波段-匹配算法组合，由此得到最佳波段和最佳匹配算法；
其中，对于任意的第f波段和第g匹配算法，其中，f∈(1，2...a)；g∈(1，2...b)，其所对应的特征匹配描述子平均值基于以下公式计算得到：
其中，Pi代表使用第f波段和第g匹配算法，对第i组影像对进行特征点提取与匹配计算时，得到的第i特征匹配描述子；i∈(1，2...R)；
其中：Pi＝Viri；
其中：Vi＝C(Ii-1,Ii-2)/t；Ii-1和Ii-2分别代表第i组影像对所包含的两个影像；C(Ii-1,Ii-2)代表第i组影像对成功匹配的点对数量；t代表匹配第i组影像对所消耗的时间；Vi代表单位数量匹配点对所耗费的时间，可以对特征匹配的时间消耗进行定量评价；
ri＝C(Ii-1,Ii-2)/mean(Ii-1,Ii-2)；其中，mean(Ii-1,Ii-2)代表第i组影像对提取到的特征点数量的平均值，即：第Ii-1影像提取到的特征点数量和第Ii-2影像提取到的特征点数量的平均值；ri代表第i组影像对的匹配成功率，可反应影像特征匹配的效果；
本步骤的思路可描述为：为了在多光谱影像的多个波段中挑选最佳波段，以及挑选与最佳波段对应的最佳匹配算法，从而用于外方位元素的计算处理。因此，本发明提出了一个特征匹配描述子对影像波段与匹配算法结合的效率和稳定性进行综合评估。但是，特征匹配描述子的计算依赖于在影像中找到的特征点数和成功匹配的点对数量，所以需要在匹配之后进行。而发明人希望在开始匹配之前得到特征匹配描述子的结果，为此引入了一个两阶段的匹配流程，即：第一个阶段，采用少量影像，如只使用航向，即前后相邻影像进行分波段匹配，参考图2，从而得到特征匹配描述子；第二段阶段，使用具有最佳特征匹配描述子的波段和匹配算法，对剩余影像进行特征匹配，例如，进行航带间的匹配。由此实现所有影像的特征匹配。
步骤3，使用步骤2确定的最佳波段和最佳匹配算法，对剩余的P-R组存在潜在重叠区域的影像对进行特征匹配；此处，可参考图3，采用航带间匹配得到的存在潜在重叠区域的影像对。
基于步骤2确定的最佳波段和最佳匹配算法，查找步骤2.3的特征匹配结果记录，得到与最佳波段和最佳匹配算法对应的R组存在潜在重叠区域的影像对的特征匹配结果；
由此得到全部P组影像对在最佳波段和最佳匹配算法下的特征匹配结果；
步骤4，对所述特征匹配结果进行分析，识别出所有的同名点；
步骤5，对于待拼接的n张多光谱遥感影像，分别记为：影像1、影像2...影像n；
步骤5.1，令i＝1；
步骤5.2，根据航迹数据，获得影像i对应的原始的相机投影中心和航向角κi的值；还获得对应的航高Hi的值，其中，航高Hi为固定值；由此得到每张影像Vi对应的λi的值；其中，λi为相机焦距fi与航高Hi的比值，为固定值；
步骤5.3，获得影像i的相邻且未处理过的相邻影像，假设相邻且未处理过的相邻影像的数量为wi，将wi个相邻影像依次记为：相邻影像1、相邻影像2...相邻影像wi；
步骤5.4，令v＝1；
步骤5.5，根据步骤4的结果，获得影像i和相邻影像v匹配得到的mv个同名点对的像点坐标，将mv个同名点对的像点坐标按顺序从1开始编号；
步骤5.6，令j＝1；
步骤5.7，按下式计算第j个同名点对所对应的两个地面三维点坐标之间的距离：
其中：代表影像i上的像点所对应的地面二维坐标，并且，该像点与相邻影像v上的某个像点为同名点，并且，该同名点为影像i和相邻影像v之间的第j个同名点；
代表相邻影像v上的像点所对应的地面二维坐标，并且，该像点为与影像i上的某个像点为同名点，并且，该同名点为影像i和相邻影像v之间的第j个同名点；
代表影像i和相邻影像v之间第j个同名点对所对应的地面点对的距离；
表示影像i的像点坐标，并且，该像点与相邻影像v上的某个像点为同名点，并且，该同名点为影像i和相邻影像v之间的第j个同名点；
表示相邻影像v上与对应的同名点的像点坐标；
ki和kv分别为影像i与相邻影像v的航向角；
分别为影像i与相邻影像v的相机投影中心坐标；
步骤5.8，判断j是否等于mv，如果等于，则执行步骤5.9；否则，令j＝j+1，返回步骤5.7；
步骤5.9，判断v是否等于wi，如果等于，则执行步骤5.10；否则，令v＝v+1，返回步骤5.5；
步骤5.10，判断i是否等于n，如果等于，则执行步骤5.11；否则，令i＝i+1，返回步骤5.2；
步骤5.11，对步骤1到步骤5.10所计算得到的地面三维点坐标之间的距离进行求和计算，由此得到如下的目标函数：
对于本步骤的目标函数，为方便对本发明进行理解，以及，为阐述目标函数的推理过程，进行以下说明：
1)简化的光束法计算模型
无人机在户外采集数据的飞行过程中，相机每拍摄一次影像，其位置与姿态信息，均会在无人机的自驾系统进行记录，即所谓的航迹数据。在无人机影像拼接过程中，一般会采用无人机的航迹数据确定影像之间初始的姿态与位置关系，所以通常在完成影像的特征点匹配之后，直接采用光束法平差计算影像之间精确的空间关系并求解匹配点的三维坐标。一般使用透视变换模型为：
其中：(Xc,Yc,Zc)为相机投影中心在物方空间的坐标，(Xp,Yp,Zp)为地面点的物方空间坐标，(x、y、z)为地面点经K与R变换到像方空间的三维坐标。其中，K矩阵为相机的内参数矩阵；R为旋转矩阵；
设相机在横滚、俯仰、航向三个轴上的旋转角分别为ω和κ，则旋转矩阵R为：
其中：Rω、Rκ分别为横滚、俯仰、航向三个轴上的旋转矩阵。
当同一个地面点在多张影像上存在投影时，使用该透视变换模型则可建立多张影像之间的线性组合关系，通过最小二乘原理，可求解得到精确的影像间旋转关系与匹配点三维坐标。
由于多光谱影像的一个非常重要的应用领域是农作业生长监测与病害的防治，如前所述，如果使用公式(1)透视变换模型，则需要计算大量三维地面点，进行正射影像生成，由此则会引入大量重采样，会给原始的光谱信息带来失真。在顾及大范围农田所处地形一般位于相对平坦的地区，发明人对该透视变换模型进行如下简化：
首先，对于多光谱传感器而言，当其搭载在无人机上时，必须采用云台稳定装置，在天气条件较好的情况下，云台使其保持垂直于地面，此时横滚轴和俯仰轴的角度可以取近似0值，即有ω＝0，这时的旋转矩阵R可近似为：
此时透视变换模型可写为：
由于获取影像的目标区域以较为平常的农田区域为主，考虑到无人机在影像获取过程中，一般其航高设置为固定值，所以，在这种情况下，Zp-Zc可用航高H取代。因此，模型可简化为：
为了精确处理多光谱影像，一般在航拍前会对多光谱相机进行内参数校正，由于其像点坐标相对于焦距值一般很小，同时在顾及对航高方向所做的近似处理，故此时相机内参数矩阵K可近似为：
将Rk和K的详细表达式分别代入式(5)，经整理可得如下表达式：
其中λ＝f/H，其描述的是影像中的坐标系和现实世界坐标系间的尺度的关系。
显然该简化后的模型，计算时需要求解的未知数，对于每张影像的外方位元素而言，由6个(即ω、κ、Xc、Yc、Zc)已经下降至3个(即κ、Xc和Yc)，从而可以极大地提升计算效率。
2)由于在无人机航摄过程中，航迹数据已经记录了每张影像的投影中心位置与航向角等参数，虽然这些值存在误差，但与真实值的差异较小，可作为初始值，输入到目标函数中。因此，在完成前面的所有匹配工作后，不再使用匹配点计算影像之间的相对空间关系(包括相对位移与相对旋转角)，而是直接采用光束法平差，使用公式(7)进行优化计算，所需优化的未知数包括：每张影像的相机投影中心与航向角，由于航高固定为H，因此，地面三维坐标实际为二维坐标(Xp,Yp)。对于地面三维点(Xp,Yp,H)的初始值，由前面公式(7)可计算得到，对公式(7)，利用影像上某个点的像点(x1,y1)可变换为：
由于旋转矩阵为正定阵，其逆矩阵为其转置，故有公式(8)的表达形式，其中相机投影中心坐标(Xc,Yc)与航向角κ的初始值，均取自无人机航迹中的记录数据，λ为相机焦距f与航高H的比值，从而由(8)式中可直接求得未知数(Xp,Yp)。也就是说，通过公式8，可求得每个影像上每个像素点对应的地面三维坐标值。
3)由于任一影像可能具有多达八张具有相邻关系的影像，因而，同一地面点可能会投影到多张影像中，从而构成多个同名点对，使用不同的同名点对，由公式(8)计算得到对应的三维点坐标，由于误差存在的原因，三维点坐标相互之间存在一定的差异，但在理论上，这些三维点坐标应该完全一致，所以，为求得最优的计算结果，应该确保属于同一组同名点的对应的三维坐标点之间的差异达到最小值甚至接近0值。
具体的，假设任一影像M，存在w张相邻影像(w的取值一般为3～8)，以第1张相邻影像为例，假设影像M与第1张相邻影像具有m对同名点(m的取值一般应大于等于6)，设(Xpj1,Ypj1)、(Xpj2,Ypj2)为任意一对同名点计算得到的三维点坐标的平面分量，影像M存在相机投影中心(Xc、Yc)、航向角κ的最佳估值，因此，同一组同名点对所对应的地面三维点坐标的平面分量(Xpj1,Ypj1)、(Xpj2,Ypj2)之间的差异(也即距离)和为最小，即有：
将公式(8)代入到公式(9)中，得到以下更为详细的目标表达式：
步骤6，以每张影像i对应的相机投影中心和航向角κi为优化参数，对步骤5.11的目标函数进行迭代求解，求得每张影像对应的最佳相机投影中心坐标和最佳航向角的值；
步骤7，利用步骤6得到的每张影像的最佳相机投影中心和最佳航向角，对每张影像的所有波段统一进行平移与旋转处理，使每张影像均平移旋转到最佳的位置与姿态；
步骤8，对经过步骤7处理后得到的影像进行融合处理，从而得到最终拼接形成的完整影像。
本步骤具体为：
当每张影像均平移旋转到最佳的位置与姿态后，对于任意两张存在重叠区域的影像，分别记为影像V1和影像V2；影像V1和影像V2分别具有G个波段，G个波段记为：第1波段、第2波段...第G波段；因此，影像V1和影像V2采用各个对应波段逐一进行融合的方法，即：影像V1的第1波段和影像V2的第1波段进行融合处理；影像V1的第2波段和影像V2的第2波段进行融合处理...依此类推，直到影像V1的第G波段和影像V2的第G波段进行融合处理；
其中，对于任意一个第h波段，其中，h∈(1、2...G)，影像V1的第h波段和影像V2的第h波段采用以下方法进行融合：
步骤8.1，对于影像V1的第h波段和影像V2的第h波段之间的重叠区域，在该重叠区域采样到多个采样点，任意一个采样点Oc的平面坐标为(x,y)；设影像V1的第h波段的中心位置点为OV1，影像V2的第h波段的中心位置点为OV2；
步骤8.2，如果采样点Oc到OV1的直线距离小于采样点Oc到OV2的直线距离，则：Or(x,y)点的像素值为：对影像V1的第h波段在(x,y)位置的像素点进行灰度采样，从而得到Or(x,y)点的像素值；其中，Or(x,y)点为：影像V1的第h波段和影像V2的第h波段进行拼接处理后得到影像V3，影像V3在(x,y)点的表达；
反之，如果采样点Oc到OV1的直线距离大于等于采样点Oc到OV2的直线距离，则：Or(x,y)点的像素值为：对影像V2的第h波段在(x,y)位置的像素点进行灰度采样，从而得到Or(x,y)点的像素值；
步骤8.3，因此，多个采样点Oc对应得到多个Or(x,y)点的像素值，多个Or(x,y)点的像素值即形成拼接影像V3在重叠区域的像素，由此实现影像融合。
本步骤进行图像融合的思路可描述为：
在所有的影像均平移旋转到最佳的位置与姿态后，由于影像之间重叠度的存在，相邻影像之间会存在至少两张影像共同覆盖同一个区域的情况，由于拼接后的影像必须是完整的一张影像，因此，需要对多张影像重叠区域，经采样处理到一张完整的拼接影像中去，该过程即为影像的融合。
在融合过程中，为了避免因融合计算产生的原始光谱的失真，不再对重叠区域的影像像素进行加权计算或类似的处理，而是对所有波段逐一进行拼接处理，设拼接后的影像为F，则F上任一点a(x,y)，其像素值g(x,y)取覆盖(x,y)位置的原始影像上的像素点，为确保尽可能少的光谱失真，仅从覆盖当前点位的多张影像中选取一张，且从所选取的影像中进行灰度值采样时，采用较好的双线性采样方法。而原始影像选择，遵从以下方法：
假设当前位置a(x,y)存在2张覆盖该点的影像，分别为影像V1和影像V2，其中心位置点分别为：OV1和OV2；在顾及镜头畸变从影像中心点向外增大的因素，选取影像中心距离a(x,y)位置最小的影像用于像素值的采样。
例如，影像V1和影像V2存在重叠区域，其拼接后影像为拼接影像V3，因此，拼接影像V3上任一点a(x,y)，其像素值g(x,y)取影像V1或影像V2在对应位置的像素点，影像V1或影像V2即为原始影像。而至于取影像V1还是影像V2在对应位置的像素点时，选择原则为：影像中心距离a(x,y)位置最小的影像用于像素值的采样。
综上所述，本发明提供的一种无人机多光谱影像快速拼接方法，具有以下优点：
(1)本发明提出最佳波段描述子，即：从多光谱影像中挑选最佳波段和最佳匹配算法，从而组合为最佳波段描述子；然后，采用最佳波段描述子用于其他影像的特征点提取和后续的匹配计算；此外，将该最佳波段描述子计算得到的外方位元素直接赋值给其它波段影像，用于各波段独立的融合处理，从而保持多光谱影像所有波段拼接后的一致性。
(2)考虑到在农业遥感中，对于地面植被的生长状况分析等完全依赖于影像的光谱信息，所以，在多光谱影像拼接处理过程中，必须优先确保光谱信息尽可能不受拼接过程的影像。由于大范围影像拼接计算过程中，大量重采样计算主要发生在基于三维点云或三维数字地形模型的正射影像生成阶段，在顾及大范围农业环境，一般具有相对平整的地面，所以，在本发明中，将现有三维计算模型通过几种假设简化为二维模型，由此使整个拼接过程，不仅简化了计算，同时也避免了过多的重采样计算，确保光谱信息最少的损失。
因此，本发明以单波段替代所有波段的特征点提取与匹配，同时回避了三维点云的生成与正射影像的重采样，从而极大地简化并加快了计算处理过程；此外，将单波段计算的外方位元素用于所有波段的融合处理，确保了所有波段拼接结果的一致性；此外，采用简化模型极大地降低了重采样计算，从而确保影像在融合处理过程中光谱信息与原始影像尽可能一致。
验证例
为了验证本发明提供的多光谱影像快速拼接方法的可行性和优异的效果性，使用AF25b无人直升机搭载Tetracam Mini-MCA6多光谱相机对新疆北屯的两个不同区域进行数据采集，其中，该相机具有6个波段，每一波段具有1280*1024像素，水平视角为38.26度，得到数据集1和数据集2。其中，数据集1主要内容为城市、裸土和部分水体，共441张影像，航高为400m。数据集2主要内容为农田和公路，共有222张影像，航高为200m。
使用本发明提供的拼接方法对数据集1和数据集2分别进行拼接处理，结果如图4和图6所示，作为对比，使用软件Pix4UAV分别对数据集1和数据集2进行拼接处理，其中，对于数据集1，其在拼接过程中生成了三维点云，基于点云生成了正射影像，拼接结果如图5所示。而使用软件Pix4UAV对数据集2处理时，无法得到完整的拼接结果。
观察图4和图6可以看出，采用本发明拼接方法时，两组数据集所对应的拼接影像没有明显的大范围变形和错位，且图4与图5在视角上基本一致，说明本发明方法在两组数据集上都得到视觉上令人满意的效果。此外，由于本发明采用了较为简单的模型，本发明方法可以更有效的保留影像边角部分的数据，使得在最终结果中可以尽可能提高原始数据的利用率。在数据集1利用Pix4UAV拼接的结果中，观察图5可以看出，影像左下角的区域产生了轻微错位，可能是由于三维重建的不稳定特性造成。对于数据集2，利用Pix4UAV拼接时，处理记录显示，无法成功进行影像拼接的原因是：在平差过程中因为误差过大等原因丢弃了大量的影像，主要原因在于数据集2的影像几乎都是绿色农田，特征点不明显，难以得到良好的匹配结果。在这一方面，本发明拼接方法明显优于Pix4UAV。这也显示本发明拼接方法相对于传统的三维投影模型方法，具有更好的稳定性，在面对相对恶劣的原始数据时仍能稳定处理。
为了验证本发明拼接方法在处理速度方面的优势，通过增加影像数据的方式进行了测试，结果如下表1所示。
表1 数据集1的处理速度比较
影像数据(张) 本发明所用时间 Pix4UAV所用时间
147 421s 960s
254 998s 3720s
362 2047 11400s
表1的数据显示，随着影像数量的增长，本发明和已有商业软件需要的时间均有不同程度的增长。然而，已有商业软件的时间消耗增长剧烈，而本发明方法时间消耗基本为线性增长。因此，在数据量较多时，本发明方法可以显著节省处理时间。
为了验证本发明方法的光谱精度，采用了一种基于特征匹配的对比方法。首先对于影像中的单个点，通过比较这个点在每个波段上与原始影像上相同点的数值差异，从而评价这个点相比原始影像的光谱失真。而对于整幅影像，使用人工方式选点比较需要耗费大量的工作，而通过个别点配准的方式对齐影像无法做到逐像素精确对齐，所以也难以用于大面积的影像比较。所以，发明人挑选了多张原始影像，将原始影像与拼接结果进行特征匹配，然后对比匹配点之间的不同波段数值差异，这样就可以得到相对较多的采样点。由于存在原始影像与本发明方法的匹配结果/光束法平差生成的正射影像两组匹配，为了去除两组匹配点数不一样造成的结果差异，在对比时只取两组结果中对应原始影像上相同点的点对。测试区域如图7所示，共有7个采样区域，光谱失真对比如表2所示：
表2：光谱失真对比表
从表2可以看出，本发明方法在农田等植被覆盖区域表现明显优于商业软件基于光束法平差的三维重建方式。本发明在农田和植被为主的区域保持了低于基于光束法平差的正射影像的平均方均根误差，可以说明本发明方法确实可以实现更低的光谱失真，从而为进一步的遥感分析打下良好的基础。因此，本发明方法非常适合对面向农田等平坦区域进行植被分析。
综上所述，与传统的无人机影像拼接方法相比，本发明具有以下优点：
(1)充分考虑了面向农田等平坦地形场景的多光谱遥感应用的特点，降低对空间精度的要求，在地形起伏影响较小的情况下将地面视为平面，从而可以使用基于简化的投影模型实现简单的计算处理，并且避免了精确还原三维地形带来的大量不规则重采样，相对于现有方法，最大限度地保留了多光谱遥感影像的原始光谱信息。
(2)通过引用特征描述子组合，优选少量具有最佳特征信息的影像进行特征点提取与影像空间关系计算，并将这一计算结果用于其它波段影像的直接融合，从而避免了传统方法对所有波段的单独处理，以及传统方法分组处理的缺点，使得拼接后的大幅面影像仍然保留了所有波段，并可实现波段间的随意组合分析。
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。

Claims (4)

1.一种无人机多光谱影像快速拼接方法，其特征在于，包括以下步骤：
步骤1，获取无人机采集的待拼接的n张多光谱遥感影像，同时获取无人机采集所述n张多光谱遥感影像时的航迹数据；
根据所述航迹数据，分析得到各张多光谱遥感影像之间的相邻位置关系，从而确定所有存在潜在重叠区域的影像对，假设存在潜在重叠区域的影像对共有P组；
步骤2，确定最佳波段和最佳匹配算法，步骤如下：
步骤2.1，从P组存在潜在重叠区域的影像对中，挑选出R组存在潜在重叠区域的影像对，分别记为：第1组影像对、第2组影像对...第R组影像对；其中，R远小于P；
步骤2.2，确定需要筛选的波段集合和匹配算法集合，假设波段集合共有a个波段，分别为：第1波段、第2波段...第a波段；假设匹配算法集合共有b个匹配算法，分别为：第1匹配算法、第2匹配算法...第b匹配算法；
步骤2.3，使用不同波段和不同匹配算法的组合进行正交试验，即：使用每种波段-匹配算法组合对挑选出的每组影像对进行特征匹配，并对特征匹配结果进行记录，得到每种波段-匹配算法组合所对应的特征匹配描述子平均值最高值的特征匹配描述子平均值所对应的波段-匹配算法组合即为最佳的波段-匹配算法组合，由此得到最佳波段和最佳匹配算法；
其中，对于任意的第f波段和第g匹配算法，其中，f∈(1，2...a)；g∈(1，2...b)，其所对应的特征匹配描述子平均值基于以下公式计算得到：
<mrow> <mover> <mi>P</mi> <mo>&amp;OverBar;</mo> </mover> <mo>=</mo> <mrow> <mo>(</mo> <msubsup> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>R</mi> </msubsup> <msub> <mi>P</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>/</mo> <mi>R</mi> </mrow>
其中，Pi代表使用第f波段和第g匹配算法，对第i组影像对进行特征点提取与匹配计算时，得到的第i特征匹配描述子；i∈(1，2...R)；
其中：Pi＝Viri；
其中：Vi＝C(Ii-1,Ii-2)/t；Ii-1和Ii-2分别代表第i组影像对所包含的两个影像；C(Ii-1,Ii-2)代表第i组影像对成功匹配的点对数量；t代表匹配第i组影像对所消耗的时间；Vi代表单位数量匹配点对所耗费的时间；
ri＝C(Ii-1,Ii-2)/mean(Ii-1,Ii-2)；其中，mean(Ii-1,Ii-2)代表第i组影像对提取到的特征点数量的平均值，即：第Ii-1影像提取到的特征点数量和第Ii-2影像提取到的特征点数量的平均值；ri代表第i组影像对的匹配成功率；
步骤3，使用步骤2确定的最佳波段和最佳匹配算法，对剩余的P-R组存在潜在重叠区域的影像对进行特征匹配；
基于步骤2确定的最佳波段和最佳匹配算法，查找步骤2.3的特征匹配结果记录，得到与最佳波段和最佳匹配算法对应的R组存在潜在重叠区域的影像对的特征匹配结果；
由此得到全部P组影像对在最佳波段和最佳匹配算法下的特征匹配结果；
步骤4，对所述特征匹配结果进行分析，识别出所有的同名点；
步骤5，对于待拼接的n张多光谱遥感影像，分别记为：影像1、影像2...影像n；
步骤5.1，令i＝1；
步骤5.2，根据航迹数据，获得影像i对应的原始的相机投影中心和航向角κi的值；还获得对应的航高Hi的值，其中，航高Hi为固定值；由此得到每张影像Vi对应的λi的值；其中，λi为相机焦距fi与航高Hi的比值，为固定值；
步骤5.3，获得影像i的相邻且未处理过的相邻影像，假设相邻且未处理过的相邻影像的数量为wi，将wi个相邻影像依次记为：相邻影像1、相邻影像2...相邻影像wi；
步骤5.4，令v＝1；
步骤5.5，根据步骤4的结果，获得影像i和相邻影像v匹配得到的mv个同名点对的像点坐标，将mv个同名点对的像点坐标按顺序从1开始编号；
步骤5.6，令j＝1；
步骤5.7，按下式计算第j个同名点对所对应的两个地面三维点坐标之间的距离：
<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <mo>(</mo> <msup> <mrow> <mo>(</mo> <msubsup> <mi>X</mi> <mrow> <mi>v</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>p</mi> <mn>1</mn> </mrow> <mi>i</mi> </msubsup> <mo>-</mo> <msubsup> <mi>X</mi> <mrow> <mi>p</mi> <mn>2</mn> </mrow> <mrow> <mi>v</mi> <mo>,</mo> <mi>j</mi> </mrow> </msubsup> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msubsup> <mi>Y</mi> <mrow> <mi>v</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>p</mi> <mn>1</mn> </mrow> <mi>i</mi> </msubsup> <mo>-</mo> <msubsup> <mi>Y</mi> <mrow> <mi>p</mi> <mn>2</mn> </mrow> <mrow> <mi>v</mi> <mo>,</mo> <mi>j</mi> </mrow> </msubsup> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>=</mo> <mo>(</mo> <mo>(</mo> <msubsup> <mi>x</mi> <mrow> <mi>v</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mn>1</mn> </mrow> <mi>i</mi> </msubsup> <mi>cos</mi> <mi> </mi> <msub> <mi>k</mi> <mi>i</mi> </msub> <mo>+</mo> <msubsup> <mi>y</mi> <mrow> <mi>v</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mn>1</mn> </mrow> <mi>i</mi> </msubsup> <mi>sin</mi> <mi> </mi> <msub> <mi>k</mi> <mi>i</mi> </msub> <mo>-</mo> <msubsup> <mi>v</mi> <mrow> <mi>v</mi> <mo>,</mo> <mn>2</mn> </mrow> <mi>j</mi> </msubsup> <mi>cos</mi> <mi> </mi> <msub> <mi>k</mi> <mi>v</mi> </msub> <mo>-</mo> <msubsup> <mi>y</mi> <mrow> <mi>v</mi> <mo>,</mo> <mn>2</mn> </mrow> <mi>j</mi> </msubsup> <mi>sin</mi> <mi> </mi> <msub> <mi>k</mi> <mi>v</mi> </msub> <mo>+</mo> <msubsup> <mi>X</mi> <mi>c</mi> <mi>i</mi> </msubsup> <mo>-</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msubsup> <mi>X</mi> <mi>c</mi> <mi>v</mi> </msubsup> <msup> <mo>)</mo> <mn>2</mn> </msup> <mo>+</mo> <mrow> <mo>(</mo> <msup> <mrow> <mo>(</mo> <mrow> <msubsup> <mi>y</mi> <mrow> <mi>v</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mn>1</mn> </mrow> <mi>i</mi> </msubsup> <mi>cos</mi> <mi> </mi> <msub> <mi>k</mi> <mi>i</mi> </msub> <mo>-</mo> <msubsup> <mi>x</mi> <mrow> <mi>v</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mn>1</mn> </mrow> <mi>i</mi> </msubsup> <mi>sin</mi> <mi> </mi> <msub> <mi>k</mi> <mi>i</mi> </msub> <mo>-</mo> <msubsup> <mi>y</mi> <mrow> <mi>v</mi> <mo>,</mo> <mn>2</mn> </mrow> <mi>j</mi> </msubsup> <mi>cos</mi> <mi> </mi> <msub> <mi>k</mi> <mi>v</mi> </msub> <mo>+</mo> <msubsup> <mi>x</mi> <mrow> <mi>v</mi> <mo>,</mo> <mn>2</mn> </mrow> <mi>j</mi> </msubsup> <mi>sin</mi> <mi> </mi> <msub> <mi>k</mi> <mi>v</mi> </msub> <mo>+</mo> <msubsup> <mi>Y</mi> <mi>c</mi> <mi>i</mi> </msubsup> <mo>-</mo> <msubsup> <mi>Y</mi> <mi>c</mi> <mi>v</mi> </msubsup> </mrow> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>)</mo> </mrow> <mo>)</mo> <mo>;</mo> </mrow> </mtd> </mtr> </mtable> </mfenced>
其中：代表影像i上的像点所对应的地面二维坐标，并且，该像点与相邻影像v上的某个像点为同名点，并且，该同名点为影像i和相邻影像v之间的第j个同名点；
代表相邻影像v上的像点所对应的地面二维坐标，并且，该像点为与影像i上的某个像点为同名点，并且，该同名点为影像i和相邻影像v之间的第j个同名点；
代表影像i和相邻影像v之间第j个同名点对所对应的地面点对的距离；
表示影像i的像点坐标，并且，该像点与相邻影像v上的某个像点为同名点，并且，该同名点为影像i和相邻影像v之间的第j个同名点；
表示相邻影像v上与对应的同名点的像点坐标；
ki和kv分别为影像i与相邻影像v的航向角；
分别为影像i与相邻影像v的相机投影中心坐标；
步骤5.8，判断j是否等于mv，如果等于，则执行步骤5.9；否则，令j＝j+1，返回步骤5.7；
步骤5.9，判断v是否等于wi，如果等于，则执行步骤5.10；否则，令v＝v+1，返回步骤5.5；
步骤5.10，判断i是否等于n，如果等于，则执行步骤5.11；否则，令i＝i+1，返回步骤5.2；
步骤5.11，对步骤1到步骤5.10所计算得到的地面三维点坐标之间的距离进行求和计算，由此得到如下的目标函数：
<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <mi>m</mi> <mi>i</mi> <mi>n</mi> <msubsup> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </msubsup> <msubsup> <mo>&amp;Sigma;</mo> <mrow> <mi>v</mi> <mo>=</mo> <mn>1</mn> </mrow> <msup> <mi>w</mi> <mi>i</mi> </msup> </msubsup> <msubsup> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <msup> <mi>m</mi> <mi>v</mi> </msup> </msubsup> <mrow> <mo>(</mo> <mo>(</mo> <msubsup> <mi>x</mi> <mrow> <mi>v</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mn>1</mn> </mrow> <mi>i</mi> </msubsup> <mi>cos</mi> <mi> </mi> <msub> <mi>k</mi> <mi>i</mi> </msub> <mo>+</mo> <msubsup> <mi>y</mi> <mrow> <mi>v</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mn>1</mn> </mrow> <mi>i</mi> </msubsup> <mi>sin</mi> <mi> </mi> <msub> <mi>k</mi> <mi>i</mi> </msub> <mo>-</mo> <msubsup> <mi>x</mi> <mrow> <mi>v</mi> <mo>,</mo> <mn>2</mn> </mrow> <mi>j</mi> </msubsup> <mi>cos</mi> <mi> </mi> <msub> <mi>k</mi> <mi>v</mi> </msub> <mo>-</mo> <msubsup> <mi>y</mi> <mrow> <mi>v</mi> <mo>,</mo> <mn>2</mn> </mrow> <mi>j</mi> </msubsup> <mi>sin</mi> <mi> </mi> <msub> <mi>k</mi> <mi>v</mi> </msub> <mo>+</mo> <msubsup> <mi>X</mi> <mi>c</mi> <mi>i</mi> </msubsup> <mo>-</mo> <msubsup> <mi>X</mi> <mi>c</mi> <mi>v</mi> </msubsup> <mo>)</mo> </mrow> <mn>2</mn> <mo>+</mo> <mo>(</mo> <mo>(</mo> <msubsup> <mi>y</mi> <mrow> <mi>v</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mn>1</mn> </mrow> <mi>i</mi> </msubsup> <mi>cos</mi> <mi> </mi> <msub> <mi>k</mi> <mi>i</mi> </msub> <mo>-</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msubsup> <mi>x</mi> <mrow> <mi>v</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mn>1</mn> </mrow> <mi>i</mi> </msubsup> <mi>sin</mi> <mi> </mi> <msub> <mi>k</mi> <mi>i</mi> </msub> <mo>-</mo> <msubsup> <mi>y</mi> <mrow> <mi>v</mi> <mo>,</mo> <mn>2</mn> </mrow> <mi>j</mi> </msubsup> <mi>cos</mi> <mi> </mi> <msub> <mi>k</mi> <mi>v</mi> </msub> <mo>+</mo> <msubsup> <mi>x</mi> <mrow> <mi>v</mi> <mo>,</mo> <mn>2</mn> </mrow> <mi>j</mi> </msubsup> <mi>sin</mi> <mi> </mi> <msub> <mi>k</mi> <mi>v</mi> </msub> <mo>+</mo> <msubsup> <mi>Y</mi> <mi>c</mi> <mi>i</mi> </msubsup> <mo>-</mo> <msubsup> <mi>Y</mi> <mi>c</mi> <mi>v</mi> </msubsup> <msup> <mo>)</mo> <mn>2</mn> </msup> <mo>)</mo> <mo>)</mo> <mo>;</mo> </mrow> </mtd> </mtr> </mtable> </mfenced>
步骤6，以每张影像i对应的相机投影中心和航向角κi为优化参数，对步骤5.11的目标函数进行迭代求解，求得每张影像对应的最佳相机投影中心坐标和最佳航向角的值；
步骤7，利用步骤6得到的每张影像的最佳相机投影中心和最佳航向角，对每张影像的所有波段统一进行平移与旋转处理，使每张影像均平移旋转到最佳的位置与姿态；
步骤8，对经过步骤7处理后得到的影像进行融合处理，从而得到最终拼接形成的完整影像。
2.根据权利要求1所述的无人机多光谱影像快速拼接方法，其特征在于，步骤2.1中，R组存在潜在重叠区域的影像对为：顺序匹配得到的存在潜在重叠区域的影像对。
3.根据权利要求1所述的无人机多光谱影像快速拼接方法，其特征在于，步骤3中，剩余的P-R组存在潜在重叠区域的影像对是指：航带间匹配得到的存在潜在重叠区域的影像对。
4.根据权利要求1所述的无人机多光谱影像快速拼接方法，其特征在于，步骤8中，所述融合处理具体为：
当每张影像均平移旋转到最佳的位置与姿态后，对于任意两张存在重叠区域的影像，分别记为影像V1和影像V2；影像V1和影像V2分别具有G个波段，G个波段记为：第1波段、第2波段...第G波段；因此，影像V1和影像V2采用各个对应波段逐一进行融合的方法，即：影像V1的第1波段和影像V2的第1波段进行融合处理；影像V1的第2波段和影像V2的第2波段进行融合处理...依此类推，直到影像V1的第G波段和影像V2的第G波段进行融合处理；
其中，对于任意一个第h波段，其中，h∈(1、2...G)，影像V1的第h波段和影像V2的第h波段采用以下方法进行融合：
步骤8.1，对于影像V1的第h波段和影像V2的第h波段之间的重叠区域，在该重叠区域采样到多个采样点，任意一个采样点Oc的平面坐标为(x,y)；设影像V1的第h波段的中心位置点为OV1，影像V2的第h波段的中心位置点为OV2；
步骤8.2，如果采样点Oc到OV1的直线距离小于采样点Oc到OV2的直线距离，则：Or(x,y)点的像素值为：对影像V1的第h波段在(x,y)位置的像素点进行灰度采样，从而得到Or(x,y)点的像素值；其中，Or(x,y)点为：影像V1的第h波段和影像V2的第h波段进行拼接处理后得到影像V3，影像V3在(x,y)点的表达；
反之，如果采样点Oc到OV1的直线距离大于等于采样点Oc到OV2的直线距离，则：Or(x,y)点的像素值为：对影像V2的第h波段在(x,y)位置的像素点进行灰度采样，从而得到Or(x,y)点的像素值；
步骤8.3，因此，多个采样点Oc对应得到多个Or(x,y)点的像素值，多个Or(x,y)点的像素值即形成拼接影像V3在重叠区域的像素，由此实现影像融合。
CN201710552651.5A 2017-07-07 2017-07-07 一种无人机多光谱影像快速拼接方法 Active CN107274380B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710552651.5A CN107274380B (zh) 2017-07-07 2017-07-07 一种无人机多光谱影像快速拼接方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710552651.5A CN107274380B (zh) 2017-07-07 2017-07-07 一种无人机多光谱影像快速拼接方法

Publications (2)

Publication Number Publication Date
CN107274380A true CN107274380A (zh) 2017-10-20
CN107274380B CN107274380B (zh) 2019-10-11

Family

ID=60073228

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710552651.5A Active CN107274380B (zh) 2017-07-07 2017-07-07 一种无人机多光谱影像快速拼接方法

Country Status (1)

Country Link
CN (1) CN107274380B (zh)

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108734685A (zh) * 2018-05-10 2018-11-02 中国矿业大学(北京) 一种无人机载高光谱线阵列遥感影像的拼接方法
CN111754603A (zh) * 2020-06-23 2020-10-09 自然资源部四川测绘产品质量监督检验站(四川省测绘产品质量监督检验站) 一种无人机影像连接图构建方法及系统
CN112991186A (zh) * 2021-04-27 2021-06-18 湖南大学 一种无人机大视场高光谱图像生成方法及系统
CN113125351A (zh) * 2021-03-25 2021-07-16 国家卫星气象中心(国家空间天气监测预警中心) 一种多时次遥感影像优化合成方法及系统
CN109447902B (zh) * 2018-10-15 2023-04-25 广州地理研究所 一种图像拼接方法、装置、储存介质及设备
CN116977868A (zh) * 2023-06-07 2023-10-31 珠江水利委员会珠江水利科学研究院 一种基于特征匹配的影像乘积融合方法、系统及存储介质
US11964762B2 (en) 2020-02-11 2024-04-23 Raytheon Company Collaborative 3D mapping and surface registration

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101866482A (zh) * 2010-06-21 2010-10-20 清华大学 基于摄像设备自标定技术的全景图拼接方法和装置
US20120237137A1 (en) * 2008-12-15 2012-09-20 National Tsing Hua University (Taiwan) Optimal Multi-resolution Blending of Confocal Microscope Images
CN103516995A (zh) * 2012-06-19 2014-01-15 中南大学 一种基于orb特征的实时全景视频拼接方法和装置
CN105184736A (zh) * 2015-09-09 2015-12-23 山东大学 一种窄重叠双视场高光谱成像仪的图像配准的方法
CN105844587A (zh) * 2016-03-17 2016-08-10 河南理工大学 一种低空无人机载高光谱遥感影像自动拼接方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20120237137A1 (en) * 2008-12-15 2012-09-20 National Tsing Hua University (Taiwan) Optimal Multi-resolution Blending of Confocal Microscope Images
CN101866482A (zh) * 2010-06-21 2010-10-20 清华大学 基于摄像设备自标定技术的全景图拼接方法和装置
CN103516995A (zh) * 2012-06-19 2014-01-15 中南大学 一种基于orb特征的实时全景视频拼接方法和装置
CN105184736A (zh) * 2015-09-09 2015-12-23 山东大学 一种窄重叠双视场高光谱成像仪的图像配准的方法
CN105844587A (zh) * 2016-03-17 2016-08-10 河南理工大学 一种低空无人机载高光谱遥感影像自动拼接方法

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108734685A (zh) * 2018-05-10 2018-11-02 中国矿业大学(北京) 一种无人机载高光谱线阵列遥感影像的拼接方法
CN108734685B (zh) * 2018-05-10 2022-06-03 中国矿业大学(北京) 一种无人机载高光谱线阵列遥感影像的拼接方法
CN109447902B (zh) * 2018-10-15 2023-04-25 广州地理研究所 一种图像拼接方法、装置、储存介质及设备
US11964762B2 (en) 2020-02-11 2024-04-23 Raytheon Company Collaborative 3D mapping and surface registration
CN111754603A (zh) * 2020-06-23 2020-10-09 自然资源部四川测绘产品质量监督检验站(四川省测绘产品质量监督检验站) 一种无人机影像连接图构建方法及系统
CN111754603B (zh) * 2020-06-23 2024-02-13 自然资源部四川测绘产品质量监督检验站(四川省测绘产品质量监督检验站) 一种无人机影像连接图构建方法及系统
CN113125351A (zh) * 2021-03-25 2021-07-16 国家卫星气象中心(国家空间天气监测预警中心) 一种多时次遥感影像优化合成方法及系统
CN113125351B (zh) * 2021-03-25 2022-11-29 国家卫星气象中心(国家空间天气监测预警中心) 一种多时次遥感影像优化合成方法及系统
CN112991186A (zh) * 2021-04-27 2021-06-18 湖南大学 一种无人机大视场高光谱图像生成方法及系统
CN116977868A (zh) * 2023-06-07 2023-10-31 珠江水利委员会珠江水利科学研究院 一种基于特征匹配的影像乘积融合方法、系统及存储介质
CN116977868B (zh) * 2023-06-07 2024-03-01 珠江水利委员会珠江水利科学研究院 一种基于特征匹配的影像乘积融合方法、系统及存储介质

Also Published As

Publication number Publication date
CN107274380B (zh) 2019-10-11

Similar Documents

Publication Publication Date Title
CN107274380B (zh) 一种无人机多光谱影像快速拼接方法
CN104156968B (zh) 大面积地形复杂区域无人机序列影像快速无缝拼接方法
US9798928B2 (en) System for collecting and processing aerial imagery with enhanced 3D and NIR imaging capability
CN103822616B (zh) 一种图分割与地形起伏约束相结合的遥感影像匹配方法
CN106127697B (zh) 无人机机载成像高光谱几何校正方法
CN106373088B (zh) 大倾斜低重叠率航空图像的快速拼接方法
CN107918927A (zh) 一种匹配策略融合及低误差的快速图像拼接方法
CN104732482A (zh) 一种基于控制点的多分辨率图像拼接方法
CN101930603B (zh) 中高速传感器网络图像数据融合的方法
CN106485751B (zh) 应用于基桩检测中的无人机摄影成像及数据处理方法及系统
EP3193136B1 (en) Method and system for geometric referencing of multi-spectral data
Guo et al. Mapping crop status from an unmanned aerial vehicle for precision agriculture applications
Li et al. A study on automatic UAV image mosaic method for paroxysmal disaster
CN107192376A (zh) 基于帧间连续性的无人机多帧图像目标定位校正方法
CN105931185A (zh) 一种多视角图像自动拼接方法
CN104318583A (zh) 一种可见光宽带光谱图像配准方法
Ribera et al. Estimating phenotypic traits from UAV based RGB imagery
CN110097498A (zh) 基于无人机航迹约束的多航带图像拼接与定位方法
CN110084743A (zh) 基于多航带起始航迹约束的图像拼接与定位方法
CN112862683A (zh) 一种基于弹性配准和网格优化的邻接图像拼接方法
CN113610905A (zh) 基于子图像匹配的深度学习遥感图像配准方法及应用
CN114581307A (zh) 用于目标追踪识别的多图像拼接方法、系统、设备及介质
CN116228539A (zh) 一种无人机遥感图像拼接的方法
CN106204507A (zh) 一种无人机图像拼接方法
CN117173224A (zh) 基于Delauny算法和Kruskal最大生成树算法的大田影像拼接方法

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
GR01 Patent grant
GR01 Patent grant