基于GPU并行计算的无人机影像拉花区域自动探测方法

摘要

本发明公开了一种基于GPU并行计算的无人机影像拉花区域自动探测方法，包括CPU端加载原始无人机影像、影像内方位元素、外方位元素和数字高程模型DEM，并根据中心投影构像方程计算正射纠正后正射影像的大小和范围；CPU端创建GPU线程格网、开辟相应的显存，并拷入数据；GPU端按线程格网从纠正后影像像素反算其对应原始影像上的像素坐标；GPU端进行拉花像素判断；GPU端按线程格网对探测结果二值图像进行图像腐蚀和图像膨胀处理；CPU端将拉花检测后的二值图像结果从GPU全局存储器拷到内存，并将二值图像矢量化，获得拉花变形区域矢量数据等步骤。其显著效果是：实现了无人机单片正射影像中拉花变形区域探测的自动化，极大提高了无人机单片正射影像拉花区域探测效率。

说明书

技术领域

本发明涉及到无人机航摄影像处理技术领域，具体地说，是一种基于GPU并行计算的无人机影像拉花区域自动探测方法。

背景技术

数字正射影像图(DOM)可以通过无人机航拍影像结合地形数据(数字高程模型DEM)根据中心投影构像方程利用反解法数字微分纠正原理逐像素计算地面点与原始影像上像素的关系，并进行灰度重采样得到的既有正确的位置信息又有丰富的纹理信息的影像图。无人机具有机动灵活、成本小、影像分辨率高等优势在数字正射影像制作更新工作中发挥着越来越的作用。而无人机影像正射纠正是正射影像制作中的重要一环，正射纠正是利用无人机原始影像及其内外方位元素、数字高程模型DEM根据中心投影构像方程共同消除原始无人机影像中各种畸变(如投影差)而得到一幅既有地理坐标位置信息又有纹理信息的新影像的过程。由于地形起伏和中心投影的成像方式和影像拍摄瞬间相机的姿态位置等原因，不能保证所有地面点都能在无人机影像上成像，如坡势较陡的山坡可能会被山顶所遮挡。无人机影像正射纠正是根据地面点位置计算原始影像上对应的像素，然后进行灰度重采样，在被遮挡的区域或是成像信息匮乏的区域重采样将会过于稠密或是重复采样，这样就会造成正射纠正后影像出现拉伸现象，如果拉伸过度就会造成纹理失真的现象，我们称之为“拉伸变形”，纹理失真的区域我们称之为“拉花变形区域”。拉伸变形造成的纹理失真直接影响着数字正射影像图的质量，特别是山地区域拉伸变形现象特别严重，直接影响着数字正射影像图制作的效率。

在正常的数字正射影像制作中，如果出现拉花变形，需要人工目视辨别和寻找，然后再通过修改DEM再纠正的方式进行处理。目前，无人机正射影像拉花区域自动探测的方法较少，已有方法效率低下无法满足大序列无人机影像拉花区域自动探测的需求，同时现有方法漏检率也非常高，无法保证拉花区域自动探测结果的质量。

发明内容

针对现有技术的不足，本发明的目的是提供一种基于GPU并行计算的无人机影像拉花区域自动探测方法，利用无人机原始影像及其内、外方位元素和测区数字高程模型DEM，基于GPU并行处理技术快速准确地检测出纠正后单片正射影像上的拉花变形区域，达到拉花变形自动、快速探测的目的。

为达到上述目的，本发明采用的技术方案如下：

一种基于GPU并行计算的无人机影像拉花区域自动探测方法，其关键在于按照以下步骤进行处理：

步骤1：CPU端加载原始无人机影像、影像内方位元素、外方位元素和数字高程模型DEM，根据中心投影构像方程计算正射纠正后正射影像的大小和范围，计算方法如下：

步骤1-1：将原始无人机影像、影像内方位元素、外方位元素和数字高程模型DEM等数据读入内存中，计算原始无人机影像四个角点(x₁,y₁),(x₂,y₂),(x₃,y₃),(x₄,y₄)对应地面点的坐标(X₁,Y₁),(X₂,Y₂),(X₃,Y₃),(X₄,Y₄)，计算公式为：

其中，(X,Y)为地面点坐标，(x,y)为像点坐标，Xs,Ys,Zs为影像外方位元素中的三个线元素，a₁,b₁,c₁,a₂,b₂,c₂,a₃,b₃,c₃为影像外方位元素中三个角元素计算得到的旋转矩阵中的9个参数，f为相机焦距，Z₀为测区平均高程值；

