CN106326846B - 无人机图像的林木植株并行提取方法 - Google Patents

Abstract

本发明公开了一种无人机图像的林木植株并行提取方法，先采用GPU并行处理的方式，对无人机拍摄的林区图像进行基于尺度空间技术的二进制大型对象检测，得到检测点，再用CPU串行方式，将不含绿色的检测点删除，剩余的检测点即为提取结果。所采用的尺度空间技术可以在各种尺度下识别不同冠层大小的林木，精度较高，又由于采用GPU并行处理的方式进行对象检测，使得本提取方法具有较高的效率和速度，为林木植株的调查与研究提供了支持。

Description

无人机图像的林木植株并行提取方法

技术领域

本发明涉及林木植株的调查与研究技术领域，具体涉及一种无人机图像的林木植株并行提取方法。

背景技术

无人机被动光学遥感在时空分辨率、可获取性、成本方面较卫星被动光学遥感具有较大优势，当前越来越多的被应用于林木植株的调查与研究中。对于林木(主要是果园、人工林等)植株提取而言，主要是利用数字图像处理技术，从无人机拍摄的图像中识别出目标植株个体并进行统计或分析。

然而，目前的数字图像处理技术的精确度普遍较低，不同冠层的林木难以同时识别。同时，由于当前无人机图像分辨率普遍较高，使得数字图像处理需要进行大量的运算，导致算法耗时较长。

发明内容

针对现有技术的不足，本发明的目的在于提供一种无人机图像的林木植株并行提取方法，以提高图像提取的精度和速度。

为了实现上述目的，本发明采取的技术方案是：

一种无人机图像的林木植株并行提取方法，包括步骤：

获取无人机拍摄的林区图像；

采用GPU并行处理的方式，对林区图像进行基于尺度空间技术的二进制大型对象检测，得到初始检测点；

对初始检测点进行筛选，将不含绿色的检测点删除；

将筛选后的检测点作为提取结果。

与现有技术相比，本发明的有益效果在于：

本发明采用尺度空间技术对无人机获取的林区图像进行检测，可以在各种尺度下识别不同冠层大小的林木，由于加入了对象分析步骤，精度较高。同时，对关键步骤采用GPU并行处理的方式，使得本提取方法具有较高的效率和速度，为林木植株的调查与研究提供了支持。

附图说明

图1为本发明无人机图像的林木植株并行提取方法的流程示意图。

具体实施方式

下面结合具体实施方式对本发明作进一步的说明。

本发明无人机图像的林木植株并行提取方法，如图1所示，包括步骤：

步骤s101、获取无人机拍摄的林区图像；

步骤s102、采用GPU并行处理的方式，对林区图像进行基于尺度空间技术的二进制大型对象检测，得到检测点；

步骤s103、对初始检测点进行筛选，将不含绿色的检测点删除；

步骤s104、将筛选后的检测点作为提取结果。

利用PyCUDA平台，实现无人机图像的林木植株并行提取算法。CUDA(ComputeUnified Device Architecture)技术为NVIDIA公司开发的类C语言GPU编程平台。而PyCUDA是CUDA的Python语言封装，该封装提供对象自动清理功能，并通过结合Numpy、Scipy和GDAL等Python科学计算库，能够方便的开发复杂应用。

上述算法可以分为：1)输入；2)基于尺度空间技术的二进制大型对象(BinaryLarge Objects，BLOBS)检测；3)对象分析；4)输出，共四个部分。其中第3部分为图像处理，能够进行GPU并行处理，并且可以取得效率提升，是本发明的核心部分，采用并行算法编程，其余部分均采用串行编程。

其中，BLOBS检测部分包含3个核心步骤：

1)原始图像的多尺度LoG(Laplacian of Gaussian，LoG，拉普拉斯-高斯)滤波；

2)LoG滤波尺度空间图像的局部最大值滤波；

3)检测点的基于面积的组合筛选。

下面详细介绍上述步骤。

算法总共包含5个参数，分别为m_grey,min_sigma,max_sigma,num_sigma,threshold,overlap。(注：Sigma代表高斯滤波的尺度，其值越高模糊程度越高)其中，m_grey是输入的灰度图像，min_sigma,max_sigma,num_sigma分别对应拉普拉斯-高斯滤波的最小尺度、最大尺度和尺度数量，threshold代表BLOB识别灰度阈值，overlap代表BLOB筛选面积重叠阈值。

