基于CUDA加速的高光谱图像异常目标探测方法

基于cuda加速的高光谱图像异常目标探测方法
技术领域
1.本发明涉及异常目标检测领域，更具体的说是基于cuda加速的高光谱图像异常目标探测方法。


背景技术：

2.高光谱成像传感器通常在图像的每个空间位置测量数十或数百个窄光谱波段中接收光的能量，由于不同物质通常会表现出差异化的光谱特性，利用高光谱成像技术可对观测场景中的物体进行分类与识别，相关技术已经广泛应用于环境检测、农业生产、军事探测等领域。而当目标光谱特性未知时，一般通过寻找偏离图像典型光谱的像元来做目标探测，这种方法通常被称为异常目标探测。
3.目前针对高光谱图像的异常检测技术大多源于统计学中的异常值检测方法，在目标占整体图像较小比例的条件下，图像中的光谱特征由背景主导，此时异常值检测可以通过预设定的度量形式，计算待检测样本与背景分布间的偏离程度，以检测不太可能是背景的对象。高光谱图像中的每个像元都可以表示为由采样光谱组成的高维矢量，因此基于上述原理的异常检测方法通常需度量每个矢量的偏离值，以判别其异常程度。
4.在当前针对高光谱图像的异常检测方法中，最经典是基于广义似然检验的rx算法，rx算法在将图像建模为高斯分布的基础上通过比较当前位置与背景位置的马氏距离来寻找异常值。后续又有人针对rx算法建模简单，没有考虑谱间相关性等问题进行了改进，衍生出了基于核函数的rx，背景联合稀疏表示的方法等，这些方法计算复杂度都高于基准rx算法。随着成像技术的快速发展，高光谱数据的波段数和空间尺寸在不断增大，在常用的cpu平台上，即使是运算复杂度最低的rx算法也无法满足实时运算的需求。近年来随着半导体工艺的进步，gpu计算能力的快速提升，gpu由于优越的并行计算能力被越来越多地应用于大规模图像矩阵的密集计算。本专利对当前检测算法各步骤进行深入分析，将其中串行执行的计算密集型步骤转换为更高效的矩阵计算模式，并设计多个基于cuda的kernel函数提高算法的并行粒度，合理调配计算资源，发挥cpu和gpu平台各自计算优势，实现对高光谱图像中异常目标的实时检测。


技术实现要素：