步骤1-2：按照公式计算地面点坐标(X₁,Y₁),(X₂,Y₂),(X₃,Y₃),(X₄,Y₄)在X、Y方向的极值X_min,X_max,Y_min,Y_max；

步骤1-3：根据四个极值计算得到正射纠正后正射影像的列数width_dst和行数height_dst，计算公式为：

其中，M是正射纠正后影像的分辨率；

步骤2：CPU端根据纠正后正射影像大小创建GPU线程格网，开辟相应的显存，并将相应数据从内存拷贝入对应显存；

步骤3：GPU端按线程格网从纠正后影像像素反算其对应原始影像上的像素坐标；

步骤4：GPU端按照线程格网统计当前像素与周边像素重叠次数，并判断其是否为拉花像素；

步骤5：GPU端按线程格网对探测结果二值图像进行图像腐蚀与图像膨胀处理；

步骤6：CPU端将拉花检测后的二值图像结果从全局存储器拷到内存，并将二值图像矢量化，获得拉花变形区域矢量数据。

进一步的，步骤3中所述原始影像上的像素坐标求解步骤为：

步骤3-1：GPU端所有线程同时执行，根据当前线程索引计算当前像素的地面点坐标；

步骤3-2：根据地面点坐标和数字高程模型DEM，采用双线性内插法内插出该点的高程值；

步骤3-3：计算当前地面点在原始影像上的像点坐标，并判断该坐标是否在原始影像内部；

步骤3-4：等待所有线程计算完毕后，将整幅正射影像反算得到的每个像素对应原始影像像点位置存储在全局存储器中。

进一步的，步骤4中拉花像素的判断步骤为：

步骤4-1：按照线程索引以当前像素为中心，建立一个适当大小的方形窗口；

步骤4-2：遍历和统计方形窗口中所有像素对应的原始影像上像点坐标和当前中心像素对应原始影像像点坐标的重叠次数；

步骤4-3：当重叠次数大于阈值时，在检测结果影像上将当前中心像素标记为拉花像素，否则标记为非拉花像素；

步骤4-4：等待所有线程计算完毕，并将拉花像素判定结果存储到全局存储器中。

进一步的，步骤5中对二值图像进行图像形态学处理的步骤为：

步骤5-1：按照线程索引以当前像素为中心，建立一个适当大小的方形腐蚀窗口，若当前中心像素灰度值为g1时进入步骤5-2；

步骤5-2：遍历腐蚀窗口中每一个像素，统计灰度值为g1的像素个数sum，若sum小于腐蚀窗口内像素总个数的一半时，将当前中心像素赋值为0；

步骤5-3：按照线程索引以当前像素为中心，建立一个适当大小的方形膨胀窗口；

步骤5-4：将当前像素对应的灰度值为与灰度阈值g2进行比较，若相等则遍历膨胀窗口中每一个像素，并将其对应的灰度值赋值为灰度阈值g2；

步骤5-5：等待所有线程计算完毕，并将计算结果存储到全局存储器中。

进一步的，所述灰度阈值g1＝g2＝255。

进一步的，步骤6中所述矢量数据的获得步骤为：

步骤6-1：CPU端将拉花检测后的二值图像结果从全局存储器拷贝到内存中，同时释放已经开辟的GPU存储器；

步骤6-2：对二值图像中拉花变形区域进行边缘检测及提取，并存储为矢量数据。

本方案中，首先CPU端读取原始无人机影像及其内、外方位元素、测区DEM，并计算正射纠正后对应正射影像的大小和范围；然后CPU端根据纠正后正射影像大小分配线程格网开辟多种GPU存储器，并将DEM等数据绑定在相应存储器；接着GPU端根据中心投影构像方程利用反解法计算正射影像上每个像素在原始影像上的像点坐标；再由GPU端统计正射影像上每个像素与周边像素所对应原始影像像点的重叠次数判断其是否为变形像素，生成变形区域二值图像，并进行图像腐蚀和膨胀处理；最后CPU端将最终拉花区域检测结果二值图像从GPU显存拷贝到内存，在内存中对二值图像进行矢量化处理，得到单片无人机正射影像上变形区域的矢量数据。

本发明的显著效果是：利用无人机原始影像及其内、外方位元素和测区数字高程模型DEM等数据，不仅实现了无人机单片正射影像中拉花变形区域探测的自动化，解决了传统人工目视辨别查找费时费力和人为遗漏的问题；还将GPU并行计算技术成功运用于无人机单片正射影像拉花区域自动探测中，探测效率是传统CPU串行算法的40倍左右，极大提高了无人机单片正射影像拉花区域探测效率。

附图说明

