Hyperspectral Image Anomaly Target Detection Method Based on CUDA Acceleration

Abstract

Hyperspectral Image Anomaly Target Detection Method Accelerated by CUDA. This invention relates to the field of anomaly target detection. The purpose of this invention is to solve the problem of poor real-time performance in detecting anomaly targets using hyperspectral images. A detection method accelerated by CUDA is proposed, which mainly includes the following processes: 1. The CPU host obtains and preprocesses the hyperspectral image through data reconstruction; 2. Copy the preprocessed data from the memory of the CPU host to the video memory of the GPU device; 3. Design matrix mean and difference kernel functions for data de-centralization, and design matrix multiplication kernel functions, generalized matrix inversion kernel functions, Hadamard product kernel functions, and row summation kernel functions for the calculation of covariance matrices, inverse covariance matrices, and anomaly detection results; 4. Perform threshold segmentation on the CPU host and save the detection results to the memory. This invention is implemented based on a heterogeneous programming model and can meet the real-time processing requirements of hyperspectral image data.

Description

Technical Field

[0001] The present invention relates to the field of abnormal target detection, and more specifically to a hyperspectral image abnormal target detection method based on CUDA acceleration. Background Art

[0002] Hyperspectral imaging sensors typically measure the energy of light received in dozens or hundreds of narrow spectral bands at each spatial position of an image. Since different substances usually exhibit different spectral characteristics, hyperspectral imaging technology can be used to classify and identify objects in the observed scene, and related technologies have been widely applied in fields such as environmental detection, agricultural production, and military detection. When the spectral characteristics of the target are unknown, target detection is generally performed by finding pixels that deviate from the typical spectrum of the image, and this method is usually referred to as abnormal target detection.

[0003] Currently, most of the abnormal detection techniques for hyperspectral images are derived from the outlier detection methods in statistics. Under the condition that the target accounts for a small proportion of the overall image, the spectral features in the image are dominated by the background. At this time, outlier detection can calculate the deviation degree between the sample to be detected and the background distribution through a preset metric form to detect objects that are unlikely to be the background. Each pixel in a hyperspectral image can be represented as a high-dimensional vector composed of sampled spectra. Therefore, the abnormal detection methods based on the above principle usually need to measure the deviation value of each vector to determine its degree of abnormality.

[0004] Among the current abnormal detection methods for hyperspectral images, the most classic one is the RX algorithm based on the generalized likelihood ratio test. The RX algorithm finds outliers by comparing the Mahalanobis distance between the current position and the background position on the basis of modeling the image as a Gaussian distribution. Subsequently, some people improved the RX algorithm for problems such as simple modeling and no consideration of spectral correlation, and derived methods such as the RX based on kernel function and the background joint sparse representation method. The computational complexity of these methods is higher than that of the benchmark RX algorithm. With the rapid development of imaging technology, the number of bands and the spatial size of hyperspectral data are constantly increasing. On the commonly used CPU platform, even the RX algorithm with the lowest computational complexity cannot meet the requirements of real-time operation. In recent years, with the progress of semiconductor technology, the computing power of GPUs has been rapidly improved. Due to its superior parallel computing ability, GPUs are increasingly applied to the intensive computing of large-scale image matrices. This patent deeply analyzes each step of the current detection algorithm, converts the computationally intensive steps executed serially into a more efficient matrix calculation mode, designs multiple CUDA-based kernel functions to improve the parallel granularity of the algorithm, rationally allocate computing resources, and give full play to the respective computing advantages of the CPU and GPU platforms to achieve real-time detection of abnormal targets in hyperspectral images. Summary of the Invention

[0005] The object of the present invention is to solve the problem that the existing hyperspectral image anomaly target detection method has a slow operation speed when implemented on a hardware platform and cannot meet the real-time requirement, and to propose a hyperspectral image anomaly target detection method based on CUDA acceleration. For the hyperspectral image processing task, the heterogeneous computing advantage of CPU+GPU in the CUDA model is fully utilized. The serial processing instructions are run on the CPU host side, and the operation instructions containing a large number of intensive numerical operations are run on the GPU device side with more computing units, so as to achieve high-speed operation.