5.本发明的目的是解决现有高光谱图像异常目标探测方法在硬件平台实现时存在的运算速度慢，无法满足实时性的问题，提出基于cuda加速的高光谱图像异常目标探测方法。针对高光谱图像处理任务，充分利用cuda模型中cpu+gpu的异构运算优势，将串行处理指令放到cpu主机端运行，将包含大量密集数值操作的运算指令放到具有更多计算单元的gpu设备端运行，从而实现高速运算。
6.为实现上述目的，本发明通过以下技术方案来实现：
7.基于cuda加速的高光谱图像异常目标探测方法，包括以下步骤：
8.步骤一、高光谱图像数据获取与重构；
9.cpu主机端从设备存储器中读取高光谱图像数据dh×w×b，将输入数据dh×w×b中的元素按照波段维度优先的顺序重排列，重构为二维数据矩阵in×b，该指令放到在具有更高内存带宽的cpu主机端执行，其中h为图像高度，w为图像宽度，b为图像波段数，n＝h
×
w为单谱段图像的像元数量。
10.步骤二、将去中心化后的差分矩阵rn×b从cpu主机端内存拷贝到gpu设备端显存中。
11.步骤三、二维数据矩阵去中心化；
12.(3.1)首先计算二维数据矩阵in×b在光谱维度的均值向量计算公式为：
[0013][0014]
其中i
ij
表示二维矩阵在第i个像元位置，第j个波段的数值，设备端计算均值的kernel函数采用规约求和的思想，将单波段中包含的n个数据分配给不同数据块，先将属于不同线程块的数据对应加和到一个线程块中，再对线程块内数据规约求和，并按照公式(1)计算数据在光谱维度的均值向量。
[0015]
(3.2)然后用数据矩阵减去对应波段均值，得到去中心化后的差分矩阵rn×b，该矩阵定义如下：
[0016][0017]
真实高光谱数据的波段数b通常为几百个，一般的gpu设备能够保证有足够线程资源分配给每个波段数据，此处的差分kernel函数将单个线程设置为对应单个数据点，函数按照公式(2)将原数据对应波段中的所有数据值减去对应波段均值。
[0018]
步骤四、计算高光谱图像中各像元位置的异常程度，该步骤通过在gpu设备端设计矩阵乘法kernel函数，广义矩阵求逆kernel函数，矩阵hadamard积kernel函数、以及矩阵行求和kernel函数实现，求解过程包括以下几步：
[0019]
(4.1)协方差矩阵求解，由矩阵乘法kernel函数实现，去中心化后的差分矩阵rn×b的协方差矩阵c的计算方法如下：
[0020][0021]
其中(rn×b)
t
表示矩阵rn×b的转置，该步骤中的矩阵乘法涉及大量相同操作浮点数乘加运算，对于输入大小为n
×
b的矩阵，结果矩阵中单个元素的求解包含n次浮点数相乘和n-1次浮点数相加，实际高光谱数据中的n多达到104个数量集以上，超出gpu设备支持的最大线程数，无法直接将每个浮点数相乘及相加运算对应分配给gpu的一个线程。因此在硬件资源约束下，需要将该kernel函数的输入矩阵rn×b分成若干个小的分块，每个分块作为一个单独的cuda线程块，每个线程块并行地执行计算任务，并将每个小矩阵分配给一个gpu计算单元进行计算，将所有小矩阵的计算结果合并，得到最终的矩阵乘积c。
[0022]
(4.2)逆协方差矩阵求解，由广义矩阵求逆kernel函数实现，用于计算协方差矩阵c的逆矩阵c-1
,矩阵逆的求解通过lu分解实现，lu分解将矩阵求逆运算转换为三角矩阵间的乘法运算，首先将协方差矩阵分解为两个三角矩阵：
[0023][0024]
其中包含l
ij
(1《i,j,《b)的为下三角矩阵，包含u
ij
(1《i,j,《b)的为上三角矩阵，上三角矩阵l和下三角矩阵u的逆可通过迭代公式得到，进而计算出原矩阵c的逆为：
[0025]
c-1
＝(lu)-1
＝u-1
l-1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0026]
lu分解求逆kernel函数由getrfbatched和getribatched两个子函数组成：getrfbatched函数首先将待求逆的矩阵c分块。对于每个分块，使用cuda并行编程技术进行lu分解将协方差矩阵c分解为下三角矩阵l和上三角矩阵u；由于lu分解可能存在数值不稳定的问题，getribatched首先判断lu分解是否成功，再将求解矩阵的逆问题转化为求解一系列线性方程组的问题，并每个矩阵分块使用该函数求解线性方程组，最后通过将每个分块的逆矩阵拼接在一起，得到完整的矩阵c的逆矩阵c-1
。
[0027]
(4.3)异常程度矩阵求解，异常程度矩阵通过矩阵乘法kernel函数、矩阵hadamard积kernel函数、以及矩阵行求和kernel函数得到，基准rx算法中各像元位置异常程度的表征为mahalanobis距离，即：
[0028][0029]
其中为高光谱图像中第i个像元位置的光谱向量，本方法通过公式代换将异常算子的求解从串行计算各像元mahalanobis距离优化为矩阵乘法、矩阵点乘以及矩阵行求和运算，式(6)等价于：
[0030][0031]
其中
·
为矩阵乘法运算，为矩阵hadamard积，即两个矩阵对应位置元素相乘，sum(*,1)表示对*位置处矩阵每行求和。公式(7)中矩阵乘法运算可复用步骤(4.1)中的矩阵乘法kernel函数，矩阵hadamard积运算涉及的浮点数乘法次数等于高光谱图像的像元个数与谱段数之积，即n
×
b，受gpu硬件资源限制无法为每个运算单独赋予一个线程，因此矩阵hadamard积kernel函数先将数据矩阵分块，每个分块作为一个单独的cuda线程块，完成乘法计算后再将结果合并，矩阵行求和运算则是步骤(3.1)中求矩阵均值kernel函数的一个子函数，同样可以复用，经过上述计算最终可得到结果矩阵d
rx
。
[0032]
步骤五、将结果矩阵d
rx
从gpu设备端拷贝到cpu主机端，在cpu主机端申请内存空间，存放结果矩阵数据。
[0033]
步骤六、根据每个待检测样本的异常程度统计在不同阈值条件下的检测率与虚警率，根据可接受虚警α确定对应的全局阈值t，得到最终检测结果如下：
[0034][0035]
其中rd(x)＝1代表该像元为异常像元，rd(x)＝0则表示该像元无异常，此部分涉及大量分支判断逻辑操作，且运算复杂度不高，放到cpu主机端完成，cpu完成异常判断后，将结果保存到存储器设备中。
[0036]
与现有技术相比，本发明具有以下优点：
[0037]
本发明利用gpu并行处理速度快的特点，在gpu设备端设计矩阵均值与差分kernel函数，协方差矩阵求解kernel函数，基于lu分解的广义矩阵逆求解kernel函数，矩阵乘法kernel函数在gpu设备端以多线程，共享内存的形式提高计算速度；
[0038]
上述步骤中数据均以单精度浮点数形式存储，以满足计算精度；
[0039]
通过公式代换进行异常检测算子求解方法优化，将其从串行计算各像元mahalanobis距离优化为矩阵点乘运算，并设计异常检测算子求解kernel函数，相较于现有运行在cpu主机端的基准rx方法，显著提高计算效率，能够满足高光谱图像异常目标检测的运算实时性。
附图说明
[0040]
图1为基准rx算法流程图；
[0041]
图2为本发明方法数据处理及异构编程模型设计图；
[0042]
图3为待检测高光谱图像的真值图；
[0043]
图4a为采用基准rx算法的处理结果图；
[0044]
图4b为采用本发明方法的处理结果图；
具体实施方式
[0045]
基于cuda加速的高光谱图像异常目标探测方法，包括以下步骤：
[0046]
步骤一、高光谱数据获取与重构；
[0047]
cpu主机端从设备存储器中读取高光谱图像数据dh×w×b，将输入数据dh×w×b中的元素按照波段维度优先的顺序重排列，重构为二维数据矩阵in×b，该指令放到在具有更高内存带宽的cpu主机端执行，其中h为图像高度，w为图像宽度，b为图像波段数，n＝h
×
w为单谱段图像的像元数量。
[0048]
步骤二、将去中心化后的差分矩阵rn×b从cpu主机端拷贝到gpu设备端，首先在gpu设备端申请显存空间，之后按显存空间地址拷贝输入的差分矩阵数据。
[0049]
步骤三、二维数据矩阵去中心化；
[0050]
对于二维数据矩阵in×b，计算高光谱数据在各光谱维度的均值计算方法为
[0051][0052]
其中i
ij
表示二维矩阵在第i个像元位置，第j个波段的数值，设备端计算均值的kernel函数采用规约求和的思想，将单波段中包含的n个数据分配给不同数据块，先将属于不同线程块的数据对应加和到一个线程块中，再对线程块内数据规约求和，并按照公式(1)计算数据在光谱维度的均值向量。
[0053]
去中心化第二步是用数据矩阵减去对应波段均值，得到区中心化后的差分矩阵rn×b，该矩阵定义如下：
[0054]
[0055]
真实高光谱数据的波段数b通常为几百个，一般的gpu设备能够保证有足够线程资源分配给每个波段数据，此处的差分kernel函数将单个线程设置为对应单个数据点，函数按照公式(2)将原数据对应波段中的所有数据值减去对应波段均值。
[0056]
步骤四、计算高光谱图像各像元位置异常程度，该步骤通过在gpu设备端设计矩阵乘法kernel函数，广义矩阵求逆kernel函数，矩阵hadamard积kernel函数、以及矩阵行求和kernel函数实现，求解过程包括以下几步：
[0057]
(4.1)协方差矩阵求解，由矩阵乘法kernel函数实现，去中心化后的差分矩阵rn×b的协方差矩阵c的计算方法如下：
[0058][0059]
其中(rn×b)
t
表示矩阵rn×b的转置，该步骤中的矩阵乘法涉及大量相同操作浮点数乘加运算，对于输入大小为n
×
b的矩阵，结果矩阵中单个元素的求解包含n次浮点数相乘和n-1次浮点数相加，实际高光谱数据中的n多达到104个数量集以上，超出gpu设备支持的最大线程数，无法直接将每个浮点数相乘及相加运算对应分配给gpu的一个线程。因此在硬件资源约束下，需要将该kernel函数的输入矩阵rn×b分成若干个小的分块，每个分块作为一个单独的cuda线程块，每个线程块并行地执行计算任务，并将每个小矩阵分配给一个gpu计算单元进行计算，同时将保存在gpu设备端全局内存中的数据以分块的形式拷贝到线程块的共享内存中，以提升数据访存速度，最后将所有小矩阵的计算结果合并得到最终的矩阵乘积c。
[0060]
(4.2)逆协方差矩阵求解，由广义矩阵求逆kernel函数实现，用于计算协方差矩阵c的逆矩阵c-1
,矩阵逆的求解通过lu分解实现，将矩阵求逆运算转换为三角矩阵间的乘法运算，首先将协方差矩阵分解为两个三角矩阵：
[0061][0062]
其中包含l
ij
(1《i,j,《b)的为下三角矩阵，包含u
ij
(1《i,j,《b)的为上三角矩阵，上三角矩阵l和下三角矩阵u的逆可通过迭代公式得到，进而计算出原矩阵c的逆为：
[0063]
c-1
＝(lu)-1
＝u-1
l-1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0064]
lu分解求逆kernel函数由getrfbatched和getribatched两个子函数组成：getrfbatched函数首先将待求逆的矩阵c分块。对于每个分块，使用cuda并行编程技术进行lu分解将协方差矩阵c分解为下三角矩阵l和上三角矩阵u；由于lu分解可能存在数值不稳定的问题，getribatched首先判断lu分解是否成功，再将求解矩阵逆的问题转化为求解一系列线性方程组的问题，并每个矩阵分块使用该函数求解线性方程组，最后通过将每个分块的逆矩阵拼接在一起，得到完整的矩阵c的逆矩阵c-1
。
[0065]
(4.3)异常程度矩阵求解，异常程度矩阵通过矩阵乘法kernel函数、矩阵hadamard积kernel函数、以及矩阵行求和kernel函数计算得到，基准rx算法中各像元位置异常程度的表征为mahalanobis距离，即：
[0066]
[0067]
其中为高光谱图像中第i个像元位置的光谱向量，从式(6)可以看出异常程度矩阵的计算需要串行访问高光谱数据中的每个像元位置的光谱向量，时间成本大，本方法通过公式代换将异常算子的求解从串行计算各像元mahalanobis距离优化为矩阵乘法、矩阵点乘以及矩阵行求和运算，式(6)等价于：
[0068][0069]
其中
·
为矩阵乘法运算，为矩阵hadamard积，即两个矩阵对应位置元素相乘，sum(*,1)表示对*位置处矩阵的每行求和。公式(7)中矩阵乘法运算可复用步骤(4.1)中的矩阵乘法kernel函数，矩阵hadamard积运算涉及的浮点数乘法次数等于高光谱图像的像元个数与谱段数之积，即n
×
b，受gpu硬件资源限制无法为每个运算单独赋予一个线程，因此矩阵hadamard积kernel函数先将数据矩阵分块，每个分块作为一个单独的cuda线程块，完成乘法计算后再将结果合并，矩阵行求和运算则是步骤(3.1)中求矩阵均值kernel函数的一个子函数，同样可以复用，经过上述计算最终可得到结果矩阵d
rx
。
[0070]
步骤五、将结果矩阵d
rx
从gpu设备端拷贝到cpu主机端，首先在cpu主机端申请内存空间，然后按内存地址进行结果矩阵的数据拷贝。
[0071]
步骤六、根据每个待检测样本的异常程度统计在不同阈值条件下的检测率与虚警率，根据可接受虚警α确定对应的全局阈值t，得到最终检测结果如下：
[0072][0073]
其中rd(x)＝1代表该像元为异常像元，rd(x)＝0则表示该像元无异常，之后将结果结果从cpu主机端的内存中拷贝到数据存储器。
[0074]
参照图2，是上述实施方法在处理真实高光谱数据的流程与涉及到的存储器设备，cpu主机设备，gpu设备间的关系。
[0075]
真实高光谱图像处理实验
[0076]
下面分别基于cpu串行实现异常检测算法以及基于cuda并行实现算法开展对比实验，进一步说明本发明的有效性。
[0077]
1、实验条件
[0078]
数据：真实高光谱图像数据集hydice urban dataset，图像高度h＝80，宽度w＝100，波段数b＝162；
[0079]
硬件条件：gpu型号gtx3070，cpu型号intel i7-11800h；
[0080]
2、实验内容
[0081]
基于真实高光谱图像数据集urban dataset，分别实现只运行在cpu主机端的基准rx处理程序和采用本发明方法的异构编程处理程序，并通过分段计时统计两种实现方式的运算耗时，结果如表1所示。
[0082]
表1基于cpu的串行实现方法和基于cuda的并行实现方法用时对比表
[0083]
步骤执行操作名称cpu运算耗时gpu运算耗时1数据预处理0.0240.0242数据去中心化0.0030.0013协方差矩阵计算0.5510.003
4协方差矩阵求逆0.0160.0115mahalanobis距离计算0.6360.017总 1.2300.056
[0084]
其中步骤1的数据预处理包括本发明方法中的高光谱数据获取与重构，在两种处理方法中该步骤均需要在cpu设备端执行。
[0085]
3、结果分析
[0086]
参照图3为标示待检测高光谱数据中异常目标位置的真值图。参照图4a和图4b，分别为采用cpu串行方法和gpu并行方法得到的结果图，二者计算结果一致；
[0087]
表1为基于cpu的串行实现方法和基于cuda的并行实现方法用时对比表，可以看出本发明提出的方法将运行耗时从1.230提减少到0.056，加速比达到22倍，且0.056的耗时可以满足利用高光谱图像进行异常目标实时探测的任务需求。
[0088]
本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

