An Accelerated Computational Method for PCA Dimensionality Reduction Algorithm for High-Resolution Multispectral Data Based on Heterogeneous Computation
By optimizing the PCA dimensionality reduction algorithm through a heterogeneous computing framework and leveraging the computing advantages of GPUs and CPUs, the problem of high CPU resource consumption in high-resolution remote sensing data processing is solved, and efficient data dimensionality reduction processing is achieved.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2022-01-07
- Publication Date
- 2026-03-13
AI Technical Summary
Existing technologies consume high CPU computing resources in the analysis and processing of high-resolution remote sensing data covering multiple time series and wide geographical areas, resulting in long PCA dimensionality reduction calculation times, which affects work efficiency and cost.
We adopt a heterogeneous computing framework based on GPGPU+CPU, optimize the traditional PCA dimensionality reduction algorithm, and utilize the computing resource advantages of GPU and CPU through the cyclic split window algorithm and the Chivens rotation matrix iteration method to achieve parallel and serial computing of high-resolution multispectral data.
It improves the efficiency of PCA dimensionality reduction processing of high-resolution multispectral remote sensing data, reduces computation time, and makes full use of computer hardware resources.
Smart Images
Figure CN114372524B_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the technical field of PCA dimensionality reduction algorithm for high-resolution multispectral data, specifically involving an accelerated computation method for PCA dimensionality reduction algorithm for high-resolution multispectral data based on heterogeneous computing. Background Technology
[0002] Currently, in the field of multispectral and hyperspectral remote sensing data analysis and mining, PCA principal component analysis (PCA) dimensionality reduction of multidimensional or high-dimensional data has become a common method in qualitative analysis and quantitative inversion of remote sensing ground features. With the development of satellite sensor technology and the continuous improvement of remote sensing accuracy requirements, the use of high-resolution and even ultra-high-resolution data is becoming increasingly indispensable. However, the overall scale and volume of data are also increasing, leading to higher data processing time and computational resource consumption, especially CPU resources. Particularly in the analysis and processing of high-resolution remote sensing data covering multiple time series and wide geographical areas, computers with conventional configurations take relatively long times in the PCA dimensionality reduction calculation stage, which becomes a constraint on actual work efficiency and configuration costs. Therefore, optimizing the implementation of traditional CPU-based PCA dimensionality reduction algorithms, while considering computer hardware configurations, to maximize the utilization of traditional computer computing hardware resources to improve the efficiency of PCA principal component analysis for massive datasets is an urgent problem to be solved. Summary of the Invention
[0003] To address the aforementioned technical problems, this invention provides an accelerated computation method for PCA dimensionality reduction algorithms for high-resolution multispectral data based on heterogeneous computing. Addressing the time-consuming and resource-intensive nature of traditional PCA algorithms and computer computing environments when processing massive amounts of remote sensing multispectral data, this invention fully explores and utilizes computer hardware resources. It optimizes the implementation steps of traditional PCA dimensionality reduction algorithms using a heterogeneous computing framework based on GPGPU+CPU, forming a rapid computational processing method for large volumes of high-resolution remote sensing data, thereby improving the overall efficiency of data dimensionality reduction processing.
[0004] To achieve the aforementioned objectives, an accelerated computation method for high-resolution multispectral data PCA dimensionality reduction algorithm based on heterogeneous computing is provided, specifically including the following steps:
[0005] Step 1: Read and load multidimensional data volumes of high-resolution multispectral remote sensing data;
[0006] Step 2: Solve the self-weighted average of non-background pixels in the single-dimensional data using a heterogeneous computation method based on the cyclic split window algorithm;
[0007] Step 3: Achieve centralized standardization of non-background image metadata for single-dimensional data;
[0008] Step 4: Solve for the covariance matrix between each dimension of the data volume using heterogeneous computing methods;
[0009] Step 5: Use the Chivens rotation matrix iteration method to solve for the eigenvalues and corresponding eigenvectors of the covariance matrix;
[0010] Step 6: Solve for the projection of the multidimensional cell vector onto each principal component;
[0011] Step 7: Acquire and store the result data of several principal component bands after projection calculation.
[0012] Compared with the prior art, the beneficial effects of the present invention are at least as follows:
[0013] This invention provides an accelerated computation method for PCA dimensionality reduction algorithm for high-resolution multispectral data based on heterogeneous computing. It solves the problem that existing technologies, when analyzing and processing high-resolution remote sensing data with multi-temporal and wide-area coverage, require relatively long computation times for PCA dimensionality reduction calculations on conventionally configured computers. This method fully explores and utilizes the computing resources of the computer's CPU and GPGPU, achieving accelerated optimization of PCA dimensionality reduction algorithm computation for high-resolution multispectral remote sensing data under a heterogeneous computing framework, especially for large data volumes, thereby improving the efficiency of data dimensionality reduction processing in practical work. Attached Figure Description
[0014] Figure 1 This is a flowchart illustrating the steps of an accelerated computation method for a high-resolution multispectral data PCA dimensionality reduction algorithm based on heterogeneous computing according to the present invention.
[0015] Figure 2 This is a flowchart of the steps of the cyclic split window algorithm of the present invention;
[0016] Figure 3 This is a flowchart illustrating the steps of the present invention for calculating the self-weighted average of non-background pixels in one-dimensional data;
[0017] Figure 4 This is a flowchart illustrating the steps of solving the covariance matrix between different dimensions of a data volume according to the present invention. Detailed Implementation
[0018] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.
[0019] It is understood that the terms "first," "second," etc., used in this application may be used herein to describe various elements, but unless otherwise specified, these elements are not limited by these terms. These terms are used only to distinguish one element from another. For example, without departing from the scope of this application, a first script may be referred to as a second script, and similarly, a second script may be referred to as a first script.
[0020] For reference Figure 1 As shown, this invention provides an accelerated computation method for high-resolution multispectral data PCA dimensionality reduction algorithm based on heterogeneous computing, which is implemented through the following steps:
[0021] Step 1: Read and load multidimensional data volumes of high-resolution multispectral remote sensing data;
[0022] Step 2: Solve the self-weighted average of non-background pixels in the single-dimensional data using a heterogeneous computation method based on the cyclic split window algorithm;
[0023] Step 3: Achieve centralized standardization of non-background image metadata for single-dimensional data;
[0024] Step 4: Solve for the covariance matrix between each dimension of the data volume using heterogeneous computing methods;
[0025] Step 5: Use the Chivens rotation matrix iteration method to solve for the eigenvalues and corresponding eigenvectors of the covariance matrix;
[0026] Step 6: Solve for the projection of the multidimensional cell vector onto each principal component;
[0027] Step 7: Acquire and store the result data of several principal component bands after projection calculation.
[0028] Furthermore, in step one above, the multidimensional data volume of high-resolution multispectral remote sensing data can be understood as a data cube with X, Y, and Z directions. XSize and YSize represent the scale of the multidimensional data volume in the X and Y directions, respectively, and Dimensions represents the dimension of the multidimensional data volume in the Z direction. Each dimension represents a band. In this embodiment, ESA Sentinel-2A multispectral data is used as the experimental data, and the area within the administrative boundary of a certain region is used as the research area. The four commonly used bands of blue, green, red, and near-infrared in the multispectral data are extracted as the processing bands. It should be noted that before using the accelerated calculation method of the high-resolution multispectral data PCA dimensionality reduction algorithm based on heterogeneous computing provided by this invention, it is necessary to preprocess the high-resolution multispectral remote sensing data, including data calibration, correction, band extraction, resampling, and obtaining binarized band data.
[0029] Furthermore, the aforementioned heterogeneous computing method delegates the logic control and data access processes, as well as some computationally insignificant operations, of the PCA dimensionality reduction algorithm to the CPU, while delegating the pixel-level processing of massive amounts of single graphic data to the GPU.
[0030] Specifically, heterogeneous computing is adopted because the traditional CPU chip architecture is primarily designed to handle instruction control, logic judgment, data caching, and computation. This architecture results in a relatively low proportion of computational units responsible for computation within the CPU chip. Furthermore, in terms of computation, the CPU, as a logic judgment and instruction control core, constantly accesses and retrieves data from memory and logic instructions when processing data. Therefore, its computation is better suited for serial data computation. This design architecture and computational method are often inefficient for tasks with massive computational loads, such as massive pixel-level computation or graphics pipeline rendering computation. In contrast, the design of GPU chips focuses primarily on computation. It contains a large number of computational units, but is weaker than the CPU in terms of control and instruction. Therefore, it is better suited for parallel data computation. Its initial design was mainly to handle massive graphics pixel rendering pipeline computation, and it has a relatively obvious advantage in pixel-level computation.
[0031] Further, see references such as Figure 2 As shown, the cyclic split window algorithm in step two above specifically includes the following steps:
[0032] The first step is to obtain the single-dimensional data of a certain band in the multidimensional data volume of high-resolution multispectral remote sensing data, and to represent the image data scale of the single-dimensional data in the X direction as XSize, and the image data scale of the single-dimensional data in the Y direction as YSize.
[0033] The second step is to set the size of the image metadata window within the scope of the single-dimensional data, and to represent the image metadata scale of the image metadata window in the X direction as Winx, and the image metadata scale of the image metadata window in the Y direction as Winy.
[0034] The third step is to scale down the single-dimensional data using a splitting formula. For each pixel in the scaled-down single-dimensional data, the value is the sum of its corresponding non-background pixels in the original single-dimensional data. The splitting formula is described as follows:
[0035]
[0036] Where XSize' and YSize' are the scales of the scaled single-dimensional data in the X and Y directions, respectively; XSize and YSize are the scales of the single-dimensional data in the X and Y directions before scaling down; and Winx and Winy are the scales of the image metadata window in the single-dimensional data in the X and Y directions, respectively.
[0037] Step 4: If, in Step 3 above, the scale of the single-dimensional data in the X or Y direction is insufficient to divisible the scale of the image data window in the X or Y direction, supplement the single-dimensional data with corresponding 0-pixel columns and 0-pixel rows. The number of columns (col) and rows (row) to be supplemented are calculated using the following formula:
[0038]
[0039] Where Winx and Winy are the scales of the image metadata window in the X and Y directions, respectively, and XSize and YSize are the scales of the single-dimensional data in the X and Y directions, respectively.
[0040] Step 5: Repeat the scaling down process on the single-dimensional data until the preset number of processing times is reached, or the scale of the single-dimensional data in both the X and Y directions has been reduced to below the preset scale, then stop the scaling down process.
[0041] Specifically, step two above uses GPGPU computing resources to sum the non-background data of a certain band in high-resolution multispectral remote sensing data. The inventors considered that GPU parallel computing is mainly for computing large sets of relatively independent data, that is, there is no data dependency between each data computing process. However, for data of a given band, the summation of its non-background data is a process of judging and accumulating all pixel data in sequence. Obviously, each summation is based on the previous summation and the accumulation of the next non-background pixel. This computing process is a sequential computing method, which does not satisfy the independence of data computing units. Therefore, by reducing the size of the single-dimensional data through the cyclic split window algorithm in step two, and then putting the single-dimensional data at this scale into the CPU for serial accumulation calculation, there is no increase in computing burden and time consumption. In particular, the process of calculating the single pixel value of the reduced single-dimensional data in the third step of the cyclic split window algorithm can be put into the GPGPU as a parallel computing unit for processing.
[0042] Further, see references such as Figure 3 As shown, step two, which calculates the self-weighted average of non-background pixels in the single-dimensional data, includes the following steps:
[0043] The first step is to obtain B using the cyclic split window algorithm.k Sum(DataB) of all pixels in a single-dimensional data band k );
[0044] Step 2, B k The value of each cell in the single-dimensional data of the band, DataB k Substituting (i, j) into the factor calculation formula, we obtain the intermediate calculation factor W. t (i, j), the factor calculation formula is described as follows:
[0045]
[0046] The third step is to obtain B using the formula for calculating the self-proportional weighted average. k The self-weighted average of single-dimensional data in the band, Avg(DataB) k The formula for calculating the self-proportionally weighted average is described below:
[0047]
[0048] Where xsize and ysize are B respectively k The scale of single-dimensional data in the X and Y directions of the band.
[0049] Specifically, the purpose of calculating the self-weighted average of non-background pixels in single-dimensional data is to avoid confusion and interference between zero-pixel outliers that appear after removing certain noise data such as upper-level cloud data or other noise data from high-resolution multispectral remote sensing data and background pixels during the calculation process. The execution process of the first and third steps mentioned above can be accelerated using the idea of the cyclic split window algorithm.
[0050] Furthermore, step three specifically achieves the central standardization of non-background image data of single-dimensional data through the following calculation formula:
[0051]
[0052] Among them, CenDataB k (i,j) is B k The pixel values of the single-dimensional data of the band after center normalization, DataB k (i,j) is B k The pixel values of a single-dimensional data band, Avg(DataB) k ) is B k The self-weighted average of single-dimensional data for a band.
[0053] Specifically, the calculation of a single cell CenDataB mentioned above can be performed. kThe process of (i,j) is used as a parallel computing factor and fed into the GPGPU for parallel processing, thereby accelerating the central standardization of non-background image data of single-dimensional data.
[0054] Further, see references such as Figure 4 As shown, step four, solving for the covariance matrix between each dimension of the data volume, includes the following steps:
[0055] Step 1: Calculate any two centered bands b of the data volume. m ,b n The intermediate factor TmpDataB is obtained by sequentially substituting all pixels of the single-dimensional data into the intermediate factor calculation formula. mn (i,j), the intermediate factor calculation formula is described as TmpDataB mn (i,j)=CenDataB m (i, j)×CenDataB n (i,j), where CenDataB m (i, j), CenDataB n (i,j) represent band b respectively. m ,b n The pixel values of single-dimensional data after center standardization;
[0056] The second step is to calculate the sum M of all the intermediate calculation factors, specifically by expressing the calculation process of M as follows: Where xsize and ysize are respectively b m ,b n The scale of single-dimensional data in the X and Y directions of the band;
[0057] Step 3: Calculate band b m ,b n Each pixel in the single-dimensional binary data DataB bin The sum of (i, j) and N, specifically the calculation process of N is expressed as follows: Where xsize and ysize are respectively b m ,b n The scale of the single-dimensional binary data of the band in the X and Y directions;
[0058] Step 4: Through calculation Obtain band b m ,b n The covariance of single-dimensional data;
[0059] Step 5: Repeat the above steps until all elements of the covariance matrix have been calculated.
[0060] Specifically, the calculation processes in the first, second, and third steps mentioned above, as well as the process of calculating each element in the covariance matrix, can all be accelerated by utilizing the cyclic split window algorithm, thereby improving the speed of solving the covariance matrix between different dimensions of the data volume.
[0061] Furthermore, in step five above, the eigenvalues and corresponding eigenvectors of the covariance matrix are specifically solved using the Chivens rotation matrix iteration method. This is the conventional calculation process for solving the eigenvalues and corresponding eigenvectors of a matrix in linear algebra, so this calculation process will not be described in detail in this embodiment.
[0062] Furthermore, in step six above, the projection of the multidimensional pixel vectors of the multidimensional data volume onto each principal component can be calculated as a parallel computation factor and placed into the GPGPU for parallel processing to accelerate the execution speed of step six. Specifically, in step six, the orthogonal basis vectors obtained in step five are sorted from largest to smallest according to their corresponding eigenvalues, and the projection calculation process is performed sequentially to obtain a new multidimensional data volume. Each band of this multidimensional data volume is the PCA result of the principal components from largest to smallest. Generally speaking, the first 1-3 bands are the principal component bands that retain more information, while the subsequent bands generally represent information such as noise that can be ignored. Thus, the goal of dimensionality reduction of the multidimensional data volume of high-resolution multispectral remote sensing data is ultimately achieved within the framework of heterogeneous computing optimization.
[0063] It should be understood that although the steps in the flowcharts of the various embodiments of the present invention are shown sequentially according to the arrows, these steps are not necessarily executed in the order indicated by the arrows. Unless explicitly stated herein, there is no strict order restriction on the execution of these steps, and they can be executed in other orders. Moreover, at least some steps in the various embodiments may include multiple sub-steps or multiple stages. These sub-steps or stages are not necessarily completed at the same time, but can be executed at different times. The execution order of these sub-steps or stages is not necessarily sequential, but can be performed alternately or in turn with other steps or at least a portion of the sub-steps or stages of other steps.
[0064] Those skilled in the art will understand that all or part of the processes in the methods of the above embodiments can be implemented by a computer program instructing related hardware. The program can be stored in a non-volatile computer-readable storage medium, and when executed, it can include the processes of the embodiments of the above methods. Any references to memory, storage, databases, or other media used in the embodiments provided in this application can include non-volatile and / or volatile memory. Non-volatile memory can include read-only memory (ROM), programmable ROM (PROM), electrically programmable ROM (EPROM), electrically erasable programmable ROM (EEPROM), or flash memory. Volatile memory can include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM is available in various forms, such as static RAM (SRAM), dynamic RAM (DRAM), synchronous DRAM (SDRAM), dual data rate SDRAM (DDRSDRAM), enhanced SDRAM (ESDRAM), synchronous link DRAM (SLDRAM), RAMbus direct RAM (RDRAM), direct memory bus dynamic RAM (DRDRAM), and RAMbus dynamic RAM (RDRAM), etc.
[0065] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.
[0066] The above-described embodiments are merely illustrative of several implementations of the present invention, and while the descriptions are specific and detailed, they should not be construed as limiting the scope of the present invention. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these modifications and improvements all fall within the scope of protection of the present invention. Therefore, the scope of protection of this patent should be determined by the appended claims.
[0067] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.
Claims
1. A method for accelerating computation of high-resolution multispectral data PCA dimensionality reduction algorithm based on heterogeneous computing, characterized in that, include: S1. Realize the reading and loading of multi-dimensional data volumes of high-resolution multispectral remote sensing data; S2. Solve for the self-weighted average of non-background pixels in single-dimensional data using a heterogeneous computation method based on the cyclic split window algorithm; S3. Achieve center normalization of non-background pixel data in single-dimensional data; S4. Solve for the covariance matrix between different dimensions of the data volume using a heterogeneous computation method; S5. Solve for the eigenvalues and corresponding eigenvectors of the covariance matrix using the Kevens rotation matrix iteration method; S6. Solve for the projection of multi-dimensional pixel vectors onto each principal component; S7. Acquire and store the result data of several principal component bands after projection calculation; The cyclic split window algorithm includes: Based on single-dimensional data of a certain band in a multidimensional data volume of high-resolution multispectral remote sensing data, the image data scale of the single-dimensional data in the X direction is represented as XSize, and the image data scale of the single-dimensional data in the Y direction is represented as YSize; within the range of the single-dimensional data, the size of the image data window is set, and the image data scale of the image data window in the X direction is represented as Winx, and the image data scale of the image data window in the Y direction is represented as Winy; a split formula is used to scale down the single-dimensional data. For the pixels in the scaled-down single-dimensional data, each pixel is the sum of its corresponding non-background pixels in the single-dimensional data before scaling down. The split formula is described as follows: Where XSize' and YSize' are the scales of the scaled single-dimensional data in the X and Y directions, respectively; XSize and YSize are the scales of the single-dimensional data in the X and Y directions before scaling down; Winx and Winy are the scales of the pixel data window in the single-dimensional data in the X and Y directions, respectively; and the process of calculating the individual pixel value of the scaled single-dimensional data is put into GPGPU as a parallel computing unit for processing. When the scale of a single-dimensional data point in the X or Y direction is insufficient to be divided by the scale of the image data window in the X or Y direction, the single-dimensional data is supplemented with corresponding 0-cell columns and 0-cell rows. The number of columns (col) and rows (row) to be supplemented are calculated using the following formula: Where Winx and Winy are the scales of the image metadata window in the X and Y directions, respectively, and XSize and YSize are the scales of the single-dimensional data in the X and Y directions, respectively. Repeat the scaling down process on the single-dimensional data until a preset number of processing steps are reached, or the scale of the single-dimensional data in both the X and Y directions has been reduced to below the preset scale, at which point the scaling down process stops.
2. The accelerated computation method for high-resolution multispectral data PCA dimensionality reduction algorithm based on heterogeneous computing according to claim 1, characterized in that, Heterogeneous computing methods delegate the logic control, data access process, and low-computation processing involved in the PCA dimensionality reduction algorithm to the CPU, while delegating massive pixel-level graphics processing to the GPU.
3. The accelerated computation method for high-resolution multispectral data PCA dimensionality reduction algorithm based on heterogeneous computing according to claim 2, characterized in that, S2 solves for the self-weighted average of non-background pixels in one-dimensional data, including: obtaining B using the cyclic split window algorithm. k Sum(DataB) of all pixels in a single-dimensional data band k ); B k The value of each cell in the single-dimensional data of the band, DataB k Substituting (i, j) into the factor calculation formula, we obtain the intermediate calculation factor W. t (i, j), the factor calculation formula is described as follows: B is obtained through the formula for calculating the self-proportionally weighted average. k The self-weighted average of single-dimensional data in the band, Avg(DataB) k The formula for calculating the self-proportionally weighted average is described below: Where xsize and ysize are B respectively k The scale of single-dimensional data in the X and Y directions of the band.
4. The accelerated computation method for high-resolution multispectral data PCA dimensionality reduction algorithm based on heterogeneous computing according to claim 3, characterized in that, S3 achieves the central normalization of non-background image data of single-dimensional data through the following calculation formula: Among them, CenDataB k (i,j) is B k The pixel values of the single-dimensional data of the band after center normalization, DataB k (i,j) is B k The pixel values of a single-dimensional data band, Avg(DataB) k ) is B k The self-weighted average of single-dimensional data for a band.
5. The accelerated computation method for high-resolution multispectral data PCA dimensionality reduction algorithm based on heterogeneous computing according to claim 1, characterized in that, S4 solves for the covariance matrix between each dimension of the data volume, including: centering any two bands b of the data volume. m ,b n The intermediate factor TmpDataB is obtained by sequentially substituting all pixels of the single-dimensional data into the intermediate factor calculation formula. mn (i,j), the intermediate factor calculation formula is described as TmpDataB mn (i,j)=CenDataB m (i, j)×CenDataB n (i,j), where CenDataB m (i, j), CenDataB n (i,j) represent band b respectively. m ,b n The values of the single-dimensional data after center standardization; the sum M of all intermediate computational factors is calculated, and the calculation process of M is specifically expressed as follows: Where xsize and ysize are respectively b m ,b n The scale of single-dimensional data in the X and Y directions of the band; Calculate band b m ,b n Each pixel in the single-dimensional binary data DataB bin The sum of (i, j) and N, specifically the calculation process of N is expressed as follows: Where xsize and ysize are respectively b m ,b n The scale of the single-dimensional binary data of the band in the X and Y directions; Through calculation Obtain band b m ,b n The covariance of single-dimensional data; Repeat the above steps until all elements of the covariance matrix have been calculated.