图1是本发明的方法流程图；

图2是无人机正射影像上检测出的拉花区域局部1；

图3是无人机正射影像上检测出的拉花区域局部2。

具体实施方式

下面结合附图对本发明的具体实施方式以及工作原理作进一步详细说明。

如图1所示，一种基于GPU并行计算的无人机影像拉花区域自动探测方法，具体步骤如下：

步骤1：CPU端加载原始无人机影像、影像内方位元素、外方位元素和数字高程模型DEM，根据中心投影构像方程计算正射纠正后正射影像的大小和范围，计算方法如下：

步骤1-1：将原始无人机影像、影像内外方位元素、数字高程模型DEM等数据读入到内存中，并根据公式(1)计算原始影像四角点(x₁,y₁),(x₂,y₂),(x₃,y₃),(x₄,y₄)对应地面点的坐标(X₁,Y₁),(X₂,Y₂),(X₃,Y₃),(X₄,Y₄)。

其中，(X,Y)为地面点坐标，(x,y)为像点坐标，Xs,Ys,Zs为影像外方位元素中的三个线元素，a₁,b₁,c₁,a₂,b₂,c₂,a₃,b₃,c₃为影像外方位元素中三个角元素计算得到的旋转矩阵中的9个参数，f为相机焦距，Z₀为该测区平均高程值；

步骤1-2：按照公式(2)计算四个角点对应的地面点坐标(X₁,Y₁),(X₂,Y₂),(X₃,Y₃),(X₄,Y₄)在X、Y方向的极值X_min,X_max,Y_min,Y_max，则(X_min,Y_max)为正射纠正后正射影像左上角的地理坐标；

步骤1-3：根据公式(3)，采用上述的四个极值X_min,X_max,Y_min,Y_max计算得到正射纠正后正射影像的列数width_dst和行数height_dst；

其中，M是正射纠正后影像的分辨率；

步骤2：CPU端根据纠正后正射影像大小创建GPU线程格网，开辟相应的显存，并将相应数据从内存拷贝入对应显存，具体过程如下所述：

GPU并行计算是依托GPU多线程之间的并行计算，而GPU线程则是通过线程格网进行组织和调度的。正射纠正是对影像每个像素进行重采样的过程，因此总的计算任务为width_dst*height_dst，也就是说至少需要一个包含width_dst*height_dst个线程的线程格网，线程格网中每个线程对应一个像素。线程格网是由一个个线程块组成，如果按照二维线程格网进行组织，线程块的大小为dimBlock(N,N)，即每个线程块中有N*N个线程，格网dimGrid的大小可以根据公式(4)计算得到：

GPU共有六种存储器：寄存器、局部存储器、共享存储器、全局存储器、常数存储器和纹理存储器，每种存储器的大小和访问速度均不相同，构成了GPU多层次性存储的特点。有效利用GPU多层次性存储的特点，可以进一步提高GPU的计算效率。线程格网中的每条线程拥有仅供自己访问的寄存器与局部存储器，同一个线程块中的线程能够读写块内的共享存储器，全局存储器可以供格网内的所有线程读写，而常数存储器和纹理存储器只能被格网内的所有线程读取，不能进行写操作。

因此针对不同大小的数据及其访问方式，选择不同的存储器可以有效提高数据的访问速度和处理能力，因此本方法中将供每条线程访问且数据量大的DEM存储在纹理存储器中，利用纹理存储器缓存加速的特点提高访问效率，将每条线程读取且数据量较小的原始影像内外方位元素存储在访问速度较快的常数存储器中，将大量的中间计算结果存放在全局存储器中。

步骤3：GPU端按线程格网从纠正后影像像素反算其对应原始影像上的像素坐标；求解步骤为：

步骤3-1：利用公式(5)根据当前线程索引计算当前像素的地面点坐标(X,Y)，其中M为正射影像的地面分辨率，(row,col)为当前线程索引对应像素所在的行列坐标，(X_min,Y_max)为正射纠正后正射影像的左上角点坐标，height_dst为正射纠正后影像的行数。

步骤3-2：根据地面点坐标(X,Y)和数字高程模型DEM，采用双线性内插法内插出该点的高程值Z：

首先，根据目标点的坐标(X,Y)和DEM坐标变换参数按照公式(6)得到目标点在DEM格网中的行列数(R_DEM,C_DEM)，式中(X_0DEM,Y_0DEM)是DEM数据左上角的地理坐标，Cellsize_DEM是DEM数据相邻两点之间的间隔。

