CN117495740A - 一种基于无人机航拍遥感图像自动处理方法 - Google Patents

一种基于无人机航拍遥感图像自动处理方法 Download PDF

Info

Publication number
CN117495740A
CN117495740A CN202311449087.6A CN202311449087A CN117495740A CN 117495740 A CN117495740 A CN 117495740A CN 202311449087 A CN202311449087 A CN 202311449087A CN 117495740 A CN117495740 A CN 117495740A
Authority
CN
China
Prior art keywords
image
slice
pixel
longitude
latitude
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202311449087.6A
Other languages
English (en)
Inventor
林少宇
朱定
梁德成
张远辉
黄兵祥
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Guangxi Zhuang Autonomous Region Public Information Industry Co ltd
Original Assignee
Guangxi Zhuang Autonomous Region Public Information Industry Co ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Guangxi Zhuang Autonomous Region Public Information Industry Co ltd filed Critical Guangxi Zhuang Autonomous Region Public Information Industry Co ltd
Priority to CN202311449087.6A priority Critical patent/CN117495740A/zh
Publication of CN117495740A publication Critical patent/CN117495740A/zh
Pending legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/60—Analysis of geometric attributes
    • G06T7/62—Analysis of geometric attributes of area, perimeter, diameter or volume
    • 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
    • G06T7/00—Image analysis
    • G06T7/70—Determining position or orientation of objects or cameras
    • G06T7/73—Determining position or orientation of objects or cameras using feature-based methods
    • G06T7/75—Determining position or orientation of objects or cameras using feature-based methods involving models
    • 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
    • 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/30—Subject of image; Context of image processing
    • G06T2207/30244—Camera pose

Landscapes

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

Abstract

本发明公开了一种基于无人机航拍遥感图像自动处理方法，通过网联无人机在飞行过程中就实时将拍摄图像自动上传至对象存储服务器；自动开始下载完整的任务遥感图像；启动图片处理程序合成正射影像；按坐标参数要求自动调用切图程序，并将正射影像切割为地图可显示文件，根据任务名自动将文件上传至服务器，自动解压切图文件；处理得到的地图切片全景图像，可供无人机巡检平台根据任务查看该任务瓦片，并加载到地图上。通过本方法技术极大地节省任务时间，通过合成后的整体拼接合同大图进行隐患排查，直观减少原来对飞行任务全量相片进行隐患排查的工作模式，任务时间由之前的3‑4天可缩短至1天，有效的提高工作效率和节省了人工成本。

Description