[0006] To achieve the above object, the present invention is realized through the following technical solutions:

[0007] A hyperspectral image anomaly target detection method based on CUDA acceleration includes the following steps:

[0008] Step 1: Hyperspectral image data acquisition and reconstruction;

[0009] The CPU host reads the hyperspectral image data D H×W×B from the device memory, rearranges the elements in the input data D H×W×B in the order of priority in the band dimension, and reconstructs it into a two-dimensional data matrix I N×B . This instruction is executed on the CPU host with higher memory bandwidth, where H is the image height, W is the image width, B is the number of image bands, and N = H×W is the number of pixels in a single spectral band image.

[0010] Step 2: Copy the decentralized difference matrix R N×B from the CPU host memory to the GPU device side video memory.

[0011] Step 3: Decentralize the two-dimensional data matrix;

[0012] (3.1) First, calculate the mean vector N×B of the two-dimensional data matrix I in the spectral dimension. The calculation formula is:

[0013]

[0014] where I ij represents the value of the two-dimensional matrix at the i-th pixel position and the j-th band. The kernel function for calculating the mean on the device side adopts the idea of reduction summation. The N data included in a single band are assigned to different data blocks. First, the data belonging to different thread blocks are added up to a thread block, and then the data within the thread block is reduced and summed, and the mean vector of the data in the spectral dimension is calculated according to formula (1).

[0015] (3.2) Then subtract the corresponding band mean from the data matrix to obtain the decentralized difference matrix R N×B, the matrix is defined as follows:

[0016]

[0017] The number of bands B of real hyperspectral data is usually several hundred. General GPU devices can ensure that there are sufficient thread resources to be allocated to each band of data. Here, the differential kernel function sets a single thread to correspond to a single data point, and the function subtracts all the data values in the corresponding band of the original data from the mean value of the corresponding band according to formula (2).

[0018] Step 4: Calculate the anomaly degree of each pixel position in the hyperspectral image. This step is implemented by designing matrix multiplication kernel functions, generalized matrix inversion kernel functions, matrix Hadamard product kernel functions, and matrix row summation kernel functions on the GPU device side. The solution process includes the following steps:

[0019] (4.1) Solve the covariance matrix, which is implemented by the matrix multiplication kernel function. The calculation method of the covariance matrix C of the decentralized difference matrix R N×B is as follows:

[0020]

[0021] where (R N×B ) T represents the transpose of matrix R N×B . The matrix multiplication in this step involves a large number of identical floating-point multiply-add operations. For a matrix with an input size of N×B, the solution of a single element in the result matrix contains N floating-point multiplications and N - 1 floating-point additions. In actual hyperspectral data, N can reach more than 10 4 quantity sets, exceeding the maximum number of threads supported by the GPU device. It is impossible to directly allocate each floating-point multiplication and addition operation to a single thread of the GPU. Therefore, under the constraint of hardware resources, the input matrix R N×B of this kernel function needs to be divided into several small blocks, each block is used as a separate CUDA thread block, each thread block executes the calculation task in parallel, and each small matrix is allocated to a GPU computing unit for calculation. The calculation results of all small matrices are combined to obtain the final matrix product C.

[0022] (4.2) Solve the inverse covariance matrix, which is implemented by the generalized matrix inversion kernel function and is used to calculate the inverse matrix C -1 of the covariance matrix C. The solution of the matrix inverse is realized through LU decomposition. LU decomposition converts the matrix inverse operation into a multiplication operation between triangular matrices. First, the covariance matrix is decomposed into two triangular matrices:

[0023]

[0024] It contains L ij For \(1 \lt i,j \lt B\), it is a lower triangular matrix, containing U ij For \(1 \lt i,j \lt B\), it is an upper triangular matrix. The inverses of the upper triangular matrix L and the lower triangular matrix U can be obtained through iterative formulas, and then the inverse of the original matrix C is calculated as:

[0025] C -1 =(LU) -1 =U -1 L -1 (5)