然后，根据(R_DEM,C_DEM)得到DEM格网中最邻近的四个高程点Z_zsh,Z_ysh,Z_yx,Z_zx，然后根据公式(7)双线性内插出该点的高程值Z，其中dx,dy为该点距离高程点Z_zsh在x,y方向上距离。

步骤3-3：根据公式(8)计算地面点(X,Y,Z)在原始影像上的像点坐标(s,l)：其中，Xs,Ys,Zs为影像外方位元素中影像拍摄瞬间，摄影中心的坐标，a₁,b₁,c₁,a₂,b₂,c₂,a₃,b₃,c₃为影像外方位元素中三个角元素计算得到的旋转矩阵中9个参数，f为相机焦距。

并根据公式(9)判断坐标(s,l)是否在原始影像内部，如果在外部则将s,l的值其设置为-1，width_src和height_src分别为原始影像的列数和行数。

步骤3-4：等待所有线程计算完毕后，将整幅正射影像反算得到的每个像素对应原始影像像点位置存储在全局存储器中。

步骤4：GPU端按照线程格网统计当前像素与周边像素重叠次数，并判断其是否为拉花像素，在处理计算时GPU端根据线程格网进行并行计算，每条线程对应一个像素，所有线程并行执行：

步骤4-1：按照线程索引以当前像素(row,col)为中心，建立一个适当大小的方形窗口Win₁[20,20]，当前像素对应原始影像上的像点坐标为(S,L)；

步骤4-2：设重叠像素个数num＝0，遍历方形窗口Win₁[20,20]中每一个像素对应的原始影像上像点坐标(r_i,c_i)∈Win₁[20,20]与当前像素对应的原始影像上像点坐标(S,L)对比，若满足公式(10)则令num＝num+1，此实例中δ＝0.8；

步骤4-3：当重叠次数num大于阈值2时，将该像素视为拉花像素，在探测结果二值图像上将其标记为拉花像素(灰度值为255)，否则标记为非拉花像素(灰度值为0)；

步骤4-4：等待所有线程计算完毕，并将拉花像素判定结果存储到全局存储器中。

步骤5：GPU端按线程格网对探测结果二值图像进行图像腐蚀处理和图像膨胀处理：

首先进行图像腐蚀处理步骤：

步骤5-1：按照线程索引以当前像素(row,col)为中心，建立一个大小为3*3的方形腐蚀窗口Win₂[3,3]，若当前中心像素灰度值为g1＝255时进入步骤5-2；

步骤5-2：遍历腐蚀窗口中每一个像素(r_i,c_i)∈Win₂[3,3]，统计灰度值为g1的像素个数sum，若sum小于4，则该中心像素赋值为0；

然后进行图像膨胀处理步骤：

步骤5-3：按照线程索引以当前像素(row,col)为中心，建立一个大小为5*5的方形膨胀窗口Win₃[5,5]，当前像素对应的灰度值为gray₀；

步骤5-4：将当前像素对应的灰度值为与灰度阈值g2＝255进行比较，若gray₀＝g2则遍历膨胀窗口中每一个像素(r_i,c_i)∈Win₃[5,5]，并将其对应的灰度值赋值为灰度阈值g2＝255；

步骤5-5：等待所有线程计算完毕，并将计算结果存储到全局存储器中。

最后进入步骤6：CPU端将拉花检测后的二值图像结果从全局存储器拷到内存，并将二值图像矢量化，获得拉花变形区域矢量数据，获得步骤为：

步骤6-1：CPU端将拉花检测后的二值图像结果从全局存储器拷贝到内存中，同时释放已经开辟的GPU存储器；

步骤6-2：对二值图像中拉花变形区域进行边缘检测及提取，并存储为矢量数据。

本方案中，首先CPU端读取原始无人机影像及其内、外方位元素、测区DEM，并计算正射纠正后对应正射影像的大小和范围；然后CPU端根据纠正后正射影像大小分配线程格网开辟多种GPU存储器，并将DEM等数据绑定在相应存储器；接着GPU端根据中心投影构像方程利用反解法计算正射影像上每个像素在原始影像上的像点坐标；再由GPU端统计正射影像上每个像素与周边像素所对应原始影像像点的重叠次数，判断其是否为变形像素生成变形区域二值图像，并进行图像腐蚀和膨胀处理；最后CPU端将最终拉花区域检测结果二值图像从GPU显存拷贝到内存，在内存中对二值图像进行矢量化处理，得到单片无人机正射影像上变形区域的矢量数据，如图2与图3白色实线框中所示，从而不仅实现了无人机单片正射影像中拉花变形区域探测的自动化，极大提高了无人机单片正射影像拉花区域探测效率。