一种基于无人机航拍遥感图像自动处理方法
技术领域
本发明属于无人机航拍遥感图像自动处理技术领域，具体涉及一种基于无人机航拍遥感图像自动处理方法。
背景技术
随着科技的发展，无人机的使用范围逐渐扩大，当前无人机主要应用在城市管理、巡检等领域，行业用户已经不仅仅关注于无人机本身，而是对无人机的功能提出了整体化需求，亟待解决。现阶段存在的问题主要是无人机航拍任务结束后，需人工将拍摄的图片导入上传，再将图片进行处理，人工成本高。并且人工干预图片处理容易受到知识技术等原因的限制，导致错误率上升。为解决以上问题，本发明旨在解决无人机航拍图形自动处理技术方面的缺失，助力无人机行业快速发展；解决了当前无人机市场在航拍任务结束后需人工手动进行图传、人工处理图片、准确率低，任务执行过程不透明等问题。
发明内容
本发明针对现有技术的不足，提供一种无需人工干预并且有效提高工作效率和节约人力成本的基于无人机航拍遥感图像自动处理方法。
本发明的技术方案如下：
一种基于无人机航拍遥感图像自动处理方法，包括以下步骤：
步骤一、通过网联无人机在飞行过程中就实时将拍摄图像自动上传至对象存储服务器；
步骤二、飞行任务完成后自动开始下载完整的任务遥感图像；
步骤三、启动图片处理程序，自动合成正射影像；
步骤四、按国家自然资源局地图2000的坐标参数要求自动调用切图程序，并将正射影像切割为地图可显示文件，根据任务名自动将文件上传至服务器，自动解压切图文件；
步骤五、以上步骤处理得到的地图切片全景图像，可供无人机巡检平台根据任务查看该任务瓦片，并加载到地图上。
对本发明的进一步说明，所述方法还包括使用无人机拍摄遥感图像结合自动化流程对遥感图像进行自动化处理；读取遥感图像同归对比图像中的参考物，将具有相同图像的像素移除，拼接不同的图像，合成正射影像；再将影像投影至web墨卡托，获取遥感图像经纬范围，计算影像像素分辨率，坐标转wgs84经纬度；将原始影影像地理范围与指定缩放级别指定行列号的切片交集，计算交集的像素信息，将切片数据写入文件。
对本发明的进一步说明，所述图片处理程序自动合成正射影像具体包括：
（1）根据任务号从对象存储服务器下载相应任务遥感图像；
（2）对图像进行去畸变处理；
（3）特征提取；
（4）跟踪当前帧，并添加关键帧；
（5）环路检测：提取不同特征点的集合中的特征点和特征描述子；通过特征匹配算法对特征点进行匹配，并根据匹配结果进行几何一致性验证；最后，根据验证结果判断是否存在环路连接，如果存在，则进行局部优化处理；
所述局部优化通过优化相机位姿和特征点的三维坐标，调整关键帧之间的对齐关系，使用非线性优化算法Levenberg-Marquardt来最小化重投影误差，使得特征点在各个关键帧中的投影位置更加一致；
所述经过环路检测和局部优化后的关键帧集合，通过特征点的三维坐标和相机位姿，进行平面拟合，找到能够最好拟合关键帧上特征点的平面，使用最小二乘法算法，得到平面的参数方程；最后，将图像序列中的特征点投影到拟合的平面上，获得更好的对齐效果和视觉连续性；
（6）图像变换：根据相机位姿和平面参数，对图像进行透视变换，将图像投影到平面上，得到正交视图；
（7）图像拼接、融合：根据相机位姿和特征点的三维坐标，进行图像拼接，将多个图像拼接成全景图；所述进行图像融合：通过加权平均、多重曝光，将拼接后的图像进行平滑过渡和细节保留；
通过以上处理得到二维正射影像用于后续处理。
对本发明的进一步说明，所述对图像进行去畸变处理具体为：
（1）创建目标点集合dst；
（2）从相机矩阵 cameraMatrix 和畸变系数 distortionCoeff 中获取相机参数：焦距（fx、fy），光心坐标（ux、uy），径向畸变系数（k1、k2、k3）和切向畸变系数（p1、p2）；
（3）遍历输入的源点坐标数组 src ；
（4）对于每个源点坐标，首先进行坐标转换，将其转换为相对于光心的归一化坐标：xDistortion = (p.x - ux) / fx，yDistortion = (p.y - uy) / fy；
（5）使用迭代的方式进行去畸变的求解。通过设定初始值（x0、y0），根据畸变模型的公式进行迭代计算 ：
畸变公式为：r2 = xDistortion*xDistortion + yDistortion*yDistortion
distRadialA = 1 / (1. + k1 * r2 + k2 * r2 * r2 + k3 * r2 * r2 * r2)
distRadialB = 1. + k4 * r2 + k5 * r2 * r2 + k6 * r2 * r2 * r2
deltaX = 2. * p1 * xDistortion * yDistortion + p2 * (r2 + 2. *xDistortion * xDistortion)
deltaY = p1 * (r2 + 2. * yDistortion * yDistortion) + 2. * p2 *xDistortion * yDistortion
xCorrected = (x0 - deltaX)* distRadialA * distRadialB
yCorrected = (y0 - deltaY)* distRadialA * distRadialB
xDistortion = xCorrected
yDistortion = yCorrected
其中，xDistortion和yDistortion表示归一化坐标系下的畸变坐标，xCorrected和yCorrected表示去畸变后的坐标，ux和uy表示相机光心坐标，fx和fy表示相机的焦距，k1、k2、k3、k4、k5、k6、p1和p2是相机的畸变系数，这些参数通过cameraMatrix和distortionCoeff传递给函数，最终得到xDistortion和yDistortion畸变坐标，用于下一步坐标变换；
（6）进行坐标变换，将去畸变后的坐标转换为相机坐标系下的像素坐标：xCorrected = xCorrected * fx + ux，yCorrected = yCorrected * fy + uy；
（7）将去畸变后的坐标（xCorrected、yCorrected）添加到目标点坐标数组 dst中。
（8）循环结束后， dst 中将存储了所有源点坐标经过去畸变处理后的结果。
通过以上处理得到了去畸变后的像素坐标，用于特征提取。
对本发明的进一步说明，所述特征提取具体为：
（1）获取像素坐标矩阵A的行数和列数，分别赋值给变量r和c；
（2）计算样本矩阵A的均值，将结果存储在变量meanVec中；
（3）计算协方差矩阵的转置covMatT：通过将样本矩阵A减去均值向量meanVec并进行矩阵乘法得到；
（4）使用eigs函数计算covMatT的前k个特征值和对应的特征向量，将结果分别存储在变量V和D中；所述eigs函数的主要参数包括输入矩阵、所需计算的特征值数量（k值）以及可选参数；通过调用eigs函数，可以得到指定数量（k个）的特征值和对应的特征向量；
（5）将协方差矩阵(covMatT)'的特征向量V乘以样本矩阵A减去均值向量meanVec的转置Z，得到变换后的特征向量；
（6）对变换后的特征向量进行归一化处理，使其成为单位特征向量；
（7）使用线性变换（投影）将样本矩阵A减去均值向量meanVec的转置Z乘以变换后的特征向量V，得到降维至k维的样本特征向量组成的矩阵pcaA。
对本发明的进一步说明，所述跟踪当前帧，并添加关键帧具体为：
（1）提取矩阵区域的特征点作为初始特征点集合；
（2）特征点跟踪：将初始特征点作为输入，使用cv2.calcOpticalFlowPyrLK()函数在下一帧中跟踪这些特征点。该函数会返回跟踪到的特征点的新位置；
（3）更新目标区域：根据跟踪到的新特征点的位置，更新目标区域的位置和大小；
（4）关键帧判断：根据特征点数量的变化、目标区域的变化来判断当前帧是否满足添加关键帧的条件；
（5）添加关键帧：如果当前帧被判断为关键帧，将当前帧作为关键帧，并将当前帧的特征点集合添加到关键帧集合中；
（6）局部优化：对关键帧集合进行局部优化，通过优化特征点的位置和相机位姿，提高整体的几何一致性；
重复执行步骤（2）到步骤（6），以连续跟踪目标并添加关键帧，将得到具有不同特征点的集合。
对本发明的进一步说明，所述图像拼接完成后使用PIL去除正射影像黑边，此时图中坐标丢失；
所述合成正射影像的tif图片存在黑色边界背景，切片完成后在地图显示也会出现黑色边界，需要进行处理，具体为：
（1）使用python的PIL库设置允许处理大尺寸图像并加载被截断的图像；
（2）打开了指定路径的图像文件，并将其转换为RGBA模式；
（3）获取图像的宽度L和高度H，获取图像左上角像素的颜色color_0，然后对图像的每个像素进行遍历，获取像素的颜色color_1；如果color_1与color_0相同，则将color_1的透明度设置为0，并将修改后的颜色重新赋值给该像素；
（4）保存处理后的图像到指定路径，并返回输出路径。
对本发明的进一步说明，所述还包括使用gdal库还原坐标，将带有完整坐标正射影像的仿射矩阵、地图投影信息写入去除黑边后的正射影像；所述除去黑边后导致丢失地理坐标信息的正射影像，通过以下处理还原完整坐标正射影像的仿射矩阵、地图投影信息，具体为：
（1）根据去除黑边后返回的输出路径打开栅格数据集文件；
（2）读取栅格数据集的数组表示（img_array）；
（3）根据数组的数据类型，确定输出文件的数据类型（datatype）；
（4）确定图像的波段数（img_bands）以及图像的高度和宽度（im_height和im_width）；
（5）创建输出文件驱动（driver）和数据集（dataset），并设置输出文件的大小、波段数、数据类型信息；
（6）设置输出文件的仿射变换参数（img_transf）和投影信息（img_proj）（此信息通过读取原始正射影像文件中带有的完整坐标正射影像的仿射矩阵、地图投影信息获得）；
（7）将图像数据写入输出文件的相应波段；
（8）返回输出文件路径（outPath）；
所述除去黑边后导致丢失地理坐标信息的正射影像，通过以上处理还原完整坐标正射影像的仿射矩阵、地图投影信息。
对本发明的进一步说明，所述正射影像处理完成后使用gdal进行切片，包括以下步骤：
（1）将原始影像数据集重投影到Web墨卡托投影坐标系，并将结果保存到临时文件dataset中即打开指定路径的TIFF文件，并将其重投影到Web墨卡托投影坐标系，具体为：
S1创建一个临时文件路径tmpFile 用于存储GTiff格式的输出数据集；
S2创建一个SpatialReference对象，用于存储空间参考信息；
S3使用spatialReference的ImportFromEPSG方法导入Web墨卡托投影坐标系的EPSG代码（WEB_MERCATOR_EPSG_CODE）；
S4使用gdal.AutoCreateWarpedVRT方法创建一个虚拟数据集webMercatorDs，将原始影像数据集（srcDs）重投影到Web墨卡托投影坐标系；
S5使用gdal.GetDriverByName("GTiff").CreateCopy方法创建一个与webMercatorDs相同的GTiff格式的输出数据集，并将其保存到tmpFile路径中；
S6删除原始影像数据集（srcDs）和虚拟数据集（webMercatorDs）；
（2）获取重投影后影像的地理范围和像素分辨率，并计算经纬度范围，获取原始影像在WGS84经纬度坐标系下的地理范围；具体为：
S1使用dataset.GetGeoTransform()方法获取原始影像的地理变换参数（geoTransform），包括地理坐标系的起始经度、像素的宽度、像素的高度、地理坐标系的起始纬度、像素的宽度和像素的高度；
S2获取原始影像的像素分辨率，通过dataset.getRasterXSize()和dataset.getRasterYSize()方法获取原始影像的宽度和高度；
S3 根据地理变换参数计算原始影像的经纬度范围，经度范围由起始经度（lngMin）和宽度乘以像素宽度（xSize * geoTransform[1]）以及高度乘以像素高度（ySize* geoTransform[2]）得到；纬度范围由起始纬度（latMax）和宽度乘以像素宽度（xSize *geoTransform[4]）以及高度乘以像素高度（ySize * geoTransform[5]）得到；
S4 得到的结果是原始影像的经纬度范围（lngMin, latMax, lngMax, latMin）；
（3）将原始影像的地理范围从Web墨卡托投影坐标系转换为WGS84经纬度坐标系，以便切片完成的瓦片可在天地图上展示，具体为：
S1使用EPSG代码解码获取源坐标参考系统（sourceCRS）和目标坐标参考系统（targetCRS）；其中，源坐标参考系统为Web墨卡托投影坐标系，目标坐标参考系统为WGS84经纬度坐标系；
S2 使用源坐标参考系统和目标坐标参考系统，通过CRS.findMathTransform方法获取坐标转换对象（transform）；
S3 创建Coordinate对象，用于存储经纬度的最大值和最小值；
S4 使用JTS.transform方法，将原始影像的地理范围的最大经纬度（lngMax,latMax）和最小经纬度（lngMin, latMin）从源坐标参考系统转换为目标坐标参考系统，并存储到新的Coordinate对象（lngLatMax, lngLatMin）中；
（4）根据指定的缩放级别，计算切片的行列号，具体为：
S1计算切片的x坐标：
int x = (int) Math.floor((lng + 180) / 360 * (1 << zoom));
公式中，lng为经度，zoom为缩放级别；首先将经度lng转换为切片内的相对经度，然后乘以切片数量(1 << zoom)，再向下取整得到切片的x坐标；
S2计算切片的y坐标：
int y = (int) Math.floor((1 - Math.log(Math.tan(Math.toRadians(lat))+ 1 / Math.cos(Math.toRadians(lat))) / Math.PI) / 2 * (1 << zoom));
公式中，lat为纬度，zoom为缩放级别；首先将纬度lat转换为弧度，然后使用数学公式计算切片内的相对纬度，再乘以切片数量(1 << zoom)，最后向下取整得到切片的y坐标；
S3对坐标进行边界处理：如果x坐标小于0，则将其设置为0；如果x坐标大于等于切片数量(1 << zoom)，则将其设置为切片数量减1；同样的，如果y坐标小于0或大于等于切片数量(1 << zoom)，也进行相应的边界处理；
最后得到行列号用于计算经纬度范围；
（5）根据切片的行列号和缩放级别，计算切片的经纬度范围，具体为：
S1接收两个参数：切片的列号x和缩放级别z；
S2使用公式 x / Math.pow(2.0, z) * 360.0 - 180 将切片的列号x转换为对应的经度值；
- Math.pow(2.0, z) 表示2的z次方，即切片在当前缩放级别下的总列数；
- x / Math.pow(2.0, z) 将切片的列号x转换为相对于总列数的比例；
- x / Math.pow(2.0, z) * 360.0 将相对比例转换为经度范围，取值范围为0到360；
- x / Math.pow(2.0, z) * 360.0 - 180 将经度范围从0到360转换为-180到180，得到最终的经度值；
S3返回转换后的经度值；
（6）根据切片的地理范围、像素分辨率以及原始影像的地理范围，计算出原始影像在切片内的像素大小和偏移像素坐标，用于后续的切片生成和裁剪操作，具体为：
根据原始影像的地理范围与切片的交集，确定切片的像素分辨率和范围。
S1创建了切片的地理范围（tileBound）和原始影像的地理范围（imageBound），都使用了WGS84经纬度坐标系作为参考坐标系；
S2使用tileBound和imageBound计算它们的交集（intersect），得到切片和原始影像共同的地理范围；
S3计算切片的像素分辨率：
- 切片的东西方向像素分辨率（dstWePixelResolution）：切片的经度范围除以切片的像素宽度（256）；
- 切片的南北方向像素分辨率（dstNsPixelResolution）：切片的纬度范围除以切片的像素高度（256）；
S4计算交集的像素信息：
- offsetX和offsetY：计算在切片地理范围内，原始影像的起始点像素坐标。通过将交集起始点的经度和纬度减去原始影像的起始经度和纬度，并除以原始影像的像素分辨率得到；
- blockXSize和blockYSize：计算在切片地理范围内，原始影像的像素大小。通过将交集的经度和纬度范围除以原始影像的像素分辨率得到；
- imageXBuf和imageYBuf：计算原始影像在切片内的像素大小；通过将交集的经度和纬度范围除以切片的像素分辨率得到，并向上取整；
- imageOffsetX和imageOffsetY：计算原始影像在切片中的偏移像素坐标。通过将交集起始点的经度和纬度减去切片的起始经度和纬度，并除以切片的像素分辨率得到。如果计算结果小于0，则将其设为0；
S5得到的结果是切片的像素分辨率和原始影像在切片内的像素大小和偏移像素坐标；
（7）使用GDAL的ReadRaster方法读取原始影像在切片范围内的像素数据到对应的缓冲区数组，具体为：
S1获取原始影像数据集的第1、2、3个波段（Band）对象，分别赋值给inBand1、inBand2和inBand3变量；
S2创建用于存储波段数据的缓冲区数组：
- band1BuffData数组用于存储第1个波段的数据，数组大小为TILE_SIZE *TILE_SIZE * gdalconst.GDT_Int32；
- band2BuffData数组用于存储第2个波段的数据，数组大小同上；
- band3BuffData数组用于存储第3个波段的数据，数组大小同上；
S3使用inBand1的ReadRaster方法从原始影像的第1个波段中读取像素数据，并将数据存储到band1BuffData数组中。读取的范围是从原始影像的偏移像素坐标(offsetX,offsetY)开始，读取blockXSize * blockYSize大小的像素块，并将其写入到imageXBuf *imageYBuf大小的缓冲区中；
S4 同样的方式，使用inBand2的ReadRaster方法从原始影像的第2个波段中读取像素数据，并将数据存储到band2BuffData数组中；
S5同样的方式，使用inBand3的ReadRaster方法从原始影像的第3个波段中读取像素数据，并将数据存储到band3BuffData数组中；
（8）创建内存驱动器，并将图像数据写入内存中；通过创建了一个内存数据集（msmDS），并获取了该数据集的各个波段对象，以便后续的像素数据写入和处理操作，具体为：
S1 获取名为"MEM"的GDAL驱动器（Driver）对象，并将其赋值给memDriver变量；
S2使用memDriver的Create方法创建一个名为"msmDS"的内存数据集（Dataset），该数据集的大小为256x256像素，有4个波段；
S3获取msmDS数据集的第1、2、3个波段（Band）对象，分别赋值给dstBand1、dstBand2和dstBand3变量；
（9）为图像数据设置alpha通道，实现背景透明效果；通过将原始影像的各个波段数据和alpha波段数据写入到内存数据集msmDS中，并将生成的切片文件保存到指定目录中，以便后续在Web地图上展示，具体为：
S1获取msmDS数据集的第4个波段（alphaBand）对象，并将其赋值给alphaBand变量；
S2创建一个用于存储alpha波段数据的缓冲区数组alphaData，数组大小为256x256xgdalconst.GDT_Int32；
S3使用循环遍历alphaData数组的每个元素，如果对应的band1BuffData数组的元素大于0，则将alphaData数组的元素设置为255，否则保持为0；
S4使用dstBand1、dstBand2、dstBand3和alphaBand的WriteRaster方法将相应的波段数据写入到msmDS数据集中；写入的范围是从imageOffsetX、imageOffsetY开始，写入imageXBuf x imageYBuf大小的像素块；
S5创建一个目录对象dir，用于存储切片文件；
S6创建切片文件的路径pngPath，路径格式为outputDir/Z[zoom]/[row]/[col].png。
本发明的有益效果：
通过本方法技术极大地节省任务时间，通过合成后的整体拼接合同大图进行隐患排查，直观减少原来对飞行任务全量相片进行隐患排查的工作模式，任务时间由之前的3-4天可缩短至1天，有效的提高工作效率和节省了人工成本。
本技术所具备的优点：
1、平台自动化能力提升：只需进行飞行任务规划，无人机执行完飞行任务，自动对图片进行处理，为企业节省人工成本；
2、进度可视化：自动化任务执行时可通过后台管理系统查看实时进度；
3、准确率提升：无需人工干预，避免了因人员知识技术等原因导致的错误；
4、实时图传：实时4/5G网联传图，提升工作效率的同时，提升通讯领域在无人机应用的占比。
附图说明
图1为本发明方法的流程示意图。
实施方式
下面结合具体实施例对本发明做进一步说明。
实施例
一种基于无人机航拍遥感图像自动处理方法，包括以下步骤：
步骤一、通过网联无人机在飞行过程中就实时将拍摄图像自动上传至对象存储服务器；
步骤二、飞行任务完成后自动开始下载完整的任务遥感图像；
步骤三、启动图片处理程序，自动合成正射影像；
步骤四、按国家自然资源局地图2000的坐标参数要求自动调用切图程序，并将正射影像切割为地图可显示文件，根据任务名自动将文件上传至服务器，自动解压切图文件；
步骤五、以上步骤处理得到的地图切片全景图像，可供无人机巡检平台根据任务查看该任务瓦片，并加载到地图上。
对本发明的进一步说明，所述方法还包括使用无人机拍摄遥感图像结合自动化流程对遥感图像进行自动化处理；读取遥感图像同归对比图像中的参考物，将具有相同图像的像素移除，拼接不同的图像，合成正射影像；再将影像投影至web墨卡托，获取遥感图像经纬范围，计算影像像素分辨率，坐标转wgs84经纬度；将原始影影像地理范围与指定缩放级别指定行列号的切片交集，计算交集的像素信息，将切片数据写入文件。
所述图片处理程序自动合成正射影像具体包括：
（1）根据任务号从对象存储服务器下载相应任务遥感图像；
（2）对图像进行去畸变处理；
（3）特征提取；
（4）跟踪当前帧，并添加关键帧；
（5）环路检测：提取不同特征点的集合中的特征点和特征描述子；通过特征匹配算法对特征点进行匹配，并根据匹配结果进行几何一致性验证；最后，根据验证结果判断是否存在环路连接，如果存在，则进行局部优化处理；
所述局部优化通过优化相机位姿和特征点的三维坐标，调整关键帧之间的对齐关系，使用非线性优化算法Levenberg-Marquardt来最小化重投影误差，使得特征点在各个关键帧中的投影位置更加一致；
所述经过环路检测和局部优化后的关键帧集合，通过特征点的三维坐标和相机位姿，进行平面拟合，找到能够最好拟合关键帧上特征点的平面，使用最小二乘法算法，得到平面的参数方程；最后，将图像序列中的特征点投影到拟合的平面上，获得更好的对齐效果和视觉连续性；
（6）图像变换：根据相机位姿和平面参数，对图像进行透视变换，将图像投影到平面上，得到正交视图；
（7）图像拼接、融合：根据相机位姿和特征点的三维坐标，进行图像拼接，将多个图像拼接成全景图；所述进行图像融合：通过加权平均、多重曝光，将拼接后的图像进行平滑过渡和细节保留；
通过以上处理得到二维正射影像用于后续处理。
所述对图像进行去畸变处理具体为：
（1）创建目标点集合dst；
（2）从相机矩阵 cameraMatrix 和畸变系数 distortionCoeff 中获取相机参数：焦距（fx、fy），光心坐标（ux、uy），径向畸变系数（k1、k2、k3）和切向畸变系数（p1、p2）；
（3）遍历输入的源点坐标数组 src ；
（4）对于每个源点坐标，首先进行坐标转换，将其转换为相对于光心的归一化坐标：xDistortion = (p.x - ux) / fx，yDistortion = (p.y - uy) / fy；
（5）使用迭代的方式进行去畸变的求解。通过设定初始值（x0、y0），根据畸变模型的公式进行迭代计算 ：
畸变公式为：r2 = xDistortion*xDistortion + yDistortion*yDistortion
distRadialA = 1 / (1. + k1 * r2 + k2 * r2 * r2 + k3 * r2 * r2 * r2)
distRadialB = 1. + k4 * r2 + k5 * r2 * r2 + k6 * r2 * r2 * r2
deltaX = 2. * p1 * xDistortion * yDistortion + p2 * (r2 + 2. *xDistortion * xDistortion)
deltaY = p1 * (r2 + 2. * yDistortion * yDistortion) + 2. * p2 *xDistortion * yDistortion
xCorrected = (x0 - deltaX)* distRadialA * distRadialB
yCorrected = (y0 - deltaY)* distRadialA * distRadialB
xDistortion = xCorrected
yDistortion = yCorrected
其中，xDistortion和yDistortion表示归一化坐标系下的畸变坐标，xCorrected和yCorrected表示去畸变后的坐标，ux和uy表示相机光心坐标，fx和fy表示相机的焦距，k1、k2、k3、k4、k5、k6、p1和p2是相机的畸变系数，这些参数通过cameraMatrix和distortionCoeff传递给函数，最终得到xDistortion和yDistortion畸变坐标，用于下一步坐标变换；
（6）进行坐标变换，将去畸变后的坐标转换为相机坐标系下的像素坐标：xCorrected = xCorrected * fx + ux，yCorrected = yCorrected * fy + uy；
（7）将去畸变后的坐标（xCorrected、yCorrected）添加到目标点坐标数组 dst中。
（8）循环结束后， dst 中将存储了所有源点坐标经过去畸变处理后的结果。
通过以上处理得到了去畸变后的像素坐标，用于特征提取。
所述特征提取具体为：
（1）获取像素坐标矩阵A的行数和列数，分别赋值给变量r和c；
（2）计算样本矩阵A的均值，将结果存储在变量meanVec中；
（3）计算协方差矩阵的转置covMatT：通过将样本矩阵A减去均值向量meanVec并进行矩阵乘法得到；
（4）使用eigs函数计算covMatT的前k个特征值和对应的特征向量，将结果分别存储在变量V和D中；所述eigs函数的主要参数包括输入矩阵、所需计算的特征值数量（k值）以及可选参数；通过调用eigs函数，可以得到指定数量（k个）的特征值和对应的特征向量；
（5）将协方差矩阵(covMatT)'的特征向量V乘以样本矩阵A减去均值向量meanVec的转置Z，得到变换后的特征向量；
（6）对变换后的特征向量进行归一化处理，使其成为单位特征向量；
（7）使用线性变换（投影）将样本矩阵A减去均值向量meanVec的转置Z乘以变换后的特征向量V，得到降维至k维的样本特征向量组成的矩阵pcaA。
所述跟踪当前帧，并添加关键帧具体为：
（1）提取矩阵区域的特征点作为初始特征点集合；
（2）特征点跟踪：将初始特征点作为输入，使用cv2.calcOpticalFlowPyrLK()函数在下一帧中跟踪这些特征点。该函数会返回跟踪到的特征点的新位置；
（3）更新目标区域：根据跟踪到的新特征点的位置，更新目标区域的位置和大小；
（4）关键帧判断：根据特征点数量的变化、目标区域的变化来判断当前帧是否满足添加关键帧的条件；
（5）添加关键帧：如果当前帧被判断为关键帧，将当前帧作为关键帧，并将当前帧的特征点集合添加到关键帧集合中；
（6）局部优化：对关键帧集合进行局部优化，通过优化特征点的位置和相机位姿，提高整体的几何一致性；
重复执行步骤（2）到步骤（6），以连续跟踪目标并添加关键帧，将得到具有不同特征点的集合。
所述图像拼接完成后使用PIL去除正射影像黑边，此时图中坐标丢失；
所述合成正射影像的tif图片存在黑色边界背景，切片完成后在地图显示也会出现黑色边界，需要进行处理，具体为：
（1）使用python的PIL库设置允许处理大尺寸图像并加载被截断的图像；
（2）打开了指定路径的图像文件，并将其转换为RGBA模式；
（3）获取图像的宽度L和高度H，获取图像左上角像素的颜色color_0，然后对图像的每个像素进行遍历，获取像素的颜色color_1；如果color_1与color_0相同，则将color_1的透明度设置为0，并将修改后的颜色重新赋值给该像素；
（4）保存处理后的图像到指定路径，并返回输出路径。
所述还包括使用gdal库还原坐标，将带有完整坐标正射影像的仿射矩阵、地图投影信息写入去除黑边后的正射影像；所述除去黑边后导致丢失地理坐标信息的正射影像，通过以下处理还原完整坐标正射影像的仿射矩阵、地图投影信息，具体为：
（1）根据去除黑边后返回的输出路径打开栅格数据集文件；
（2）读取栅格数据集的数组表示（img_array）；
（3）根据数组的数据类型，确定输出文件的数据类型（datatype）；
（4）确定图像的波段数（img_bands）以及图像的高度和宽度（im_height和im_width）；
（5）创建输出文件驱动（driver）和数据集（dataset），并设置输出文件的大小、波段数、数据类型信息；
（6）设置输出文件的仿射变换参数（img_transf）和投影信息（img_proj）（此信息通过读取原始正射影像文件中带有的完整坐标正射影像的仿射矩阵、地图投影信息获得）；
（7）将图像数据写入输出文件的相应波段；
（8）返回输出文件路径（outPath）；
所述除去黑边后导致丢失地理坐标信息的正射影像，通过以上处理还原完整坐标正射影像的仿射矩阵、地图投影信息。
所述正射影像处理完成后使用gdal进行切片，包括以下步骤：
（1）将原始影像数据集重投影到Web墨卡托投影坐标系，并将结果保存到临时文件dataset中即打开指定路径的TIFF文件，并将其重投影到Web墨卡托投影坐标系，具体为：
S1创建一个临时文件路径tmpFile 用于存储GTiff格式的输出数据集；
S2创建一个SpatialReference对象，用于存储空间参考信息；
S3使用spatialReference的ImportFromEPSG方法导入Web墨卡托投影坐标系的EPSG代码（WEB_MERCATOR_EPSG_CODE）；
S4使用gdal.AutoCreateWarpedVRT方法创建一个虚拟数据集webMercatorDs，将原始影像数据集（srcDs）重投影到Web墨卡托投影坐标系；
S5使用gdal.GetDriverByName("GTiff").CreateCopy方法创建一个与webMercatorDs相同的GTiff格式的输出数据集，并将其保存到tmpFile路径中；
S6删除原始影像数据集（srcDs）和虚拟数据集（webMercatorDs）；
（2）获取重投影后影像的地理范围和像素分辨率，并计算经纬度范围，获取原始影像在WGS84经纬度坐标系下的地理范围；具体为：
S1使用dataset.GetGeoTransform()方法获取原始影像的地理变换参数（geoTransform），包括地理坐标系的起始经度、像素的宽度、像素的高度、地理坐标系的起始纬度、像素的宽度和像素的高度；
S2获取原始影像的像素分辨率，通过dataset.getRasterXSize()和dataset.getRasterYSize()方法获取原始影像的宽度和高度；
S3 根据地理变换参数计算原始影像的经纬度范围，经度范围由起始经度（lngMin）和宽度乘以像素宽度（xSize * geoTransform[1]）以及高度乘以像素高度（ySize* geoTransform[2]）得到；纬度范围由起始纬度（latMax）和宽度乘以像素宽度（xSize *geoTransform[4]）以及高度乘以像素高度（ySize * geoTransform[5]）得到；
S4 得到的结果是原始影像的经纬度范围（lngMin, latMax, lngMax, latMin）；
（3）将原始影像的地理范围从Web墨卡托投影坐标系转换为WGS84经纬度坐标系，以便切片完成的瓦片可在天地图上展示，具体为：
S1使用EPSG代码解码获取源坐标参考系统（sourceCRS）和目标坐标参考系统（targetCRS）；其中，源坐标参考系统为Web墨卡托投影坐标系，目标坐标参考系统为WGS84经纬度坐标系；
S2 使用源坐标参考系统和目标坐标参考系统，通过CRS.findMathTransform方法获取坐标转换对象（transform）；
S3 创建Coordinate对象，用于存储经纬度的最大值和最小值；
S4 使用JTS.transform方法，将原始影像的地理范围的最大经纬度（lngMax,latMax）和最小经纬度（lngMin, latMin）从源坐标参考系统转换为目标坐标参考系统，并存储到新的Coordinate对象（lngLatMax, lngLatMin）中；
（4）根据指定的缩放级别，计算切片的行列号，具体为：
S1计算切片的x坐标：
int x = (int) Math.floor((lng + 180) / 360 * (1 << zoom));
公式中，lng为经度，zoom为缩放级别；首先将经度lng转换为切片内的相对经度，然后乘以切片数量(1 << zoom)，再向下取整得到切片的x坐标；
S2计算切片的y坐标：
int y = (int) Math.floor((1 - Math.log(Math.tan(Math.toRadians(lat))+ 1 / Math.cos(Math.toRadians(lat))) / Math.PI) / 2 * (1 << zoom));
公式中，lat为纬度，zoom为缩放级别；首先将纬度lat转换为弧度，然后使用数学公式计算切片内的相对纬度，再乘以切片数量(1 << zoom)，最后向下取整得到切片的y坐标；
S3对坐标进行边界处理：如果x坐标小于0，则将其设置为0；如果x坐标大于等于切片数量(1 << zoom)，则将其设置为切片数量减1；同样的，如果y坐标小于0或大于等于切片数量(1 << zoom)，也进行相应的边界处理；
最后得到行列号用于计算经纬度范围；
（5）根据切片的行列号和缩放级别，计算切片的经纬度范围，具体为：
S1接收两个参数：切片的列号x和缩放级别z；
S2使用公式 x / Math.pow(2.0, z) * 360.0 - 180 将切片的列号x转换为对应的经度值；
- Math.pow(2.0, z) 表示2的z次方，即切片在当前缩放级别下的总列数；
- x / Math.pow(2.0, z) 将切片的列号x转换为相对于总列数的比例；
- x / Math.pow(2.0, z) * 360.0 将相对比例转换为经度范围，取值范围为0到360；
- x / Math.pow(2.0, z) * 360.0 - 180 将经度范围从0到360转换为-180到180，得到最终的经度值；
S3返回转换后的经度值；
（6）根据切片的地理范围、像素分辨率以及原始影像的地理范围，计算出原始影像在切片内的像素大小和偏移像素坐标，用于后续的切片生成和裁剪操作，具体为：
根据原始影像的地理范围与切片的交集，确定切片的像素分辨率和范围。
S1创建了切片的地理范围（tileBound）和原始影像的地理范围（imageBound），都使用了WGS84经纬度坐标系作为参考坐标系；
S2使用tileBound和imageBound计算它们的交集（intersect），得到切片和原始影像共同的地理范围；
S3计算切片的像素分辨率：
- 切片的东西方向像素分辨率（dstWePixelResolution）：切片的经度范围除以切片的像素宽度（256）；
- 切片的南北方向像素分辨率（dstNsPixelResolution）：切片的纬度范围除以切片的像素高度（256）；
S4计算交集的像素信息：
- offsetX和offsetY：计算在切片地理范围内，原始影像的起始点像素坐标。通过将交集起始点的经度和纬度减去原始影像的起始经度和纬度，并除以原始影像的像素分辨率得到；
- blockXSize和blockYSize：计算在切片地理范围内，原始影像的像素大小。通过将交集的经度和纬度范围除以原始影像的像素分辨率得到；
- imageXBuf和imageYBuf：计算原始影像在切片内的像素大小；通过将交集的经度和纬度范围除以切片的像素分辨率得到，并向上取整；
- imageOffsetX和imageOffsetY：计算原始影像在切片中的偏移像素坐标。通过将交集起始点的经度和纬度减去切片的起始经度和纬度，并除以切片的像素分辨率得到。如果计算结果小于0，则将其设为0；
S5得到的结果是切片的像素分辨率和原始影像在切片内的像素大小和偏移像素坐标；
（7）使用GDAL的ReadRaster方法读取原始影像在切片范围内的像素数据到对应的缓冲区数组，具体为：
S1获取原始影像数据集的第1、2、3个波段（Band）对象，分别赋值给inBand1、inBand2和inBand3变量；
S2创建用于存储波段数据的缓冲区数组：
- band1BuffData数组用于存储第1个波段的数据，数组大小为TILE_SIZE *TILE_SIZE * gdalconst.GDT_Int32；
- band2BuffData数组用于存储第2个波段的数据，数组大小同上；
- band3BuffData数组用于存储第3个波段的数据，数组大小同上；
S3使用inBand1的ReadRaster方法从原始影像的第1个波段中读取像素数据，并将数据存储到band1BuffData数组中。读取的范围是从原始影像的偏移像素坐标(offsetX,offsetY)开始，读取blockXSize * blockYSize大小的像素块，并将其写入到imageXBuf *imageYBuf大小的缓冲区中；
S4 同样的方式，使用inBand2的ReadRaster方法从原始影像的第2个波段中读取像素数据，并将数据存储到band2BuffData数组中；
S5同样的方式，使用inBand3的ReadRaster方法从原始影像的第3个波段中读取像素数据，并将数据存储到band3BuffData数组中；
（8）创建内存驱动器，并将图像数据写入内存中；通过创建了一个内存数据集（msmDS），并获取了该数据集的各个波段对象，以便后续的像素数据写入和处理操作，具体为：
S1 获取名为"MEM"的GDAL驱动器（Driver）对象，并将其赋值给memDriver变量；
S2使用memDriver的Create方法创建一个名为"msmDS"的内存数据集（Dataset），该数据集的大小为256x256像素，有4个波段；
S3获取msmDS数据集的第1、2、3个波段（Band）对象，分别赋值给dstBand1、dstBand2和dstBand3变量；
（9）为图像数据设置alpha通道，实现背景透明效果；通过将原始影像的各个波段数据和alpha波段数据写入到内存数据集msmDS中，并将生成的切片文件保存到指定目录中，以便后续在Web地图上展示，具体为：
S1获取msmDS数据集的第4个波段（alphaBand）对象，并将其赋值给alphaBand变量；
S2创建一个用于存储alpha波段数据的缓冲区数组alphaData，数组大小为256x256xgdalconst.GDT_Int32；
S3使用循环遍历alphaData数组的每个元素，如果对应的band1BuffData数组的元素大于0，则将alphaData数组的元素设置为255，否则保持为0；
S4使用dstBand1、dstBand2、dstBand3和alphaBand的WriteRaster方法将相应的波段数据写入到msmDS数据集中；写入的范围是从imageOffsetX、imageOffsetY开始，写入imageXBuf x imageYBuf大小的像素块；
S5创建一个目录对象dir，用于存储切片文件；
S6创建切片文件的路径pngPath，路径格式为outputDir/Z[zoom]/[row]/[col].png。