技术内容将按照算法处理步骤分别阐述：

(1)输入：通过GDAL读入无人机可见光RGB(Red,Green,Blue，红绿蓝)图像。利用公式1将RGB波段转换为灰度图像，Y表示灰度值。

Y＝0.2125*R+0.7154*G+0.0721*B (公式1)

(2)基于尺度空间技术的二进制大型对象(Binary Large Objects，BLOBS)检测：

准备步骤：根据参数min_sigma,max_sigma,num_sigma，分别计算num_sigma个等距的sigma值，记为sigma_list

1)原始图像的多尺度LoG滤波：针对原始图像m_grey，分别采用sigma_list中的各个sigma进行LoG滤波，得到LoG滤波尺度空间图像。其中，LoG滤波的公式为：

其中，x,y分别为图像中的x,y坐标，σ(sigma)为高斯滤波尺度。

实现流程为：

由于处理无人机图像需要较大的sigma(≤50)，而当前CUDA设备允许的一个块(Block)中的线程(Thread)数不能超过1024个(32*32)，当sigma为8.0是窗口大小已达到32。因此，二维滤波器无法满足需求，必须采用两个一维LoG滤波器进行组合，以便sigma较大时进行LoG滤波。组合方法为，针对某一sigma：

①对原始图像进行列方向LoG滤波，然后进行行方向高斯滤波，得到M_logx；

②对原始图像进行行方向LoG滤波，然后进行列方向高斯滤波，得到M_logy；

③M_log＝(M_logx+M_logy)*(-sigma²)。

最终针对sigma_list中的每个sigma进行滤波，并合并为三维图像数组。

注：Kernel函数采用常数内存(Constant Memory)存储滤波器参数；而对于每行或列采用一维数组索引共享内存(Shared Memory)，用以加快运算速度，即每个Block内的所有Thread(1024个)共用一份缓存的图像数据。由于滤波器尺寸较大，为避免边缘数据浪费，边界处理方式为反射(Reflect)，即M(x-i)＝M(x+i-1)。第③步由于不涉及滤波，直接在主显示内存(Global Memory)操作。

2)LoG滤波尺度空间图像的局部最大值滤波

实现流程为：

①窗口尺寸为3(半径为1)个像素，分别求取三维LoG滤波图像组在x,y,z方向的最大值，并生成最大值立方体；

②对比大值立方体和LoG立方体，对于满足公式的像素作为潜在的BLOBS。

Image_LoG(x,y,z)＝Image_max(x,y,z)并且Image_LoG(x,y,z)>Threshold (公式2)

③输出所有潜在BLOBS的x,y坐标，以及z轴对应的sigma值

注：对于每行或列采用一维数组索引共享内存(Shared Memory)，用以加快运算速度，即每个Block内的所有Thread(1024个)共用一份缓存的图像数据。而对于z轴，由于数量小于1024，与x轴共同分配1024个Threads，采用二维数组索引共享内存。为避免滤波器在边缘时窗口像元数量减少导致过多的像素被选取，边界处理方式为常数(Constant)，即M(x-i)＝0。

3)检测点的基于面积的组合筛选

输入潜在的BLOBS，为Python语言的List格式，每个单元包含x,y,sigma三元素。本步将对所有BLOBS进行两两组合，对比BLOBS的面积，如果某两个BLOB的重叠度小于阈值overlap，则删除面积较小的BLOB。本步骤利用GPU实现的难点在于，由于BLOBS长度未知，难以采用二维数组用Thread索引组合中的两个BLOB元素。本发明采用一维数组索引组合序号，然后通过组合序号反推BLOB元素。原理为代码1：

实现流程为：

①根据BLOBS数量N_BLOBS，利用公式3在CPU端计算两两组合的总数N_combination

N_combination＝N_BLOBS*(N_BLOBS-1)/2 (公式3)

②启动一维数组，为每个Block,Thread计算其对应的组合一维序号，采用迭代法计算其对应的两个BLOBS序号；

③利用公式4计算重叠度OL：

R₁＝sigma₁*2^1/2,R₂＝sigma₂*2^1/2

Distance＝((x₁–x₂)²+(y₁–y₂)²)^1/2

如果Distance>R₁+R₂