技术特征：
1.基于cuda加速的高光谱图像异常目标探测方法，其特征在于，包括以下步骤：步骤一、cpu主机端从存储器设备中读取高光谱图像数据d
h
×
w
×
b
，并在主机端将其重构为二维数据矩阵i
n
×
b
，其中h为图像高度，w为图像宽度，b为图像波段数，n＝h
×
w为单波段图像像素数；步骤二、将去中心化后的差分矩阵r
n
×
b
从cpu主机端拷贝到gpu设备端；步骤三、二维数据矩阵去中心化，对于二维数据矩阵i
n
×
b
，计算高光谱数据在光谱维度的均值μ
i1
×
b
，并用数据矩阵减去对应波段均值，得到差分矩阵r
b
×
n
；步骤四、高光谱图像各像素位置异常程度求解，根据差分矩阵r
n
×
b
，在gpu端配置多线程kernel函数计算整体协方差矩阵c，并用lu分解方法计算逆协方差矩阵c-1
，改进基准rx算法中通过遍历访问各像元光谱向量计算mahalanobis距离的过程，优化为矩阵乘法和矩阵点乘运算，并在设备端建立对应kernel函数并行处理得到异常检测结果；步骤五、将结果矩阵从gpu设备端拷贝到cpu主机端；步骤六、根据每个待检测样本的异常程度统计在不同阈值条件下的检测率与虚警率，根据可接受虚警α确定对应的全局阈值t，得到最终检测结果如下：2.根据权利要求1所述的基于cuda加速的高光谱图像异常目标探测方法，其特征在于，步骤三中的数据矩阵去中心化在设备端建立矩阵均值与差分kernel函数，根据数据大小与硬件资源约束分配线程，并行计算得出结果。3.根据权利要求1所述的基于cuda加速的高光谱图像异常目标探测方法，其特征在于，步骤四优化计算各像元mahalanobis距离的计算过程，将其从串行访问各像元的光谱向量优化为矩阵乘法运算、矩阵点乘运算和矩阵行求和运算。4.根据权利要求1所述的基于cuda加速的高光谱图像异常目标探测方法，其特征在于，步骤四中为求得到异常检测结果，在设备端建立协方差矩阵求解kernel函数，基于lu分解的广义逆求解kernel函数，矩阵乘法kernel函数、矩阵hadamard积kernel函数以及矩阵行求和kernel函数。5.根据权利要求1所述的基于cuda加速的高光谱图像异常目标探测方法，其特征在于，步骤六中根据差异度量结果分割目标与背景像素的方法为：根据每个待检测样本的异常程度统计在不同阈值条件下的检测率与虚警率，根据可接受虚警α确定对应的全局阈值t，得到最终检测结果如下：其中rd(x)＝1代表该像元为异常像元，rd(x)＝0则表示该像元位置无异常。

技术总结
基于CUDA加速的高光谱图像异常目标探测方法，本发明涉及异常目标检测领域。本发明的目的是解决当前利用高光谱图像进行异常目标探测实时性差的问题。提出了基于CUDA加速的探测方法，该方法主要包括以下过程：1、CPU主机端获取并通过数据重构预处理高光谱图像；2、将预处理数据从CPU主机端的内存拷贝到GPU设备端的显存；3、设计矩阵均值与差分kernel函数进行数据去中心化，设计矩阵乘法kernel函数、广义矩阵求逆kernel函数、Hadamard积kernel函数、以及行求和kernel函数进行协方差矩阵、逆协方差矩阵以及异常检测结果的计算；4、在CPU主机端进行阈值分割并将检测结果保存到存储器端。本发明基于异构编程模型实现，能够满足对高光谱图像数据的实时处理需求。谱图像数据的实时处理需求。谱图像数据的实时处理需求。


技术研发人员：王文正 武润培 唐林波 聂晓风 高鹏程
受保护的技术使用者：北京理工大学前沿技术研究院
技术研发日：2023.03.12
技术公布日：2023/7/12