CN116758223A - 基于双角度多光谱图像的三维多光谱点云重建方法、系统及设备 - Google Patents
基于双角度多光谱图像的三维多光谱点云重建方法、系统及设备 Download PDFInfo
- Publication number
- CN116758223A CN116758223A CN202310823147.XA CN202310823147A CN116758223A CN 116758223 A CN116758223 A CN 116758223A CN 202310823147 A CN202310823147 A CN 202310823147A CN 116758223 A CN116758223 A CN 116758223A
- Authority
- CN
- China
- Prior art keywords
- image
- multispectral
- point
- matching
- point cloud
- 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 96
- 238000012937 correction Methods 0.000 claims abstract description 40
- 230000002159 abnormal effect Effects 0.000 claims abstract description 23
- 238000000605 extraction Methods 0.000 claims abstract description 18
- 238000001228 spectrum Methods 0.000 claims abstract description 17
- 238000001914 filtration Methods 0.000 claims abstract description 12
- 230000001131 transforming effect Effects 0.000 claims abstract description 10
- 239000011159 matrix material Substances 0.000 claims description 39
- 238000010586 diagram Methods 0.000 claims description 32
- 238000002310 reflectometry Methods 0.000 claims description 32
- 238000004422 calculation algorithm Methods 0.000 claims description 28
- 230000008569 process Effects 0.000 claims description 24
- 230000003595 spectral effect Effects 0.000 claims description 19
- 230000009466 transformation Effects 0.000 claims description 18
- 238000012545 processing Methods 0.000 claims description 10
- 238000004364 calculation method Methods 0.000 claims description 9
- 238000013519 translation Methods 0.000 claims description 8
- 239000002131 composite material Substances 0.000 claims description 6
- 230000003287 optical effect Effects 0.000 claims description 4
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 3
- YIXLDRQOKWESBI-UHFFFAOYSA-N 4-hydroxy-6-pentylpyran-2-one Chemical compound CCCCCC1=CC(O)=CC(=O)O1 YIXLDRQOKWESBI-UHFFFAOYSA-N 0.000 claims description 3
- 230000000903 blocking effect Effects 0.000 claims description 3
- 238000012216 screening Methods 0.000 claims description 3
- 238000007619 statistical method Methods 0.000 claims description 3
- 230000002194 synthesizing effect Effects 0.000 claims description 2
- 230000000694 effects Effects 0.000 description 9
- 238000005516 engineering process Methods 0.000 description 7
- 230000008859 change Effects 0.000 description 6
- 238000004458 analytical method Methods 0.000 description 5
- 238000012800 visualization Methods 0.000 description 4
- 238000012544 monitoring process Methods 0.000 description 3
- 238000005070 sampling Methods 0.000 description 3
- 230000007547 defect Effects 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 230000018109 developmental process Effects 0.000 description 2
- 230000002708 enhancing effect Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 239000002028 Biomass Substances 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000000701 chemical imaging Methods 0.000 description 1
- 229930002875 chlorophyll Natural products 0.000 description 1
- 235000019804 chlorophyll Nutrition 0.000 description 1
- ATNHDLDRLWWWCB-AENOIHSZSA-M chlorophyll a Chemical compound C1([C@@H](C(=O)OC)C(=O)C2=C3C)=C2N2C3=CC(C(CC)=C3C)=[N+]4C3=CC3=C(C=C)C(C)=C5N3[Mg-2]42[N+]2=C1[C@@H](CCC(=O)OC\C=C(/C)CCC[C@H](C)CCC[C@H](C)CCCC(C)C)[C@H](C)C2=C5 ATNHDLDRLWWWCB-AENOIHSZSA-M 0.000 description 1
- 238000004040 coloring Methods 0.000 description 1
- 238000004590 computer program Methods 0.000 description 1
- 230000007123 defense Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 230000006870 function Effects 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 230000011218 segmentation Effects 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
-
- 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
-
- 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
-
- 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/80—Fusion, i.e. combining data from various sources at the sensor level, preprocessing level, feature extraction level or classification level
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- General Physics & Mathematics (AREA)
- Multimedia (AREA)
- Software Systems (AREA)
- Health & Medical Sciences (AREA)
- Artificial Intelligence (AREA)
- Computing Systems (AREA)
- Databases & Information Systems (AREA)
- Evolutionary Computation (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Geometry (AREA)
- Computer Graphics (AREA)
- Image Processing (AREA)
Abstract
基于双角度多光谱图像的三维多光谱点云重建方法、系统及设备，属于多光谱图像三维重建技术领域。为了解决现有三维重建方法不能直接从非配准相机中重建完整波段点云的问题。本发明首先针对设置角度不同的相机采集的图像进行纹理增强；相机为多光谱相机；针对倾斜图像进行倾斜图像投影校正和提取特征点，再将特征点坐标反变换回到倾斜图像中，从而确定倾斜图像中的特征点，基于特征点对原始图像进行匹配像对搜索和图像匹配；然后对图像执行多视图密集匹配得到深度图点云并得到多光谱深度图点云，将同一个相机所采集多光谱深度图点云拼接，再将不同的多光谱相机对应的包含不同波段光谱的点云进行拼接得到全波段多光谱点云，最后滤除异常点。
Description
技术领域
本发明属于多光谱图像三维重建技术领域，具体涉及一种基于双角度多光谱图像的三维多光谱点云重建方法、系统、存储介质及设备。
背景技术
随着传感器技术、无人机技术的发展，无人机遥感凭借其高机动性、高实时性、避免云层遮挡等优点，在军事防御、城市管理、国土资源环境监测、农业信息化等诸多领域获得了广泛应用。无人机平台可以根据具体任务需求搭配不同的传感器，其中多光谱相机能够采集具有高空间分辨率的多波段光谱图像，同时具有低成本、高效率等诸多优点，在农业、生态、城市监测等应用中发挥着重要应用价值。
受限于图像的二维空间维度，基于多光谱图像的遥感观测手段无法准确、完整地描述观测场景的空间几何特征，具体表现在无法感知地物的真实大小、高度等，导致许多定量遥感分析应用的精度无法保证。三维重建是一种利用多张重叠图像恢复观测场景空间三维信息的技术，能够消除图像中的几何畸变、地物遮挡情况，从空间三维视角对观测场景进行分析。点云是三维重建算法的主要目标数据，它是大量数据点的集合，点云中的每个点具有空间三维信息的同时能够加入多维光谱信息，能够良好弥补图像数据的缺陷。
垂直对地观测是当前无人机多光谱图像的主要采集方式，对观测场景中地物侧立面影像获取能力不足，导致三维重建结果完整性较低。随着传感器技术的多样化发展，搭载多台倾斜相机、从不同角度对地拍摄能够减轻地物侧立面影像获取能力不足的问题，在三维重建任务中发挥着显著的作用。
基于上述问题和研究背景，本发明结合多光谱成像技术、三维重建技术与倾斜摄影技术优势，使用双角度多光谱相机采集两组5波段多光谱图像，针对双角度多光谱图像重建三维多光谱点云，对于光谱指数分析、三维场景建模、多维立体探测等定性、定量遥感都具有重要的研究意义和应用价值。
发明内容
本发明为了解决现有三维重建方法不能直接从非配准相机中重建完整波段点云的问题。
基于双角度多光谱图像的三维多光谱点云重建方法，包括以下步骤：
S1、针对设置角度不同的相机采集的图像进行纹理增强；所述相机为多光谱相机，且设置角度不同的相机对应的光谱波段不同，多光谱相机采集的图像为多光谱反射率图；
S2、基于每个多光谱图像纹理增强后的图像，针对倾斜图像，进行倾斜图像投影校正，使用尺度不变特征变换匹配算法，SIFT算法，对投影校正后的图像提取特征点，再将特征点坐标反变换回到倾斜图像中，从而确定倾斜图像中的特征点，也就是原始图像中的特征点；
S3、基于原始图像中的特征点，对原始图像进行匹配像对搜索，搜索与参考图像重叠率超过重叠率阈值的图像用于后续的匹配；
S4、针对S3搜索到的图像进行图像匹配；
S5、基于S4确定的匹配点对估计每张图像外参数、相机的内参数，对图像执行多视图密集匹配得到深度图点云；
在获得每张图像的深度图点云后，结合相机内参数和图像外参数，对多光谱反射率图去畸变并从中获得深度图点云每个点的多波段光谱值，从而得到多光谱深度图点云；
然后将同一个相机所采集多光谱深度图点云拼接得到相同波段的多光谱深度图点云，再基于空间位置一致性将设置角度不同的相机对应的不同波段光谱的点云进行拼接，生成完整的全波段多光谱点云，最后滤除异常点；
所述的全波段多光谱点云为覆盖不同相机所有波段的多光谱点云。
进一步地，针对不同的多光谱相机采集的图像进行纹理增强包括以下步骤：
首先利用不同的多光谱相机采集图像，然后利用每个相机计算归一化差分植被指数图NDVI；
定义运算规则A＝[B⊙λb]表示对于矩阵B中的每个元素，若其值⊙λb，则与B相同维度的矩阵A中对应位置元素值为1，否则为0；λb表示一个数值，⊙表示＞、≥、<、≤符号中的一种；
利用下式去除植被指数图的阴影区域，
INDVI用于表示不同多光谱相机的植被指数图NDVI，bj用于表示不同多光谱相机的中心波长为j的单波段反射率图；
然后将植被指数图中的植被像素与蓝、绿、红三波段反射率图中的非植被像素进行合成，将植被指数值小于0.3的像素视为非植被像素，计算公式如下：
式中，分别表示多光谱图像中对应红、绿、蓝波段的反射率图像；图像I为最终获得的具有强几何纹理的NDVI-RGB合成图。
进一步地，进行倾斜图像投影校正的过程包括以下步骤：
首先确定倾斜图像校正到下视视角图像对应的透视投影变换关系，将p(x,y,1)、P′(X,Y,Z)分别表示透视投影前的图像坐标、投影后的空间坐标，有
其中，M是图像透视投影矩阵，m11至m33为图像透视投影矩阵中的元素；
设p′(x′,y′)是投影后的图像坐标，由式(3)可建立如式(4)所示的方程，
令m33＝1，则方程中未知变量变为8个，需要4组投影前后的图像坐标即可求解矩阵M；使用倾斜图像的4个顶点坐标作为p，求解对应4个投影点p′：
首先将倾斜图像中心点的像素齐次坐标o(u0,v0,1)后向投影到水平地面上，设其投影点的世界坐标为根据位姿数据可求得倾斜图像在世界坐标系中的旋转矩阵R1；
然后将每个相机曝光位置的经纬高坐标转换为地心固定直角坐标，获得图像在世界坐标系中的平移向量t1；
基于图像在世界坐标系中的旋转矩阵R1，根据式(7)求得图像中心点的地面投影点：
式中，K1为倾斜图像的相机内参数；D代表相机光心到地面的高度，即投影点的深度；设初始深度D为无人机飞行高度，当时更新深度/>并利用式(7)重新计算，α1为世界坐标Z轴基准阈值，直到/>时即求得地面投影点坐标/>同时获得其深度D；
将倾斜图像的顶点坐标投影到地面得到的地面点坐标
然后设校正图像的滚转角、俯仰角均为0，航偏角与倾斜图像保持相同，即(α,β,γ)＝(0,0,γcorr)，并求得校正图像在世界坐标系下的旋转矩阵为R2；
记校正图像在世界坐标为oC(XW,YW,h)，令h＝D，得到校正图像的平移向量t2；
最后，利用式(9)将4个地面点前向投影到校正图像上即可得到对应的4个图像投影点；
p′＝K1[R2pg+t2] (9)
将该4组投影前后的图像坐标代入式(4)中可得到8个方程，求解方程获得图像投影变换矩阵M并应用式(3)对倾斜图像进行变换，将变换后的空间坐标P′(X,Y,Z)归一化为图像坐标，即至此完成了倾斜图像的投影校正。
进一步地，基于原始图像中的特征点，对原始图像进行匹配像对搜索的过程包括以下步骤：
首先基于每张原始图像的位姿数据求得各自的地心固定直角坐标同时确定原始图像该数据集中最小的经度lonref、纬度latref，令参考高度为altref＝0的空间点为参考点，并求得参考点的地心固定直角坐标为/>由式(10)得到每张原始图像在世界坐标系下的坐标，同时分别求解世界坐标系中的旋转矩阵R和平移向量t；
然后将每张图像的中心点后向投影到地面获得投影点OG(XG,YG,0)；
设某一个相机在OC位置拍摄的一张后视倾斜图像的中心点的地面投影点坐标为OG，则参考图像的邻居匹配图像投影点O′G满足式(11)：
||OG-O′G||≤rs (11)
其中，rs为搜索半径，根据图像重叠率估计得到；
搜索与参考图像重叠率超过重叠率阈值的图像用于后续的匹配。
进一步地，针对S3搜索到的图像进行图像匹配的过程包括以下步骤：
首先，使用KD-Tree算法对参考图像中的特征点在匹配图像上搜索最近邻和次近邻，设最近邻、次近邻与参考图像中的特征点的欧式距离分别为d1、d2，当二者满足式(13)时参考图像中的特征点与最近邻匹配，即，同时将r作为该匹配点对的置信度；
其中，η是匹配点筛选阈值；
针对匹配结果，使用MAGSAC++算法剔除误匹配点对，在剔除误匹配点对的过程中，采用最大化一致性原则估计单应变换矩阵；
在估计得到进行匹配的两张图像之间的单应变换矩阵Md之后，将匹配图像的特征点坐标投影到参考图像坐标系中，并对参考图像分块，每一个块为一个子图；然后遍历参考图像的每一个子图，进行局部重匹配，即对子图区域内的参考特征点与匹配特征点执行KD-Tree近邻搜索匹配，然后合并所有子图的匹配结果；
交换参考图像与匹配图像，同时利用Md -1对新的匹配图像特征点投影，对新的参考图像分块并执行局部重匹配；
取双向匹配结果的交集作为两图像的匹配结果；
最后针对匹配结果采用MAGSAC++算法去除误匹配点对，得到最终的匹配点对。
进一步地，S5中基于S4确定的匹配点对估计每张图像外参数、相机的内参数的过程是通过运动恢复结构原理实现的。
进一步地，S5所述最后滤除异常点的过程中，使用开源点云数据处理库PDAL中基于统计学的方法进行异常点识别：
对多光谱点云中的每个点搜索k个近邻点并计算距离均值和标准差σ，得到异常点判断阈值：
其中，m为倍乘系数；
当i近邻点的距离ui≥t时将被标记为异常点并滤除。
基于双角度多光谱图像的三维多光谱点云重建系统，包括图像纹理增强单元、投影校正及特征点提取单元、匹配像对搜索单元、图像匹配单元和全波段多光谱点云生成单元；其中，
图像纹理增强单元：针对不同的多光谱相机采集的图像进行纹理增强；所述不同的多光谱相机为光谱波段不同的相机，多光谱相机采集的图像为多光谱反射率图；
投影校正及特征点提取单元：基于每个多光谱图像纹理增强后的图像，针对倾斜图像，进行倾斜图像投影校正，使用SIFT算法对投影校正后的图像提取特征点，再将特征点坐标反变换回到倾斜图像中，从而确定倾斜图像中的特征点，也就是原始图像中的特征点；
匹配像对搜索单元：基于原始图像中的特征点，对原始图像进行匹配像对搜索，搜索与参考图像重叠率超过重叠率阈值的邻居图像用于后续的匹配；
图像匹配单元：针对搜索到的图像进行图像匹配；
全波段多光谱点云生成单元：首先基于确定的匹配点对估计每张图像外参数、相机的内参数，对图像执行多视图密集匹配得到深度图点云；
在获得每张图像的深度图点云后，结合相机内参数和图像外参数，对多光谱反射率图去畸变并从中获得深度图点云每个点的多波段光谱值，从而得到多光谱深度图点云；
然后将同一个相机所采集多光谱深度图点云拼接得到相同波段的多光谱深度图点云，在将不同的多光谱相机对应的包含不同波段光谱的点云进行拼接，基于空间位置一致性生成完整的全波段多光谱点云，最后滤除异常点；
所述的全波段多光谱点云为覆盖不同的多光谱相机所有波段的多光谱点云。
一种计算机存储介质，所述存储介质中存储有至少一条指令，所述至少一条指令由处理器加载并执行以实现所述的一种基于双角度多光谱图像的三维多光谱点云重建方法。
一种基于双角度多光谱图像的三维多光谱点云重建设备，所述设备包括处理器和存储器，所述存储器中存储有至少一条指令，所述至少一条指令由处理器加载并执行以实现所述的一种基于双角度多光谱图像的三维多光谱点云重建方法。
有益效果：
1、本发明中的多光谱图像纹理增强方法能够提高图像中植被像素与非植被像素的对比度，增强图像灰度纹理、减轻多波段灰度非线性变化现象，倾斜图像投影校正方法能够将大倾角图像校正到下视视角，显著减小透视几何畸变程度，匹配像对搜索与图像局部重匹配方法能够显著减少误匹配、提高正确匹配数量。以上处理算法大幅提高了倾斜多光谱图像特征点提取与匹配对数，保证后续三维重建算法的精度和完整性。
2、本发明采用倾斜双角度视角采集多光谱图像，相比于单角度下视视角，能够获取更多地物侧立面影像，结合投影校正与局部重匹配方法，能够获得更多地物侧立面重建点云，最终重建场景的空间结构具有较高完整性。
3、由于不同视角相机拍摄图像无法直接配准，现有三维重建方法不能直接从非配准相机中重建完整波段点云。本发明所设计的基于深度图的多光谱点云生成方法，首先估计每张图像的深度图点云并将相同相机的点云融合，再从全部相机点云中提取最近邻点得到每个空间点的完整光谱波段反射率值，实现完整波段多光谱点云生成。
附图说明
图1为基于双角度多光谱图像的三维多光谱点云重建方法基本流程图。
图2为不同地物的部分波段反射率可能相同的示意图。
图3为投影校正与特征提取主要流程示意图。
图4为投影校正基本原理示意图。
图5为匹配像对搜索基本原理示意图。
图6(a)为特征点匹配的误匹配情况示意图；图6(b)为特征点匹配的漏匹配情况示意图。
图7为特征点匹配流程图。
图8为局部特征点重匹配示意图。
图9为运动恢复结构算法示意图。
图10为多光谱点云生成流程示意图。
图11为基于空间位置一致性生成10波段点云。
图12为部分采集图像。
图13为纹理增强结果及局部放大对比示意图。
图14为图像投影校正前后对比图。
图15为匹配像对搜索结果可视化效果图。
图16为一组匹配结果可视化效果图。
图17为多光谱密集点云可视化与植被指数计算结果。
图18为单个地物点云提取结果。
具体实施方式
本发明基于位姿信息辅助的双角度多光谱图像特征提取与匹配，然后基于深度图的多光谱点云重建，能够针对大场景、复杂地物类型双角度多光谱图像，重建出光谱波段完整、空间结构精细的三维点云模型。
具体实施方式一：结合图1说明本实施方式，
本实施方式为基于双角度多光谱图像的三维多光谱点云重建方法，首先对本发明所使用的实验数据进行说明，主要信息见表1，位姿传感器的数据保存在每张图像的元数据区，飞行高度为80米，航向和航带重叠率均为80％。
表1实验数据主要信息
本实施方式所诉的基于双角度多光谱图像的三维多光谱点云重建方法，包括以下步骤：
S1、多光谱图像纹理增强：
多光谱图虽然能够反映地物的光谱特性，但由于地物对不同光谱波段的反射率不同，同一地物在不同波段图像上的像素亮度呈现非线性变化，且单波段反射率图的纹理较弱，导致基于梯度的特征点提取方法效果急剧下降。对于含有建筑、车辆、水体、树木、裸地等多种类型地物的场景，采用单波段反射率图的处理方式对波段的选取和地物类型变化非常敏感，例如，两个不同的地物点在某一波段的反射率值相同，导致在特征提取过程中无法区分，进而影响匹配和重建精度，如图2所示。
针对上述问题，本发明基于地物的光谱特性对图像纹理进行增强，首先利用两个相机采集图像，分别记为相机1和相机2，本实施方式中相机1和相机2为两个倾斜相机(针对被测区域，相机为倾斜设置，即倾斜相机)，然后利用每个相机的近红外或红边波段、红波段图计算归一化差分植被指数图(NDVI)，计算公式如下：
式中，bi是中心波长为i的单波段反射率图，分别是相机1和相机2图像计算得到的植被指数图，由于本发明使用的相机2没有近红外波段，因此选用中心波长相对接近的红边740波段图像代替。
定义运算规则A＝[B＞b]表示对于矩阵B中的每个元素，若其值大于b，则与B相同维度的矩阵A中对应位置元素值为1，否则为0；对于类似B＞b的其他比较关系，如大于等于、小于、小于等于，采用相同的方式定。
利用下式去除植被指数图的阴影区域，根据处理的相机不同，INDVI、bj分别为b842或/>b740。
Icorr＝INDVI×Mshd
Mshd＝[bj≥0.2]
由于图像中的植被通常具有大量重复纹理且在不同时刻其位置易发生变化，不利于图像匹配，而人造地物具有显著的点特征且位置固定不变，是图像匹配的主要依据。然而在植被指数图中非植被像素的指数值较小，图像纹理不明显，为了进一步增强几何纹理，保留人造地物的原始纹理特征，本发明将植被指数图中的植被像素与蓝、绿、红三波段反射率图中的非植被像素进行合成，将植被指数值小于0.3的像素视为非植被像素，计算公式如下：
式中分别表示多光谱图像中对应红、绿、蓝波段的反射率图像，图像I为最终获得的具有强几何纹理的NDVI-RGB合成图，后续算法将基于该合成图进行特征点提取与匹配。
S2、基于每个多光谱图像纹理增强后的图像，针对倾斜图像，进行倾斜图像投影校正：
虽然相机1和相机2为两个倾斜多光谱相机，但数据采集过程中易受到无人机飞行姿态的影响，导致采集的图像可能为无倾斜的图像；即使相机1或相机2不为倾斜相机，由于无人机飞行姿态的影响，极有可能采集到倾斜的图像；倾斜图像投影校正的过程是针对倾斜图像的，如果不是倾斜图像，是否进行投影校正均可。
如图3所示，进行投影校正与特征提取，将倾斜图像校正到下视视角图像的过程是一种透视投影变换，如式(3)所示，式中p(x,y,1)、P′(X,Y,Z)分别表示透视投影前的图像坐标、投影后的空间坐标，M是图像透视投影矩阵。
其中，m11至m33为图像透视投影矩阵中的元素。
设p′(x′,y′)是投影后的图像坐标(P′是三维空间点坐标，p′是对应的二维图像点坐标)，由式(3)可建立如式(4)所示的方程，令m33＝1，则方程中未知变量变为8个，因此只需要4组投影前后的图像坐标即可求解矩阵M。
因此使用倾斜图像的4个顶点坐标作为p，求解对应4个投影点p′：
首先根据相机位姿数据和传感器参数将倾斜图像中心点的像素齐次坐标o(u0,v0,1)后向投影到水平地面上，设其投影点的世界坐标为如图4所示。
根据位姿数据利用式(5)可求得倾斜图像在场景世界坐标系中的旋转矩阵R1：
其中，α、β、γ分别为滚转角、俯仰角、航偏角，RZ(γ)、RY(β)、RX(α)分别为绕Z、Y、X轴的二维旋转矩阵；
利用式(6)将每个相机曝光位置的经纬高坐标(lon,lat,alt)转换为地心固定直角坐标，获得该图像在该场景世界坐标系中的平移向量t1(X,Y,Z)：
式中，N为地球的曲率半径，a＝6378137.0为地球基准椭球体的长半径，f＝1/298.257223565为地球基准椭球体的极扁率。
基于图像在场景世界坐标系中的旋转矩阵R1，根据式(7)求得图像中心点的地面投影点：
式中，D代表相机光心到地面的高度，即投影点的深度。设初始深度D为无人机飞行高度，当时更新深度/>并利用式(7)重新计算，直到/>时即求得地面投影点坐标/>同时获得其深度D；K1为倾斜图像的相机内参数。
利用相同的迭代过程将倾斜图像的顶点坐标投影到地面，得到的地面点坐标记为
然后设校正图像的滚转角、俯仰角均为0，航偏角与倾斜图像保持相同，即(α,β,γ)＝(0,0,γcorr)，根据式(5)求得校正图像在该世界坐标系下的旋转矩阵为R2；令校正图像在该场景中的世界坐标为oC(XW,YW,h)，为保持投影前后图像中心区域的地面分辨率一致，令h＝D，因此校正图像的平移向量t2可由下式获得。
最后，利用式(9)将4个地面点前向投影到校正图像上即可得到对应的4个图像投影点。
p′＝K1[R2pg+t2] (9)
将该4组投影前后的图像坐标代入式(4)中可得到8个方程，求解方程获得图像投影变换矩阵M并应用式(3)对倾斜图像进行变换，将变换后的空间坐标P′(X,Y,Z)归一化为图像坐标，即至此完成了倾斜图像的投影校正。
使用SIFT(Scale Invariant Feature Transform)算法对投影校正后的图像提取特征点，再利用式(3)将特征点坐标反变换回到倾斜图像(也就是原始图像)中，确定倾斜图像中的特征点，也就是原始图像中的特征点，完成对倾斜图像的投影变换和特征点提取。
S3、匹配像对搜索：
基于原始图像中的特征点，场景中同一空间点在不同图像上检测到的特征点称为同名点，在提取到这些同名点之后需要建立对应关系以便估计图像的相对位姿，该过程即为特征点匹配。由于大场景三维重建任务中图像数量众多，所有图像全部匹配的效率极低，且由于每张图像覆盖地面范围较小，仅有少数图像之间存在重叠区域，因此本发明还利用相机初始位置姿态数据设计了一种匹配像对搜索方法，大幅减少图像匹配次数。
匹配像对搜索的基本原理如图5所示，图中大圆点(实际为灰色圆点)代表双角度多光谱相机系统曝光时刻的位置，小圆点(实际为红色、蓝色圆点)分别代表相机1、相机2所拍摄图像中心点的地面投影点位置。
首先利用每张原始图像的位姿数据和式(6)求得各自的地心固定直角坐标同时确定原始图像该数据集中最小的经度lonref、纬度latref，令参考高度为altref＝0的空间点为参考点，利用式(6)求得该场景的参考点坐标为/>由式(10)得到每张原始图像在世界坐标系下的坐标，同时用式(5)和式(8)分别求解旋转矩阵R和平移向量t。
然后利用式(7)的方法将每张图像的中心点后向投影到地面获得投影点OG(XG,YG,0)。
设相机1在OC位置拍摄的一张后视倾斜图像的中心点的地面投影点坐标为OG，则参考图像的邻居匹配图像投影点O′G满足式(11)：
||OG-O′G||≤rs (11)
其中，rs为搜索半径，根据图像重叠率估计得到，一般保证最小航向、航带重叠率均不低于50％即可。
利用该搜索方法对参考图像直接搜索重叠率在50％以上的邻居图像用于后续的匹配，匹配效率显著优于穷尽搜索匹配策略。
S4、图像局部重匹配：
为了进一步提高特征点匹配精度，利用初始匹配结果估计像对之间的单应变换关系并将待匹配图像特征点投影到参考图像上，在此基础上对参考图像分块并在分块内部进行局部重匹配。受透视投影几何畸变、光照变化辐射畸变等因素的影响，前面的处理能够在一定程度上改善这些问题，但不能彻底解决，因此提取到的图像特征点并不理想，所以利用图像局部重匹配进一步提高效果，当然也可以不进行此步骤。
在匹配过程中普遍存在如图6(a)和图6(b)所示的两种误匹配和漏匹配情况，导致最终匹配效果不佳。在特征点描述子所表示的高维空间中，当一个非匹配点是参考点的最近邻且匹配点与参考点距离过远时，将产生如图6(a)所示的误匹配情况。当两图像之间仅存在少量误匹配对时，可通过误匹配剔除算法去除，而当存在大量误匹配时，将可能导致后续的相机位姿估计发生错误，进而影响三维重建精度；当参考点的次近邻与最近邻匹配点过于接近时，由于无法有效区分两点，参考点将不会与任何点匹配，从而导致如图6(b)所示的漏匹配情况。以上两种错误情况均会导致原本能够正确匹配的点无法匹配，使得匹配点对数减少，进而导致相机位姿估计稳健性降低。为改善上述问题，本发明还提出了一种基于单应变换引导的局部重匹配方法，对于后向投影搜索到的一对待匹配图像，其特征点匹配流程如图7所示。
首先，使用快速近似最近邻搜索库FLANN(Fast Library for ApproximateNearest Neighbors)中的KD-Tree算法对参考图像中的特征点在匹配图像上搜索最近邻和次近邻，搜索依据为特征点描述子之间的欧式距离，如式(12)所示，di、dj分别为两张图像上的特征点描述子向量，dEuclid为向量欧式距离。
dEuclid＝||di-dj||2 (12)
设最近邻、次近邻与参考点的欧式距离分别为d1、d2，当二者满足式(13)时参考点与最近邻匹配，同时将r作为该匹配点对的置信度，r值越小表示该匹配点对的准确度越高，
其中，η是匹配点筛选阈值，一般根据匹配对数量需求进行调整，取值范围为0-1，值越小匹配精度越高但匹配数量更少。
针对上述初始匹配过程，初始匹配结果中一般存在大量误匹配点对，使用MAGSAC++算法(改进的边缘样本共识算法)剔除误匹配点对，它通过自适应几何采样选择数据点，采用最大化一致性原则估计单应变换矩阵，能够在具有大量噪声和外点的数据中实现更加稳健的估计。
在估计得到进行匹配的两张图像之间的单应变换矩阵Md之后，参考式(3)、(4)的形式将匹配图像的特征点坐标投影到参考图像坐标系中(图像的特征点是图像的像素点，以像素坐标表示，倾斜图像校正的过程也是像素点坐标的变换，二者本质相同，实现过程是根据前面说明的单应变换矩阵并参考公式(3)(4)实现)，并对参考图像分块，分块方式为按照虚线对参考图像长、宽均等分为5×5个子图，如图8所示。然后遍历参考图像的每一个子图，对子图区域内的参考特征点与匹配特征点执行KD-Tree近邻搜索匹配，然后合并所有子图的匹配结果；
为了减小投影误差的影响，交换参考图像与匹配图像(即将参考图像作为新的匹配图像，将匹配图像作为新的参考图像)，同时利用逆变换矩阵Md -1对新的匹配图像特征点投影，对新的参考图像分块并执行局部重匹配，取双向匹配结果的交集作为两图像的匹配结果。
最后，对匹配结果执行MAGSAC++算法去除误匹配点对，得到最终的匹配点对。
S5、多光谱点云生成：
运动恢复结构(Structure from Motion,SfM)是一种从多张重叠图像中估计相机参数和位姿信息并重建出场景三维结构的算法，其核心思想是利用不同视角的图像进行特征点匹配，根据多个匹配像对的空间几何约束估计优化相机参数和位姿，进而估计匹配特征点的空间三维坐标，实现对场景稀疏点云的重建，如图9所示。
基于S4确定的匹配点对，利用开源三维重建库OpenSfM估计每张图像的精确外参数、相机的内参数(包含畸变参数)，对图像执行多视图密集匹配得到深度图点云，这些数据是多光谱点云重建的基础。因此利用前述特征点匹配结果，通过运动恢复结构原理估计图像外参数、相机内参数，再通过多视图匹配方法估计得到每张图的深度图点云。
在获得每张图像的深度图点云后，结合相机内参数和图像外参数，对多光谱反射率图去畸变并从中获得深度图点云每个点的5波段光谱值，从而得到多光谱深度图点云；
然后将同一个相机所采集图像的多光谱深度图点云拼接得到相同波段的多光谱深度图点云，本实施方式中为5组不同波段光谱的点云，然后将两个相机对应的两个包含不同波段光谱的点云进行拼接，基于空间位置一致性生成完整的10波段多光谱点云，最后滤除异常点，主要工作流程如图10所示。
具体地，
由于相机镜头普遍存在一定程度的径向、切向畸变，导致图像几何失真，为了提高深度图点云光谱恢复的准确性，首先需要对反射率图像进行畸变校正。如式(14)所示，由于畸变模型的逆变换难以求解，因此将无畸变的目标图像中的坐标(x,y)通过畸变变换投影到有畸变的原图像上，以投影位置(x',y')的像素值作为无畸变图像该点像素值，
其中，r表示像点(x,y)到畸变图像像平面中心的距离，k1、k2、k3与p1、p2分别是径向与切向畸变参数。
根据图像外参数和两相机的内参数，利用式(9)所示的方法将深度图点云中的每个空间点前向投影得到每个空间点前像素坐标，进而从畸变校正后的反射率图上获得5波段光谱反射率值。
获得全部多光谱深度图点云后，按照相机的不同将这些点云分成R、B两组，分别代表本发明使用的两台多光谱相机，将每台相机所拍摄图像对应的多光谱深度图点云分别拼接融合，得到两个波段不同的5波段多光谱点云PCR和PCB。
然后基于同名点空间位置一致性从两个5波段点云中生成完整的10波段多光谱点云，如图11所示，对于点云PCR中的每个点在PCB中搜索其最近邻点，当两点距离不超过阈值d时被认为是同名点，生成一个包含两点各自5波段光谱值的新的空间点，这些新的空间点组成了10波段多光谱点云。阈值d为图像的最大地面采样距离，根据所采集数据的实际情况调整，本发明倾斜图像的最大地面采样距离约为0.08米。由于点云数据量庞大，本发明使用大规模向量相似度计算库hnswlib进行最近邻搜索，提高算法运行速度。
最后，为了去除点云噪声、滤除空间位置异常的点，使最终点云的几何结构更加平滑，本发明使用开源点云数据处理库PDAL中基于统计学的方法进行异常点识别。对多光谱点云中的每个点搜索k个近邻点并计算距离均值和标准差σ，得到式(15)所示的异常点判断阈值，
其中m为倍乘系数，k、m根据点云质量进行适当调节；
当i近邻点的距离ui≥t时将被标记为异常点并滤除。
通过本发明提出的基于双角度多光谱图像的三维多光谱点云重建方法，能够对包含复杂地物场景的倾斜多光谱图像重建得到三维点云数据，提高多维度遥感数据表示与分析能力，实现多光谱立体遥感探测应用。本发明的优势可体现在如下方面：
(a)图像特征点匹配对数多：本发明中的多光谱图像纹理增强方法能够提高图像中植被像素与非植被像素的对比度，增强图像灰度纹理、减轻多波段灰度非线性变化现象，倾斜图像投影校正方法能够将大倾角图像校正到下视视角，显著减小透视几何畸变程度，匹配像对搜索与图像局部重匹配方法能够显著减少误匹配、提高正确匹配数量。以上处理算法大幅提高了倾斜多光谱图像特征点提取与匹配对数，保证后续三维重建算法的精度和完整性。
(b)重建点云空间结构完整性高：本发明采用倾斜双角度视角采集多光谱图像，相比于单角度下视视角，能够获取更多地物侧立面影像，结合投影校正与局部重匹配方法，能够获得更多地物侧立面重建点云，最终重建场景的空间结构具有较高完整性。
(c)重建点云光谱波段完整：由于不同视角相机拍摄图像无法直接配准，现有三维重建方法不能直接从非配准相机中重建完整波段点云。本发明所设计的基于深度图的多光谱点云生成方法，首先估计每张图像的深度图点云并将相同相机的点云融合，再从全部相机点云中提取最近邻点得到每个空间点的完整光谱波段反射率值，实现完整波段多光谱点云生成。
具体实施方式二：
本实施方式为基于双角度多光谱图像的三维多光谱点云重建系统，包括图像纹理增强单元、投影校正及特征点提取单元、匹配像对搜索单元、图像匹配单元和全波段多光谱点云生成单元；其中，
图像纹理增强单元：针对不同的多光谱相机采集的图像进行纹理增强；所述不同的多光谱相机为光谱波段不同的相机，多光谱相机采集的图像为多光谱反射率图；
投影校正及特征点提取单元：基于每个多光谱图像纹理增强后的图像，针对倾斜图像，进行倾斜图像投影校正，使用SIFT算法对投影校正后的图像提取特征点，再将特征点坐标反变换回到倾斜图像中，从而确定倾斜图像中的特征点，也就是原始图像中的特征点；
匹配像对搜索单元：基于原始图像中的特征点，对原始图像进行匹配像对搜索，搜索与参考图像重叠率超过重叠率阈值的邻居图像用于后续的匹配；
图像匹配单元：针对搜索到的图像进行图像匹配；
全波段多光谱点云生成单元：首先基于确定的匹配点对估计每张图像外参数、相机的内参数，对图像执行多视图密集匹配得到深度图点云；
在获得每张图像的深度图点云后，结合相机内参数和图像外参数，对多光谱反射率图去畸变并从中获得深度图点云每个点的多波段光谱值，从而得到多光谱深度图点云；
然后将同一个相机所采集多光谱深度图点云拼接得到相同波段的多光谱深度图点云，在将不同的多光谱相机对应的包含不同波段光谱的点云进行拼接，基于空间位置一致性生成完整的全波段多光谱点云，最后滤除异常点；
所述的全波段多光谱点云为覆盖不同的多光谱相机所有波段的多光谱点云。
本实施方式所述系统为基于双角度多光谱图像的三维多光谱点云重建方法对应的程序系统，系统的各处理单元用于分别用于执行基于双角度多光谱图像的三维多光谱点云重建方法的各处理步骤所述的方法。
具体实施方式三：
本实施方式为一种计算机存储介质，所述存储介质中存储有至少一条指令，所述至少一条指令由处理器加载并执行以实现所述的一种基于双角度多光谱图像的三维多光谱点云重建方法。
应当理解，指令包括本发明描述的任何方法对应的计算机程序产品、软件或计算机化方法；所述指令可以用于编程计算机系统，或其他电子装置。计算机存储介质可以包括其上存储有指令的可读介质，可以包括但不限于磁存储介质，光存储介质；磁光存储介质包括只读存储器ROM、随机存取存储器RAM、可擦除可编程存储器(例如，EPROM和EEPROM)以及闪存层，或者适合于存储电子指令的其他类型的介质。
具体实施方式四：
本实施方式为一种基于双角度多光谱图像的三维多光谱点云重建设备，所述设备包括处理器和存储器，应当理解，包括本发明描述的任何包括处理器和存储器的设备，设备还可以包括其他通过信号或指令进行显示、交互、处理、控制等以及其他功能的单元、模块；
所述存储器中存储有至少一条指令，所述至少一条指令由处理器加载并执行以实现所述的一种基于双角度多光谱图像的三维多光谱点云重建方法。
实施例
本实施例的基于双角度多光谱图像的三维多光谱点云重建方法，如下：
设定飞行高度为80米，航向、航带重叠率均为80％，对一处地物高度差异显著、地物类型丰富的场景采集双角度多光谱图像，部分图像如图12所示，其中后视倾斜相机(以下简称“相机1”)采集200张5波段图，前视倾斜相机(以下简称“相机2”)采集185张5波段图，两相机倾斜角度均为30°，相机同时采集曝光时刻相机的经纬度坐标和姿态角信息(即位姿数据)并保存在图像的元数据区。
应用本发明中的多光谱图像纹理增强方法处理每张图像，得到对比结果如图13所示，图13为纹理增强结果及局部放大对比示意图，其中(a)为原始DN值图，(b)为(a)中1区域放大图，(c)为(a)中2区域放大图，(d)为NDVI-RGB合成图，(e)为(d)中3区域放大图，(f)为(d)中4区域放大图。由图(a)、(b)对比可知，增强后的合成图整体对比度更高，图像中点、线等几何特征明显增多，植被与非植被像素的纹理差异显著提高；由图(b)、(e)和图(c)、(f)对比可知，由于部分人造地物和裸地在该波段的光谱反射强度与周围植被反射强度相近，导致这些地物在原始DN值图中难以区分，而在NDVI-RGB合成图中，由于利用了地物的多光谱特性，增大了不同种类地物像素的差异，达到了增强图像纹理的目的，有利于提取到更多更高质量的特征点，进而提高图像匹配的成功率。
应用本发明中的倾斜图像投影校正方法，得到对比结果如图14所示，图14为图像投影校正前后对比图，其中(a)为校正前倾斜图像，(b)为校正后图像。观察可知，校正前的倾斜图像中建筑物顶面两侧边缘不平行(虚线)，而校正后变为平行(虚线)，说明投影校正减小了图像的透视几何畸变程度。
应用本发明中的匹配像对搜索方法，对一张参考图像搜索其潜在匹配图像，并将搜索结果可视化，如图15所示，其中参考图像为IMG_0986_4.GIF。
使用SIFT特征点检测器对投影校正后的图像提取特征点并应用本发明中的图像局部重匹配方法进行特征点匹配，匹配结果可视化示意图如图16所示。
与现有开源框架的特征提取与匹配结果进行对比，如表1所示。从结果可知，本发明所设计的纹理增强、投影校正、像对搜索与局部重匹配方法能够获得更多的特征点数量，并显著提高了图像匹配点对数，且能在一定程度上减小匹配点对的重投影误差，在各项指标上均取得了最优的效果，证明了本发明所设计方法的有效性。
表1图像特征点提取与匹配结果
对图像匹配点对应用本发明中的多光谱点云生成方法，获得如图17所示的多光谱密集点云模型，图17为多光谱密集点云可视化与植被指数计算结果，其中图(a)是可见光真彩色波段显示效果，图(b)是标准伪彩色显示效果，图(c)是对该多光谱点云计算的归一化差分植被指数结果，图(d)是对该多光谱点云计算的绿色叶绿素植被指数结果。从结果可以看出，本发明重建的多光谱点云空间结构较为完整，光谱波段恢复良好。
从点云中提取单个地物用于分析其空间特性的应用，例如提取图17(c)中蓝色方框内树木、建筑如图18所示。图18为单个地物点云提取结果图。其中，(a)是单棵树木点云的归一化差分植被指数，测量该树木的冠幅长度、宽度分别为8.27米、7.27米，树高为5.76米，再结合多种植被指数可进一步估计生物量、生长阶段与健康状态；(b)是一栋建筑物点云的高程着色结果，根据点云可求得该建筑物的顶面面积约为276平方米，体积约为2794立方米。上述地物的空间结构信息是难以仅通过遥感图像获取的，而这些信息对于精准农业、自然资源监测等应用具有重要意义。
总的来说，本发明所设计的基于双角度多光谱图像的三维多光谱点云重建方法能够有效处理倾斜多光谱图像，生成几何、光谱较为完整的密集点云模型，突破了传统基于图像的分析手段的不足，为许多定量遥感分析提供高质量数据，具有广阔应用前景。
本发明的上述算例仅为详细地说明本发明的计算模型和计算流程，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。
Claims (10)
1.基于双角度多光谱图像的三维多光谱点云重建方法，其特征在于，包括以下步骤：
S1、针对设置角度不同的相机采集的图像进行纹理增强；所述相机为多光谱相机，且设置角度不同的相机对应的光谱波段不同，多光谱相机采集的图像为多光谱反射率图；
S2、基于每个多光谱图像纹理增强后的图像，针对倾斜图像，进行倾斜图像投影校正，使用尺度不变特征变换匹配算法，SIFT算法，对投影校正后的图像提取特征点，再将特征点坐标反变换回到倾斜图像中，从而确定倾斜图像中的特征点，也就是原始图像中的特征点；
S3、基于原始图像中的特征点，对原始图像进行匹配像对搜索，搜索与参考图像重叠率超过重叠率阈值的图像用于后续的匹配；
S4、针对S3搜索到的图像进行图像匹配；
S5、基于S4确定的匹配点对估计每张图像外参数、相机的内参数，对图像执行多视图密集匹配得到深度图点云；
在获得每张图像的深度图点云后，结合相机内参数和图像外参数，对多光谱反射率图去畸变并从中获得深度图点云每个点的多波段光谱值，从而得到多光谱深度图点云；
然后将同一个相机所采集多光谱深度图点云拼接得到相同波段的多光谱深度图点云，再基于空间位置一致性将设置角度不同的相机对应的不同波段光谱的点云进行拼接，生成完整的全波段多光谱点云，最后滤除异常点；
所述的全波段多光谱点云为覆盖不同相机所有波段的多光谱点云。
2.根据权利要求1所述的基于双角度多光谱图像的三维多光谱点云重建方法，其特征在于，针对不同的多光谱相机采集的图像进行纹理增强包括以下步骤：
首先利用不同的多光谱相机采集图像，然后利用每个相机计算归一化差分植被指数图NDVI；
定义运算规则A＝[B⊙λb]表示对于矩阵B中的每个元素，若其值⊙λb，则与B相同维度的矩阵A中对应位置元素值为1，否则为0；λb表示一个数值，⊙表示＞、≥、<、≤符号中的一种；
利用下式去除植被指数图的阴影区域，
INDVI用于表示不同多光谱相机的植被指数图NDVI，bj用于表示不同多光谱相机的中心波长为j的单波段反射率图；
然后将植被指数图中的植被像素与蓝、绿、红三波段反射率图中的非植被像素进行合成，将植被指数值小于0.3的像素视为非植被像素，计算公式如下：
式中，分别表示多光谱图像中对应红、绿、蓝波段的反射率图像；图像I为最终获得的具有强几何纹理的NDVI-RGB合成图。
3.根据权利要求1或2所述的基于双角度多光谱图像的三维多光谱点云重建方法，其特征在于，进行倾斜图像投影校正的过程包括以下步骤：
首先确定倾斜图像校正到下视视角图像对应的透视投影变换关系，将p(x,y,1)、P′(X,Y,Z)分别表示透视投影前的图像坐标、投影后的空间坐标，有
其中，M是图像透视投影矩阵，m11至m33为图像透视投影矩阵中的元素；
设p′(x′,y′)是投影后的图像坐标，由式(3)可建立如式(4)所示的方程，
令m33＝1，则方程中未知变量变为8个，需要4组投影前后的图像坐标即可求解矩阵M；使用倾斜图像的4个顶点坐标作为p，求解对应4个投影点p′：
首先将倾斜图像中心点的像素齐次坐标o(u0,v0,1)后向投影到水平地面上，设其投影点的世界坐标为根据位姿数据可求得倾斜图像在世界坐标系中的旋转矩阵R1；
然后将每个相机曝光位置的经纬高坐标转换为地心固定直角坐标，获得图像在世界坐标系中的平移向量t1；
基于图像在世界坐标系中的旋转矩阵R1，根据式(7)求得图像中心点的地面投影点：
式中，K1为倾斜图像的相机内参数；D代表相机光心到地面的高度，即投影点的深度；设初始深度D为无人机飞行高度，当时更新深度/>并利用式(7)重新计算，α1为世界坐标Z轴基准阈值，直到/>时即求得地面投影点坐标/>同时获得其深度D；
将倾斜图像的顶点坐标投影到地面得到的地面点坐标
然后设校正图像的滚转角、俯仰角均为0，航偏角与倾斜图像保持相同，即(α,β,γ)＝(0,0,γcorr)，并求得校正图像在世界坐标系下的旋转矩阵为R2；
记校正图像在世界坐标为oC(XW,YW,h)，令h＝D，得到校正图像的平移向量t2；
最后，利用式(9)将4个地面点前向投影到校正图像上即可得到对应的4个图像投影点；
p′＝K1[R2pg+t2] (9)
将该4组投影前后的图像坐标代入式(4)中可得到8个方程，求解方程获得图像投影变换矩阵M并应用式(3)对倾斜图像进行变换，将变换后的空间坐标P′(X,Y,Z)归一化为图像坐标，即至此完成了倾斜图像的投影校正。
4.根据权利要求3所述的基于双角度多光谱图像的三维多光谱点云重建方法，其特征在于，基于原始图像中的特征点，对原始图像进行匹配像对搜索的过程包括以下步骤：
首先基于每张原始图像的位姿数据求得各自的地心固定直角坐标同时确定原始图像该数据集中最小的经度lonref、纬度latref，令参考高度为altref＝0的空间点为参考点，并求得参考点的地心固定直角坐标为/>由式(10)得到每张原始图像在世界坐标系下的坐标，同时分别求解世界坐标系中的旋转矩阵R和平移向量t；
然后将每张图像的中心点后向投影到地面获得投影点OG(XG,YG,0)；
设某一个相机在OC位置拍摄的一张后视倾斜图像的中心点的地面投影点坐标为OG，则参考图像的邻居匹配图像投影点O′G满足式(11)：
||OG-O′G||≤rs (11)
其中，rs为搜索半径，根据图像重叠率估计得到；
搜索与参考图像重叠率超过重叠率阈值的图像用于后续的匹配。
5.根据权利要求4所述的基于双角度多光谱图像的三维多光谱点云重建方法，其特征在于，针对S3搜索到的图像进行图像匹配的过程包括以下步骤：
首先，使用KD-Tree算法对参考图像中的特征点在匹配图像上搜索最近邻和次近邻，设最近邻、次近邻与参考图像中的特征点的欧式距离分别为d1、d2，当二者满足式(13)时参考图像中的特征点与最近邻匹配，即，同时将r作为该匹配点对的置信度；
其中，η是匹配点筛选阈值；
针对匹配结果，使用MAGSAC++算法剔除误匹配点对，在剔除误匹配点对的过程中，采用最大化一致性原则估计单应变换矩阵；
在估计得到进行匹配的两张图像之间的单应变换矩阵Md之后，将匹配图像的特征点坐标投影到参考图像坐标系中，并对参考图像分块，每一个块为一个子图；然后遍历参考图像的每一个子图，进行局部重匹配，即对子图区域内的参考特征点与匹配特征点执行KD-Tree近邻搜索匹配，然后合并所有子图的匹配结果；
交换参考图像与匹配图像，同时利用Md -1对新的匹配图像特征点投影，对新的参考图像分块并执行局部重匹配；
取双向匹配结果的交集作为两图像的匹配结果；
最后针对匹配结果采用MAGSAC++算法去除误匹配点对，得到最终的匹配点对。
6.根据权利要求5所述的基于双角度多光谱图像的三维多光谱点云重建方法，其特征在于，S5中基于S4确定的匹配点对估计每张图像外参数、相机的内参数的过程是通过运动恢复结构原理实现的。
7.根据权利要求6所述的基于双角度多光谱图像的三维多光谱点云重建方法，其特征在于，S5所述最后滤除异常点的过程中，使用开源点云数据处理库PDAL中基于统计学的方法进行异常点识别：
对多光谱点云中的每个点搜索k个近邻点并计算距离均值和标准差σ，得到异常点判断阈值：
其中，m为倍乘系数；
当i近邻点的距离ui≥t时将被标记为异常点并滤除。
8.基于双角度多光谱图像的三维多光谱点云重建系统，其特征在于，包括图像纹理增强单元、投影校正及特征点提取单元、匹配像对搜索单元、图像匹配单元和全波段多光谱点云生成单元；其中，
图像纹理增强单元：针对不同的多光谱相机采集的图像进行纹理增强；所述不同的多光谱相机为光谱波段不同的相机，多光谱相机采集的图像为多光谱反射率图；
投影校正及特征点提取单元：基于每个多光谱图像纹理增强后的图像，针对倾斜图像，进行倾斜图像投影校正，使用SIFT算法对投影校正后的图像提取特征点，再将特征点坐标反变换回到倾斜图像中，从而确定倾斜图像中的特征点，也就是原始图像中的特征点；
匹配像对搜索单元：基于原始图像中的特征点，对原始图像进行匹配像对搜索，搜索与参考图像重叠率超过重叠率阈值的邻居图像用于后续的匹配；
图像匹配单元：针对搜索到的图像进行图像匹配；
全波段多光谱点云生成单元：首先基于确定的匹配点对估计每张图像外参数、相机的内参数，对图像执行多视图密集匹配得到深度图点云；
在获得每张图像的深度图点云后，结合相机内参数和图像外参数，对多光谱反射率图去畸变并从中获得深度图点云每个点的多波段光谱值，从而得到多光谱深度图点云；
然后将同一个相机所采集多光谱深度图点云拼接得到相同波段的多光谱深度图点云，在将不同的多光谱相机对应的包含不同波段光谱的点云进行拼接，基于空间位置一致性生成完整的全波段多光谱点云，最后滤除异常点；
所述的全波段多光谱点云为覆盖不同的多光谱相机所有波段的多光谱点云。
9.一种计算机存储介质，其特征在于，所述存储介质中存储有至少一条指令，所述至少一条指令由处理器加载并执行以实现如权利要求1至7任意一项所述的一种基于双角度多光谱图像的三维多光谱点云重建方法。
10.一种基于双角度多光谱图像的三维多光谱点云重建设备，其特征在于，所述设备包括处理器和存储器，所述存储器中存储有至少一条指令，所述至少一条指令由处理器加载并执行以实现如权利要求1至7任意一项所述的一种基于双角度多光谱图像的三维多光谱点云重建方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310823147.XA CN116758223A (zh) | 2023-07-06 | 2023-07-06 | 基于双角度多光谱图像的三维多光谱点云重建方法、系统及设备 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310823147.XA CN116758223A (zh) | 2023-07-06 | 2023-07-06 | 基于双角度多光谱图像的三维多光谱点云重建方法、系统及设备 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116758223A true CN116758223A (zh) | 2023-09-15 |
Family
ID=87949569
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310823147.XA Pending CN116758223A (zh) | 2023-07-06 | 2023-07-06 | 基于双角度多光谱图像的三维多光谱点云重建方法、系统及设备 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116758223A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117557617A (zh) * | 2024-01-12 | 2024-02-13 | 山东师范大学 | 一种基于平面先验优化的多视密集匹配方法、系统及设备 |
-
2023
- 2023-07-06 CN CN202310823147.XA patent/CN116758223A/zh active Pending
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117557617A (zh) * | 2024-01-12 | 2024-02-13 | 山东师范大学 | 一种基于平面先验优化的多视密集匹配方法、系统及设备 |
| CN117557617B (zh) * | 2024-01-12 | 2024-04-09 | 山东师范大学 | 一种基于平面先验优化的多视密集匹配方法、系统及设备 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109613513B (zh) | 一种顾及InSAR形变因子的光学遥感潜在滑坡自动识别方法 | |
| Liu et al. | Deep convolutional neural network training enrichment using multi-view object-based analysis of Unmanned Aerial systems imagery for wetlands classification | |
| CN109146948B (zh) | 基于视觉的作物长势表型参数量化与产量相关性分析方法 | |
| CN107451982B (zh) | 一种基于无人机影像的高郁闭度林分树冠面积获取方法 | |
| Jurado et al. | Remote sensing image fusion on 3D scenarios: A review of applications for agriculture and forestry | |
| CN103337052B (zh) | 面向宽幅遥感影像的自动几何纠正方法 | |
| Marí et al. | Sat-nerf: Learning multi-view satellite photogrammetry with transient objects and shadow modeling using rpc cameras | |
| CN108759788B (zh) | 无人机影像定位定姿方法及无人机 | |
| López et al. | An optimized approach for generating dense thermal point clouds from UAV-imagery | |
| Karsli et al. | Automatic building extraction from very high-resolution image and LiDAR data with SVM algorithm | |
| CN116758223A (zh) | 基于双角度多光谱图像的三维多光谱点云重建方法、系统及设备 | |
| Ribera et al. | Estimating phenotypic traits from UAV based RGB imagery | |
| Liu et al. | Robust radiometric normalization of multitemporal satellite images via block adjustment without master images | |
| CN111275616B (zh) | 低空航拍图像拼接方法和装置 | |
| CN113066173B (zh) | 三维模型构建方法、装置和电子设备 | |
| Morelli et al. | Deep-image-matching: a toolbox for multiview image matching of complex scenarios | |
| CN112991186B (zh) | 一种无人机大视场高光谱图像生成方法及系统 | |
| CN113516059B (zh) | 固体废弃物的识别方法、装置、电子设备及存储介质 | |
| Atik et al. | An automatic image matching algorithm based on thin plate splines | |
| Sjahputera et al. | GeoCDX: An automated change detection & exploitation system for high resolution satelite imagery | |
| Alobeid | Assessment of matching algorithms for urban DSM generation from very high resolution satellite stereo images | |
| Zhu et al. | Generation of thermal point clouds from uncalibrated thermal infrared image sequences and mobile laser scans | |
| Gonçalves | Using structure-from-motion workflows for 3D mapping and remote sensing | |
| CN115424022B (zh) | 输电走廊地面点云分割方法、装置和计算机设备 | |
| Yu et al. | An integrative object-based image analysis workflow for UAV images |
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 |