OL＝0

如果Distance≤|R₁-R₂|

OL＝1

否则：

Ratio₁＝(Distance²+R₁ ²–R₂ ²)/(2.0*Distance*R₁)

且Ratio₁＝max(Ratio₁,-1.0)，Ratio₁＝min(Ratio₁,1.0)

Acos₁＝arccos(Ratio₁)

Ratio₂＝(Distance²+R₂ ²–R₁ ²)/(2.0*Distance*R₂)

且Ratio₂＝max(Ratio₂,-1.0)，Ratio₂＝min(Ratio₂,1.0)

Acos₂＝arccos(Ratio₂)

A＝-Distance+R₂+R₁

B＝Distance-R₂+R₁

C＝Distance+R₂-R₁

D＝Distance+R₂+R₁

Area＝R₁ ²*Acos₁+R₂ ²*Acos₂-0.5*(|A*B*C*D|)^1/2

A_min＝min(R₁,R₂)

OL＝Area/(π*A_min ²) (公式4)

式中，sigma₁和sigma₂分别表示高斯滤波尺度1和高斯滤波尺度2，x₁、y₁表示高斯滤波尺度1下，潜在检测点1的坐标，x₂、y₂表示高斯滤波尺度2下，潜在检测点2的坐标，Distance表示欧氏距离；R₁、R₂、Ratio₁、Ratio₂、Acos₁、Acos₂、A、B、C、D、Area、A_min均表示中间变量。

④在CPU端，删除所有OL值小于overlap阈值的BLOB，剩余的BLOBS即为提取的BLOBS，接下来通过后续步骤分析其颜色，进一步筛选。

注：由于GPU芯片的对称性，只有当分配的Blocks*Threads值为2的整数次方时会发挥最佳性能，而GTX960m显卡允许的Blocks*Threads最大值为2097152(2²¹)，因此本发明将该值作为Kernel启动参数。

(3)对象分析：将原始图像从RGB色彩空间转换为CIE-Lab色彩空间，转换方式见公式9。根据上一步提取的BLOBS位置和面积，计算所有BLOB的a通道均值，删除值高于0的BLOBS(a值大于0表示不含绿色)，剩余点为检测结果

X＝R*0.4124+G*0.3576+B*0.1805 (公式9)

Y＝R*0.2126+G*0.7152+B*0.0722

Z＝R*0.0193+G*0.1192+B*0.9505

L＝116*f(Y)–16

a＝500*(f(X)-f(Y))

b＝200*(f(Y)-f(Z))

其中，当x>0.008856时，f(x)＝x^(1/3)，x表示X、Y、Z中的任意一个，当x<＝0.008856时，f(x)＝(7.787*x)+(16/116)。

(4)输出：输出为BLOBS列表，其中包含每个BLOBS的x,y坐标及其半径(注：半径为sigma)。

上列详细说明是针对本发明可行实施例的具体说明，该实施例并非用以限制本发明的专利范围，凡未脱离本发明所为的等效实施或变更，均应包含于本案的专利范围中。

Claims (2)

1.一种无人机图像的林木植株并行提取方法，其特征在于，包括步骤：

获取无人机拍摄的林区图像；

采用GPU并行处理的方式，对林区图像进行基于尺度空间技术的二进制大型对象检测，得到初始检测点；

对初始检测点进行筛选，将不含绿色的检测点删除；

将筛选后的检测点作为提取结果，

所述步骤获取无人机拍摄的林区图像包括：

通过GDAL读入无人机拍摄的林区的可见光RGB图像，利用下式将RGB图像转换为灰度图像，Y表示灰度值：

Y＝0.2125*R+0.7154*G+0.0721*B

所述步骤对林区图像进行基于尺度空间技术的二进制大型对象检测包括：

确定LoG滤波的最大尺度、最小尺度和尺度数量，计算在最大尺度和最小尺度之间等距分布的每个尺度，并确定检测点的识别灰度阈值和检测点的筛选面积重叠阈值；

采用GPU并行处理的方式，对林区的灰度图像进行多尺度LoG滤波，得到LoG滤波尺度空间图像；

采用GPU并行处理的方式，对LoG滤波尺度空间图像进行局部最大值滤波，得到潜在检测点；