[0026] The LU decomposition inverse kernel function consists of two sub-functions, getrfBatched and getriBatched: The getrfBatched function first partitions the matrix C to be inverted. For each partition, the CUDA parallel programming technique is used to perform LU decomposition to decompose the covariance matrix C into a lower triangular matrix L and an upper triangular matrix U; Since LU decomposition may have numerical instability problems, getriBatched first determines whether the LU decomposition is successful, then transforms the problem of solving the inverse of the matrix into the problem of solving a series of linear equations, and uses this function to solve the linear equations for each matrix partition. Finally, by piecing together the inverse matrices of each partition, the inverse matrix C of the complete matrix C is obtained -1 .

[0027] (4.3) Solving the anomaly degree matrix. The anomaly degree matrix is obtained through the matrix multiplication kernel function, the matrix Hadamard product kernel function, and the matrix row sum kernel function. In the benchmark RX algorithm, the anomaly degree of each pixel position is characterized by the Mahalanobis distance, that is:

[0028]

[0029] where is the spectral vector at the \(i\)-th pixel position in the hyperspectral image. This method optimizes the solution of the anomaly operator from serial calculation of the Mahalanobis distance of each pixel to matrix multiplication, matrix dot product, and matrix row sum operations. Equation (6) is equivalent to:

[0030]

[0031] where · represents matrix multiplication operation, is the matrix Hadamard product, that is, the elements at the corresponding positions of two matrices are multiplied, and sum(*, 1) represents the sum of each row of the matrix at the * position. The matrix multiplication operation in formula (7) can reuse the matrix multiplication kernel function in step (4.1). The number of floating-point multiplications involved in the matrix Hadamard product operation is equal to the product of the number of pixels and the number of spectral bands of the hyperspectral image, that is, N×B. Due to the limitation of GPU hardware resources, it is impossible to assign a separate thread to each operation. Therefore, the matrix Hadamard product kernel function first divides the data matrix into blocks, each block is used as a separate CUDA thread block, and after completing the multiplication calculation, the results are merged. The matrix row summation operation is a sub-function of the matrix mean kernel function in step (3.1) and can also be reused. After the above calculations, the result matrix D can be finally obtained RX 。

[0032] Step Five: Copy the result matrix D RX from the GPU device side to the CPU host side, and apply for memory space on the CPU host side to store the result matrix data.

[0033] Step Six: According to the anomaly degree of each sample to be detected, count the detection rate and false alarm rate under different threshold conditions, and determine the corresponding global threshold T according to the acceptable false alarm α. The final detection results are as follows:

[0034]

[0035] Among them, RD(x) = 1 indicates that the pixel is an abnormal pixel, and RD(x) = 0 indicates that the pixel is normal. This part involves a large number of branch judgment logic operations, and the operation complexity is not high, so it is completed on the CPU host side. After the CPU completes the anomaly judgment, the results are saved to the memory device.

[0036] Compared with the prior art, the present invention has the following advantages:

[0037] The present invention utilizes the characteristics of fast parallel processing speed of the GPU, designs matrix mean and difference kernel functions, covariance matrix solving kernel functions, and generalized matrix inverse solving kernel functions based on LU decomposition on the GPU device side. The matrix multiplication kernel function improves the calculation speed in the form of multi-threading and shared memory on the GPU device side;

[0038] In the above steps, the data is stored in the form of single-precision floating-point numbers to meet the calculation accuracy;

[0039] Optimize the anomaly detection operator solution method through formula substitution, optimize it from serial calculation of Mahalanobis distance for each pixel to matrix dot multiplication operation, and design the kernel function for solving the anomaly detection operator. Compared with the existing benchmark RX method running on the CPU host side, the calculation efficiency is significantly improved, and the real-time operation of hyperspectral image anomaly target detection can be satisfied. Description of the Drawings

[0040] Figure 1 It is the flowchart of the benchmark RX algorithm;

[0041] Figure 2 It is the data processing and heterogeneous programming model design diagram of the method of the present invention;

[0042] Figure 3 It is the ground truth map of the hyperspectral image to be detected;

[0043] Figure 4a It is the processing result map using the benchmark RX algorithm;