Claims (9)

1.一种基于无人机航拍遥感图像自动处理方法，其特征在于：包括以下步骤：
步骤一、通过网联无人机在飞行过程中就实时将拍摄图像自动上传至对象存储服务器；
步骤二、飞行任务完成后自动开始下载完整的任务遥感图像；
步骤三、启动图片处理程序，自动合成正射影像；
步骤四、按国家自然资源局地图2000的坐标参数要求自动调用切图程序，并将正射影像切割为地图可显示文件，根据任务名自动将文件上传至服务器，自动解压切图文件；
步骤五、以上步骤处理得到的地图切片全景图像，可供无人机巡检平台根据任务查看该任务瓦片，并加载到地图上。
2.根据权利要求1所述的一种基于无人机航拍遥感图像自动处理方法，其特征在于：所述方法还包括使用无人机拍摄遥感图像结合自动化流程对遥感图像进行自动化处理；读取遥感图像同归对比图像中的参考物，将具有相同图像的像素移除，拼接不同的图像，合成正射影像；再将影像投影至web墨卡托，获取遥感图像经纬范围，计算影像像素分辨率，坐标转wgs84经纬度；将原始影影像地理范围与指定缩放级别指定行列号的切片交集，计算交集的像素信息，将切片数据写入文件。
3.根据权利要求2所述的一种基于无人机航拍遥感图像自动处理方法，其特征在于：所述图片处理程序自动合成正射影像具体包括：
（1）根据任务号从对象存储服务器下载相应任务遥感图像；
（2）对图像进行去畸变处理；
（3）特征提取；
（4）跟踪当前帧，并添加关键帧；
（5）环路检测：提取不同特征点的集合中的特征点和特征描述子；通过特征匹配算法对特征点进行匹配，并根据匹配结果进行几何一致性验证；最后，根据验证结果判断是否存在环路连接，如果存在，则进行局部优化处理；
所述局部优化通过优化相机位姿和特征点的三维坐标，调整关键帧之间的对齐关系，使用非线性优化算法Levenberg-Marquardt来最小化重投影误差，使得特征点在各个关键帧中的投影位置更加一致；
所述经过环路检测和局部优化后的关键帧集合，通过特征点的三维坐标和相机位姿，进行平面拟合，找到能够最好拟合关键帧上特征点的平面，使用最小二乘法算法，得到平面的参数方程；最后，将图像序列中的特征点投影到拟合的平面上，获得更好的对齐效果和视觉连续性；
（6）图像变换：根据相机位姿和平面参数，对图像进行透视变换，将图像投影到平面上，得到正交视图；
（7）图像拼接、融合：根据相机位姿和特征点的三维坐标，进行图像拼接，将多个图像拼接成全景图；所述进行图像融合：通过加权平均、多重曝光，将拼接后的图像进行平滑过渡和细节保留；
通过以上处理得到二维正射影像用于后续处理。
4.根据权利要求3所述的一种基于无人机航拍遥感图像自动处理方法，其特征在于：所述对图像进行去畸变处理具体为：
（1）创建目标点集合dst；
（2）从相机矩阵 cameraMatrix 和畸变系数 distortionCoeff 中获取相机参数：焦距（fx、fy），光心坐标（ux、uy），径向畸变系数（k1、k2、k3）和切向畸变系数（p1、p2）；
（3）遍历输入的源点坐标数组 src ；
（4）对于每个源点坐标，首先进行坐标转换，将其转换为相对于光心的归一化坐标：xDistortion = (p.x - ux) / fx，yDistortion = (p.y - uy) / fy；
（5）使用迭代的方式进行去畸变的求解。通过设定初始值（x0、y0），根据畸变模型的公式进行迭代计算 ：
畸变公式为：r2 = xDistortion*xDistortion + yDistortion*yDistortion
distRadialA = 1 / (1. + k1 * r2 + k2 * r2 * r2 + k3 * r2 * r2 * r2)
distRadialB = 1. + k4 * r2 + k5 * r2 * r2 + k6 * r2 * r2 * r2
deltaX = 2. * p1 * xDistortion * yDistortion + p2 * (r2 + 2. *xDistortion * xDistortion)
deltaY = p1 * (r2 + 2. * yDistortion * yDistortion) + 2. * p2 *xDistortion * yDistortion
xCorrected = (x0 - deltaX)* distRadialA * distRadialB
yCorrected = (y0 - deltaY)* distRadialA * distRadialB
xDistortion = xCorrected
yDistortion = yCorrected
其中，xDistortion和yDistortion表示归一化坐标系下的畸变坐标，xCorrected和yCorrected表示去畸变后的坐标，ux和uy表示相机光心坐标，fx和fy表示相机的焦距，k1、k2、k3、k4、k5、k6、p1和p2是相机的畸变系数，这些参数通过cameraMatrix和distortionCoeff传递给函数，最终得到xDistortion和yDistortion畸变坐标，用于下一步坐标变换；
（6）进行坐标变换，将去畸变后的坐标转换为相机坐标系下的像素坐标：xCorrected =xCorrected * fx + ux，yCorrected = yCorrected * fy + uy；
（7）将去畸变后的坐标（xCorrected、yCorrected）添加到目标点坐标数组 dst 中；
（8）循环结束后， dst 中将存储了所有源点坐标经过去畸变处理后的结果；
通过以上处理得到了去畸变后的像素坐标，用于特征提取。
5.根据权利要求3所述的一种基于无人机航拍遥感图像自动处理方法，其特征在于：所述特征提取具体为：
（1）获取像素坐标矩阵A的行数和列数，分别赋值给变量r和c；
（2）计算样本矩阵A的均值，将结果存储在变量meanVec中；
（3）计算协方差矩阵的转置covMatT：通过将样本矩阵A减去均值向量meanVec并进行矩阵乘法得到；
（4）使用eigs函数计算covMatT的前k个特征值和对应的特征向量，将结果分别存储在变量V和D中；所述eigs函数的主要参数包括输入矩阵、所需计算的特征值数量（k值）以及可选参数；通过调用eigs函数，可以得到指定数量（k个）的特征值和对应的特征向量；
（5）将协方差矩阵(covMatT)'的特征向量V乘以样本矩阵A减去均值向量meanVec的转置Z，得到变换后的特征向量；
（6）对变换后的特征向量进行归一化处理，使其成为单位特征向量；
（7）使用线性变换（投影）将样本矩阵A减去均值向量meanVec的转置Z乘以变换后的特征向量V，得到降维至k维的样本特征向量组成的矩阵pcaA。
6.根据权利要求3所述的一种基于无人机航拍遥感图像自动处理方法，其特征在于：所述跟踪当前帧，并添加关键帧具体为：
（1）提取矩阵区域的特征点作为初始特征点集合；
（2）特征点跟踪：将初始特征点作为输入，使用cv2.calcOpticalFlowPyrLK()函数在下一帧中跟踪这些特征点。该函数会返回跟踪到的特征点的新位置；
（3）更新目标区域：根据跟踪到的新特征点的位置，更新目标区域的位置和大小；
（4）关键帧判断：根据特征点数量的变化、目标区域的变化来判断当前帧是否满足添加关键帧的条件；
（5）添加关键帧：如果当前帧被判断为关键帧，将当前帧作为关键帧，并将当前帧的特征点集合添加到关键帧集合中；
（6）局部优化：对关键帧集合进行局部优化，通过优化特征点的位置和相机位姿，提高整体的几何一致性；
重复执行步骤（2）到步骤（6），以连续跟踪目标并添加关键帧，将得到具有不同特征点的集合。
7.根据权利要求3所述的一种基于无人机航拍遥感图像自动处理方法，其特征在于：所述图像拼接完成后使用PIL去除正射影像黑边，此时图中坐标丢失；
所述合成正射影像的tif图片存在黑色边界背景，切片完成后在地图显示也会出现黑色边界，需要进行处理，具体为：
（1）使用python的PIL库设置允许处理大尺寸图像并加载被截断的图像；
（2）打开了指定路径的图像文件，并将其转换为RGBA模式；
（3）获取图像的宽度L和高度H，获取图像左上角像素的颜色color_0，然后对图像的每个像素进行遍历，获取像素的颜色color_1；如果color_1与color_0相同，则将color_1的透明度设置为0，并将修改后的颜色重新赋值给该像素；
（4）保存处理后的图像到指定路径，并返回输出路径。
8.根据权利要求7所述的一种基于无人机航拍遥感图像自动处理方法，其特征在于：所述还包括使用gdal库还原坐标，将带有完整坐标正射影像的仿射矩阵、地图投影信息写入去除黑边后的正射影像；所述除去黑边后导致丢失地理坐标信息的正射影像，通过以下处理还原完整坐标正射影像的仿射矩阵、地图投影信息，具体为：
（1）根据去除黑边后返回的输出路径打开栅格数据集文件；
（2）读取栅格数据集的数组表示（img_array）；
（3）根据数组的数据类型，确定输出文件的数据类型（datatype）；
（4）确定图像的波段数（img_bands）以及图像的高度和宽度（im_height和im_width）；
（5）创建输出文件驱动（driver）和数据集（dataset），并设置输出文件的大小、波段数、数据类型信息；
（6）设置输出文件的仿射变换参数（img_transf）和投影信息（img_proj）（此信息通过读取原始正射影像文件中带有的完整坐标正射影像的仿射矩阵、地图投影信息获得）；
（7）将图像数据写入输出文件的相应波段；
（8）返回输出文件路径（outPath）；
所述除去黑边后导致丢失地理坐标信息的正射影像，通过以上处理还原完整坐标正射影像的仿射矩阵、地图投影信息。
9.根据权利要求8所述的一种基于无人机航拍遥感图像自动处理方法，其特征在于：所述正射影像处理完成后使用gdal进行切片，包括以下步骤：
（1）将原始影像数据集重投影到Web墨卡托投影坐标系，并将结果保存到临时文件dataset中即打开指定路径的TIFF文件，并将其重投影到Web墨卡托投影坐标系，具体为：
S1创建一个临时文件路径tmpFile 用于存储GTiff格式的输出数据集；
S2创建一个SpatialReference对象，用于存储空间参考信息；
S3使用spatialReference的ImportFromEPSG方法导入Web墨卡托投影坐标系的EPSG代码（WEB_MERCATOR_EPSG_CODE）；
S4使用gdal.AutoCreateWarpedVRT方法创建一个虚拟数据集webMercatorDs，将原始影像数据集（srcDs）重投影到Web墨卡托投影坐标系；
S5使用gdal.GetDriverByName("GTiff").CreateCopy方法创建一个与webMercatorDs相同的GTiff格式的输出数据集，并将其保存到tmpFile路径中；
S6删除原始影像数据集（srcDs）和虚拟数据集（webMercatorDs）；
（2）获取重投影后影像的地理范围和像素分辨率，并计算经纬度范围，获取原始影像在WGS84经纬度坐标系下的地理范围；具体为：
S1使用dataset.GetGeoTransform()方法获取原始影像的地理变换参数（geoTransform），包括地理坐标系的起始经度、像素的宽度、像素的高度、地理坐标系的起始纬度、像素的宽度和像素的高度；
S2获取原始影像的像素分辨率，通过dataset.getRasterXSize()和dataset.getRasterYSize()方法获取原始影像的宽度和高度；
S3 根据地理变换参数计算原始影像的经纬度范围，经度范围由起始经度（lngMin）和宽度乘以像素宽度（xSize * geoTransform[1]）以及高度乘以像素高度（ySize *geoTransform[2]）得到；纬度范围由起始纬度（latMax）和宽度乘以像素宽度（xSize *geoTransform[4]）以及高度乘以像素高度（ySize * geoTransform[5]）得到；
S4 得到的结果是原始影像的经纬度范围（lngMin, latMax, lngMax, latMin）；
（3）将原始影像的地理范围从Web墨卡托投影坐标系转换为WGS84经纬度坐标系，以便切片完成的瓦片可在天地图上展示，具体为：
S1使用EPSG代码解码获取源坐标参考系统（sourceCRS）和目标坐标参考系统（targetCRS）；其中，源坐标参考系统为Web墨卡托投影坐标系，目标坐标参考系统为WGS84经纬度坐标系；
S2 使用源坐标参考系统和目标坐标参考系统，通过CRS.findMathTransform方法获取坐标转换对象（transform）；
S3 创建Coordinate对象，用于存储经纬度的最大值和最小值；
S4 使用JTS.transform方法，将原始影像的地理范围的最大经纬度（lngMax, latMax）和最小经纬度（lngMin, latMin）从源坐标参考系统转换为目标坐标参考系统，并存储到新的Coordinate对象（lngLatMax, lngLatMin）中；
（4）根据指定的缩放级别，计算切片的行列号，具体为：
S1计算切片的x坐标：
int x = (int) Math.floor((lng + 180) / 360 * (1 << zoom));
公式中，lng为经度，zoom为缩放级别；首先将经度lng转换为切片内的相对经度，然后乘以切片数量(1 << zoom)，再向下取整得到切片的x坐标；
S2计算切片的y坐标：
int y = (int) Math.floor((1 - Math.log(Math.tan(Math.toRadians(lat)) + 1/ Math.cos(Math.toRadians(lat))) / Math.PI) / 2 * (1 << zoom));
公式中，lat为纬度，zoom为缩放级别；首先将纬度lat转换为弧度，然后使用数学公式计算切片内的相对纬度，再乘以切片数量(1 << zoom)，最后向下取整得到切片的y坐标；
S3对坐标进行边界处理：如果x坐标小于0，则将其设置为0；如果x坐标大于等于切片数量(1 << zoom)，则将其设置为切片数量减1；同样的，如果y坐标小于0或大于等于切片数量(1 << zoom)，也进行相应的边界处理；
最后得到行列号用于计算经纬度范围；
（5）根据切片的行列号和缩放级别，计算切片的经纬度范围，具体为：
S1接收两个参数：切片的列号x和缩放级别z；
S2使用公式 x / Math.pow(2.0, z) * 360.0 - 180 将切片的列号x转换为对应的经度值；
- Math.pow(2.0, z) 表示2的z次方，即切片在当前缩放级别下的总列数；
- x / Math.pow(2.0, z) 将切片的列号x转换为相对于总列数的比例；
- x / Math.pow(2.0, z) * 360.0 将相对比例转换为经度范围，取值范围为0到360；
- x / Math.pow(2.0, z) * 360.0 - 180 将经度范围从0到360转换为-180到180，得到最终的经度值；
S3返回转换后的经度值；
（6）根据切片的地理范围、像素分辨率以及原始影像的地理范围，计算出原始影像在切片内的像素大小和偏移像素坐标，用于后续的切片生成和裁剪操作，具体为：
根据原始影像的地理范围与切片的交集，确定切片的像素分辨率和范围；
S1创建了切片的地理范围（tileBound）和原始影像的地理范围（imageBound），都使用了WGS84经纬度坐标系作为参考坐标系；
S2使用tileBound和imageBound计算它们的交集（intersect），得到切片和原始影像共同的地理范围；
S3计算切片的像素分辨率：
- 切片的东西方向像素分辨率（dstWePixelResolution）：切片的经度范围除以切片的像素宽度（256）；
- 切片的南北方向像素分辨率（dstNsPixelResolution）：切片的纬度范围除以切片的像素高度（256）；
S4计算交集的像素信息：
- offsetX和offsetY：计算在切片地理范围内，原始影像的起始点像素坐标。通过将交集起始点的经度和纬度减去原始影像的起始经度和纬度，并除以原始影像的像素分辨率得到；
- blockXSize和blockYSize：计算在切片地理范围内，原始影像的像素大小。通过将交集的经度和纬度范围除以原始影像的像素分辨率得到；
- imageXBuf和imageYBuf：计算原始影像在切片内的像素大小；通过将交集的经度和纬度范围除以切片的像素分辨率得到，并向上取整；
- imageOffsetX和imageOffsetY：计算原始影像在切片中的偏移像素坐标。通过将交集起始点的经度和纬度减去切片的起始经度和纬度，并除以切片的像素分辨率得到。如果计算结果小于0，则将其设为0；
S5得到的结果是切片的像素分辨率和原始影像在切片内的像素大小和偏移像素坐标；
（7）使用GDAL的ReadRaster方法读取原始影像在切片范围内的像素数据到对应的缓冲区数组，具体为：
S1获取原始影像数据集的第1、2、3个波段（Band）对象，分别赋值给inBand1、inBand2和inBand3变量；
S2创建用于存储波段数据的缓冲区数组：
- band1BuffData数组用于存储第1个波段的数据，数组大小为TILE_SIZE * TILE_SIZE * gdalconst.GDT_Int32；
- band2BuffData数组用于存储第2个波段的数据，数组大小同上；
- band3BuffData数组用于存储第3个波段的数据，数组大小同上；
S3使用inBand1的ReadRaster方法从原始影像的第1个波段中读取像素数据，并将数据存储到band1BuffData数组中。读取的范围是从原始影像的偏移像素坐标(offsetX,offsetY)开始，读取blockXSize * blockYSize大小的像素块，并将其写入到imageXBuf *imageYBuf大小的缓冲区中；
S4 同样的方式，使用inBand2的ReadRaster方法从原始影像的第2个波段中读取像素数据，并将数据存储到band2BuffData数组中；
S5同样的方式，使用inBand3的ReadRaster方法从原始影像的第3个波段中读取像素数据，并将数据存储到band3BuffData数组中；
（8）创建内存驱动器，并将图像数据写入内存中；通过创建了一个内存数据集（msmDS），并获取了该数据集的各个波段对象，以便后续的像素数据写入和处理操作，具体为：
S1 获取名为"MEM"的GDAL驱动器（Driver）对象，并将其赋值给memDriver变量；
S2使用memDriver的Create方法创建一个名为"msmDS"的内存数据集（Dataset），该数据集的大小为256x256像素，有4个波段；
S3获取msmDS数据集的第1、2、3个波段（Band）对象，分别赋值给dstBand1、dstBand2和dstBand3变量；
（9）为图像数据设置alpha通道，实现背景透明效果；通过将原始影像的各个波段数据和alpha波段数据写入到内存数据集msmDS中，并将生成的切片文件保存到指定目录中，以便后续在Web地图上展示，具体为：
S1获取msmDS数据集的第4个波段（alphaBand）对象，并将其赋值给alphaBand变量；
S2创建一个用于存储alpha波段数据的缓冲区数组alphaData，数组大小为256x256xgdalconst.GDT_Int32；
S3使用循环遍历alphaData数组的每个元素，如果对应的band1BuffData数组的元素大于0，则将alphaData数组的元素设置为255，否则保持为0；
S4使用dstBand1、dstBand2、dstBand3和alphaBand的WriteRaster方法将相应的波段数据写入到msmDS数据集中；写入的范围是从imageOffsetX、imageOffsetY开始，写入imageXBuf x imageYBuf大小的像素块；
S5创建一个目录对象dir，用于存储切片文件；
S6创建切片文件的路径pngPath，路径格式为outputDir/Z[zoom]/[row]/[col].png。
CN202311449087.6A 2023-11-02 2023-11-02 一种基于无人机航拍遥感图像自动处理方法 Pending CN117495740A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202311449087.6A CN117495740A (zh) 2023-11-02 2023-11-02 一种基于无人机航拍遥感图像自动处理方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202311449087.6A CN117495740A (zh) 2023-11-02 2023-11-02 一种基于无人机航拍遥感图像自动处理方法