采用GPU并行处理的方式，对潜在检测点进行基于面积的组合筛选：在两两组合的潜在检测点中，若面积重叠度小于所述筛选面积重叠阈值，则删除面积较小的潜在检测点，将剩余的潜在检测点作为初始检测点，

所述步骤对林区的灰度图像进行多尺度LoG滤波包括：采用GPU并行处理的方式，针对每个高斯滤波尺度sigma进行如下滤波，并合并为三维LoG滤波图像数组：

对林区的灰度图像进行列方向LoG滤波，然后进行行方向高斯滤波，得到M_logx；

对林区的灰度图像进行行方向LoG滤波，然后进行列方向高斯滤波，得到M_logy；

M_log＝(M_logx+M_logy)*(-sigma²)，

所述步骤对LoG滤波尺度空间图像进行局部最大值滤波，得到潜在检测点包括：

采用GPU并行处理的方式，以窗口尺寸为3个像素，分别求取三维LoG滤波图像数组在x,y,z方向的最大值，并生成最大值立方体；

对比最大值立方体和三维LoG滤波图像数组，将满足下式的像素作为潜在的检测点：

Image_LoG(x,y,z)＝Image_max(x,y,z)，且Image_LoG(x,y,z)>Threshold

式中，Threshold表示检测点识别灰度阈值；

输出所有潜在检测点的x,y坐标，以及z轴对应的sigma值，

所述步骤对潜在检测点进行基于面积的组合筛选包括：

采用GPU并行处理的方式，根据潜在检测点数量N_BLOBS，利用下式计算两两组合的总数N_combination；

N_combination＝N_BLOBS*(N_BLOBS-1)/2

启动一维数组，为每个块、线程计算其对应的组合一维序号，采用迭代法计算其对应的两个潜在检测点的序号；

利用下式计算重叠度OL：

R₁＝sigma₁*2^1/2,R₂＝sigma₂*2^1/2

Distance＝((x₁–x₂)²+(y₁–y₂)²)^1/2

如果Distance>R₁+R₂

OL＝0

如果Distance≤|R₁-R₂|

OL＝1

否则：

Ratio₁＝(Distance²+R₁ ²–R₂ ²)/(2.0*Distance*R₁)

且Ratio₁＝max(Ratio₁,-1.0)，Ratio₁＝min(Ratio₁,1.0)

Acos₁＝arccos(Ratio₁)

Ratio₂＝(Distance²+R₂ ²–R₁ ²)/(2.0*Distance*R₂)

且Ratio₂＝max(Ratio₂,-1.0)，Ratio₂＝min(Ratio₂,1.0)

Acos₂＝arccos(Ratio₂)

A＝-Distance+R₂+R₁

B＝Distance-R₂+R₁

C＝Distance+R₂-R₁

D＝Distance+R₂+R₁

Area＝R₁ ²*Acos₁+R₂ ²*Acos₂-0.5*(|A*B*C*D|)^1/2

A_min＝min(R₁,R₂)

OL＝Area/(π*A_min ²)

式中，sigma₁和sigma₂分别表示高斯滤波尺度1和高斯滤波尺度2，x₁、y₁表示高斯滤波尺度1下，潜在检测点1的坐标，x₂、y₂表示高斯滤波尺度2下，潜在检测点2的坐标，Distance表示欧氏距离；R₁、R₂、Ratio₁、Ratio₂、Acos₁、Acos₂、A、B、C、D、Area、A_min均表示中间变量，

所述步骤对初始检测点进行筛选，将不含绿色的检测点删除包括：

利用下式将无人机拍摄的林区图像从RGB色彩空间转换为CIE-Lab色彩空间，计算检测点a通道的均值，删除均值高于0的检测点，

X＝R*0.4124+G*0.3576+B*0.1805

Y＝R*0.2126+G*0.7152+B*0.0722

Z＝R*0.0193+G*0.1192+B*0.9505

L＝116*f(Y)–16

a＝500*(f(X)-f(Y))

b＝200*(f(Y)-f(Z))

其中，当x>0.008856时，f(x)＝x^(1/3)，x表示X、Y、Z中的任一中间变量，当x<＝0.008856时，f(x)＝(7.787*x)+(16/116)。

2.根据权利要求1所述的无人机图像的林木植株并行提取方法，其特征在于，

还包括步骤：删除a通道均值高于0的检测点之后，输出每个剩余检测点的坐标和半径。