[0044] Figure 4b It is the processing result map using the method of the present invention; Detailed Embodiment

[0045] The hyperspectral image anomaly target detection method based on CUDA acceleration includes the following steps:

[0046] Step 1: Hyperspectral data acquisition and reconstruction;

[0047] The CPU host reads the hyperspectral image data D from the device memory H×W×B , and rearranges the elements in the input data D H×W×B in the order of giving priority to the band dimension to reconstruct it into a two-dimensional data matrix I N×B . This instruction is executed on the CPU host with higher memory bandwidth, where H is the image height, W is the image width, B is the number of image bands, and N = H×W is the number of pixels in a single spectral band image.

[0048] Step 2: Copy the decentralized difference matrix R N×B from the CPU host to the GPU device. First, apply for video memory space on the GPU device, and then copy the input difference matrix data according to the video memory space address.

[0049] Step 3: Decentralize the two-dimensional data matrix;

[0050] For the two-dimensional data matrix I N×B , calculate the mean value of the hyperspectral data in each spectral dimension The calculation method is

[0051]

[0052] where I ij represents the value of the two-dimensional matrix at the i-th pixel position and the j-th band. The kernel function for calculating the mean at the device end adopts the idea of reduction summation, distributing the N data contained in a single band to different data blocks. First, the data belonging to different thread blocks are correspondingly summed up to a single thread block, then the data within the thread block is reduced and summed up, and the mean vector of the data in the spectral dimension is calculated according to formula (1).

[0053] The second step of de-centralization is to subtract the corresponding band mean from the data matrix to obtain the de-centralized difference matrix R N ×B , and this matrix is defined as follows:

[0054]

[0055] The number of bands B of real hyperspectral data is usually several hundred. General GPU devices can ensure that there are sufficient thread resources to be allocated to each band of data. Here, the difference kernel function sets a single thread to correspond to a single data point, and the function subtracts all the data values in the corresponding band of the original data from the corresponding band mean according to formula (2).

[0056] Step 4: Calculate the anomaly degree of each pixel position in the hyperspectral image. This step is implemented by designing a matrix multiplication kernel function, a generalized matrix inversion kernel function, a matrix Hadamard product kernel function, and a matrix row summation kernel function at the GPU device end. The solution process includes the following steps:

[0057] (4.1) Solve the covariance matrix, which is implemented by the matrix multiplication kernel function. The calculation method of the covariance matrix C of the de-centralized difference matrix R N×B is as follows:

[0058]

[0059] where (R N×B ) T represents the transpose of matrix R N×B . The matrix multiplication in this step involves a large number of identical floating-point multiply-add operations. For a matrix with an input size of N×B, the solution of a single element in the result matrix contains N floating-point multiplications and N - 1 floating-point additions. In actual hyperspectral data, N reaches more than 10 4 quantity sets, exceeding the maximum number of threads supported by the GPU device. It is impossible to directly allocate each floating-point multiplication and addition operation to a single thread of the GPU. Therefore, under the constraint of hardware resources, the input matrix R of this kernel function needs to be N×BIt is divided into several small blocks, and each block serves as an independent CUDA thread block. Each thread block executes the computing tasks in parallel, and each small matrix is assigned to a GPU computing unit for calculation. Meanwhile, the data stored in the global memory on the GPU device is copied into the shared memory of the thread block in a block form to improve the data access speed. Finally, the calculation results of all small matrices are merged to obtain the final matrix product C.

[0060] (4.2) Inverse covariance matrix solution is implemented by the generalized matrix inversion kernel function, which is used to calculate the inverse matrix C of the covariance matrix C. -1 , The solution of the matrix inverse is realized by LU decomposition, which converts the matrix inversion operation into a multiplication operation between triangular matrices. First, the covariance matrix is decomposed into two triangular matrices:

[0061]

[0062] which includes L ij (1 < i, j, < B) is the lower triangular matrix, which includes U ij (1 < i, j, < B) is the upper triangular matrix. The inverses of the upper triangular matrix L and the lower triangular matrix U can be obtained through iterative formulas, and then the inverse of the original matrix C is calculated as:

[0063] C -1 =(LU) -1 =U -1 L -1 (5)

[0064] The LU decomposition inverse kernel function consists of two sub-functions, getrfBatched and getriBatched: The getrfBatched function first divides the matrix C to be inverted into blocks. For each block, the CUDA parallel programming technology is used to perform LU decomposition to decompose the covariance matrix C into the lower triangular matrix L and the upper triangular matrix U; Since LU decomposition may have numerical instability problems, getriBatched first judges whether the LU decomposition is successful, then converts the problem of solving the matrix inverse into the problem of solving a series of linear equations, and uses this function to solve the linear equations for each matrix block. Finally, by splicing the inverse matrices of each block together, the inverse matrix C of the complete matrix C is obtained. -1 .

[0065] (4.3) Abnormality degree matrix solution: The abnormality degree matrix is calculated through the matrix multiplication kernel function, the matrix Hadamard product kernel function, and the matrix row sum kernel function. The characterization of the abnormality degree at each pixel position in the benchmark RX algorithm is the Mahalanobis distance, that is:

[0066]

[0067] Among them is the spectral vector at the \(i\)-th pixel position in the hyperspectral image. It can be seen from Equation (6) that calculating the anomaly degree matrix requires serial access to the spectral vectors at each pixel position in the hyperspectral data, which has a large time cost. This method optimizes the solution of the anomaly operator from serial calculation of the Mahalanobis distance of each pixel to matrix multiplication, matrix dot product, and matrix row summation operations. Equation (6) is equivalent to:

[0068]

[0069] where · represents matrix multiplication operation is the matrix Hadamard product, that is, the elements at the corresponding positions of two matrices are multiplied. sum(*,1) represents summing each row of the matrix at the * position. The matrix multiplication operation in Formula (7) can reuse the matrix multiplication kernel function in step (4.1). The number of floating-point multiplications involved in the matrix Hadamard product operation is equal to the product of the number of pixels and the number of spectral bands in the hyperspectral image, that is, N×B. Due to GPU hardware resource limitations, it is impossible to assign a separate thread to each operation. Therefore, the matrix Hadamard product kernel function first divides the data matrix into blocks, each block is used as a separate CUDA thread block, and the results are merged after the multiplication calculation. The matrix row summation operation is a sub-function of the matrix mean calculation kernel function in step (3.1) and can also be reused. After the above calculations, the result matrix D RX can be obtained finally.

[0070] Step Five: Copy the result matrix D RX from the GPU device side to the CPU host side. First, allocate memory space on the CPU host side, and then copy the data of the result matrix according to the memory address.

[0071] Step Six: According to the anomaly degree of each sample to be detected, count the detection rate and false alarm rate under different threshold conditions, and determine the corresponding global threshold T according to the acceptable false alarm α. The final detection result is as follows:

[0072]

[0073] where RD(x) = 1 represents that the pixel is an abnormal pixel, and RD(x) = 0 indicates that the pixel is normal. Then copy the result from the memory of the CPU host side to the data storage.

[0074] Refer to Figure 2 , which is the relationship between the above implementation method in processing real hyperspectral data and the memory device, CPU host device, and GPU device involved.

[0075] Real Hyperspectral Image Processing Experiment

[0076] The following is a comparative experiment based on the CPU serial implementation of the anomaly detection algorithm and the CUDA parallel implementation algorithm to further illustrate the effectiveness of the present invention.

[0077] 1. Experimental Conditions

[0078] Data: Real hyperspectral image dataset HYDICE urban dataset, image height H = 80, width W = 100, number of bands B = 162;

[0079] Hardware conditions: GPU model GTX3070, CPU model Intel i7-11800H;

[0080] 2. Experimental Content

[0081] Based on the real hyperspectral image dataset urban dataset, a benchmark RX processing program that only runs on the CPU host side and a heterogeneous programming processing program using the method of the present invention are respectively implemented, and the operation time of the two implementation methods is statistically counted by segmented timing. The results are shown in Table 1.

[0082] Table 1 Comparison of the time used by the CPU-based serial implementation method and the CUDA-based parallel implementation method