Publications (1)

Publication Number Publication Date
CN117495740A true CN117495740A (zh) 2024-02-02

Family

ID=89672028

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202311449087.6A Pending CN117495740A (zh) 2023-11-02 2023-11-02 一种基于无人机航拍遥感图像自动处理方法

Country Status (1)

Country Link
CN (1) CN117495740A (zh)

Similar Documents

Publication Publication Date Title
CN108564527B (zh) 基于神经网络的全景图内容补全和修复的方法及装置
CN110163064B (zh) 一种道路标志物的识别方法、装置及存储介质
CN108965742B (zh) 异形屏显示方法、装置、电子设备及计算机可读存储介质
Li et al. Integrated shadow removal based on photogrammetry and image analysis
EP2057585A2 (en) Mosaic oblique images and methods of making and using same
CN111582022A (zh) 一种移动视频与地理场景的融合方法、系统及电子设备
Peña-Villasenín et al. 3-D modeling of historic façades using SFM photogrammetry metric documentation of different building types of a historic center
CN113160053B (zh) 一种基于位姿信息的水下视频图像复原与拼接方法
CN111105347B (zh) 一种生成带深度信息的全景图的方法、装置及存储介质
CN112288628B (zh) 基于光流跟踪和抽帧映射的航拍图像拼接加速方法及系统
CN106504196A (zh) 一种基于空间球面的全景视频拼接方法及设备
Alsadik et al. Efficient use of video for 3D modelling of cultural heritage objects
Agisoft Metashape python reference
CN112055192A (zh) 图像处理方法、图像处理装置、电子设备及存储介质
CN109829421B (zh) 车辆检测的方法、装置及计算机可读存储介质
CN113066173B (zh) 三维模型构建方法、装置和电子设备
CN112529498B (zh) 一种仓储物流管理方法及系统
CN116760937B (zh) 一种基于多机位的视频拼接方法、装置、设备及存储介质
CN114693782A (zh) 用于确定三维场景模型坐标系和物理坐标系的转换关系的方法及装置
CN114898068B (zh) 三维建模方法、装置、设备及存储介质
CN117495740A (zh) 一种基于无人机航拍遥感图像自动处理方法
US10635925B2 (en) Method and system for display the data from the video camera
CN116012227A (zh) 图像处理方法、装置、存储介质及处理器
CN115222602A (zh) 图像拼接方法、装置、设备及存储介质
CN114663284A (zh) 红外热成像全景图像处理方法、系统及存储介质

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