[0083] Step Execute operation name CPU operation time consumption GPU operation time consumption 1 Data preprocessing 0.024 0.024 2 Data decentralization 0.003 0.001 3 Covariance matrix calculation 0.551 0.003 4 Covariance matrix inversion 0.016 0.011 5 Mahalanobis distance calculation 0.636 0.017 Total 1.230 0.056

[0084] Among them, the data preprocessing in step 1 includes the hyperspectral data acquisition and reconstruction in the method of the present invention, and this step needs to be executed on the CPU device side in both processing methods.

[0085] 3. Result Analysis

[0086] Refer to Figure 3 as the ground truth map for marking the positions of abnormal targets in the hyperspectral data to be detected. Refer to Figure 4a and Figure 4b , which are the result maps obtained by the CPU serial method and the GPU parallel method respectively, and the calculation results of the two are consistent;

[0087] Table 1 is a comparison of the time used by the CPU-based serial implementation method and the CUDA-based parallel implementation method. It can be seen that the method proposed by the present invention reduces the running time from 1.230 to 0.056, and the speedup ratio reaches 22 times. Moreover, the time of 0.056 can meet the task requirements of real-time detection of abnormal targets using hyperspectral images.

[0088] The present invention may also have many other embodiments. Without departing from the spirit and essence of the present invention, those skilled in the art can make various corresponding changes and modifications according to the present invention. However, these corresponding changes and modifications should fall within the protection scope of the appended claims of the present invention.

Claims

1. A hyperspectral image anomaly target detection method based on CUDA acceleration, characterized in that, Including the following steps: Step 1: The CPU host reads hyperspectral image data from the memory device and reconstructs it into a two-dimensional data matrix on the host side , where H is the image height, W is the image width, B is the number of image bands, and N = H × W is the number of pixels in a single-band image; Step 2: Copy the de-centralized difference matrix from the CPU host to the GPU device; Step 3. Decentralize the two-dimensional data matrix. For the two-dimensional data matrix , calculate the mean value of the hyperspectral data in the spectral dimension , and subtract the mean value of the corresponding band from the data matrix to obtain a difference matrix ; Step 4. Solve the anomaly degree of each pixel position in the hyperspectral image. According to the difference matrix , configure a multi-threaded kernel function on the GPU side to calculate the overall covariance matrix C, and use the LU decomposition method to calculate the inverse covariance matrix . Improve the process of calculating the Mahalanobis distance by traversing and accessing the spectral vectors of each pixel in the benchmark RX algorithm, and optimize it into matrix multiplication and matrix dot product operations. Then establish a corresponding kernel function on the device side for parallel processing to obtain the anomaly detection result matrix ; Step 5. Copy the anomaly detection result matrix from the GPU device side to the CPU host side; Step 6: According to the abnormality degree of each sample to be detected, the detection rate and false alarm rate under different threshold conditions are statistically analyzed. According to the acceptable false alarm the corresponding global threshold T is determined, and the final detection result is as follows: ； Among them indicates that this pixel is an abnormal pixel, while it means there is no abnormality at the position of this pixel.

2. The hyperspectral image anomaly target detection method based on CUDA acceleration according to claim 1, wherein, In step three, for data matrix de-centralization, a matrix mean value and a differential kernel function are established on the device side. Threads are allocated according to data size and hardware resource constraints, and the result is obtained through parallel computing.

3. The hyperspectral image abnormal target detection method based on CUDA acceleration according to claim 1, wherein, In step four, the calculation process of the Mahalanobis distance of each pixel is optimized, and it is optimized from serially accessing the spectral vectors of each pixel to matrix multiplication operation, matrix dot multiplication operation, and matrix row summation operation.

4. The hyperspectral image anomaly target detection method based on CUDA acceleration according to claim 1, characterized in that In step four, in order to obtain the anomaly detection result, a covariance matrix solution kernel function is established on the device side, a generalized inverse solution kernel function based on LU decomposition, a matrix multiplication kernel function, a matrix Hadamard product kernel function, and a matrix row summation kernel function.

Images