SAR interferometric phase denoising method and device based on BM3D

Abstract

The present application belongs to the technical field of interferometric phase map processing, and relates to a SAR interferometric phase denoising method and device based on BM3D. A two-dimensional thread grid at a pixel level is used to realize parallel implementation of similar block position index set screening, maximum likelihood estimation of pixels and decorrelation transformation. A two-dimensional thread grid at a block level is used to realize parallel implementation of biorthogonal wavelet transformation of image blocks. A three-dimensional thread grid is used for three-dimensional block group stacking, and a two-dimensional thread grid at a block level is used to realize parallel implementation of Hadamard transformation of the third dimension of the three-dimensional block group. A three-dimensional thread grid is used for hard threshold filtering of transformation coefficients in the three-dimensional block group. A two-dimensional thread grid at a block level is used to realize parallel implementation of inverse transformation of image blocks, a first layered two-dimensional thread grid is used to realize parallel implementation of weighted aggregation of spatial domain pixels, a two-dimensional thread grid at a pixel level is used to realize parallel implementation of inverse decorrelation transformation of spatial domain pixels, and a complex interferometric image after preliminary denoising is obtained, and then secondary filtering is performed to obtain a target complex interferometric image.

Description

Technical Field

[0001] This invention relates to the field of interferometric phase image processing technology, and in particular to a SAR interferometric phase denoising method and apparatus based on BM3D. Background Technology

[0002] Synthetic Aperture Radar (SAR) is an active imaging system capable of operating in all weather conditions. Through interferometric Synthetic Aperture Radar (InSAR), it can obtain high-precision information on surface elevation or deformation. However, due to speckle noise and system thermal noise, a large amount of random noise exists in the interferometric phase map, which not only degrades the quality of the interferometric fringes but also interferes with subsequent phase unwrapping and deformation extraction. Therefore, effectively suppressing noise while preserving the details of the interferometric fringes has always been a key issue in InSAR data processing.

[0003] Among numerous phase denoising algorithms, the BM3D (Block-Matching and 3D Filtering) algorithm is often overlooked in InSAR phase denoising due to its excellent ability to utilize structural redundancy. It achieves a good balance between noise suppression and detail preservation. This algorithm stacks similar image blocks into three-dimensional block groups through non-local block matching and performs collaborative filtering in the three-dimensional transform domain. However, the computational complexity of core steps such as block matching, three-dimensional transformation, and weighted aggregation is extremely high. When implemented serially, it requires a large number of loops and matrix operations, which takes too long to process high-resolution interferometric images and cannot meet the high real-time requirements of interferometric image denoising.

[0004] In summary, existing interferometric image denoising methods using the BM3D algorithm suffer from high computational complexity and excessive processing time, failing to meet high real-time performance requirements. Summary of the Invention

[0005] Therefore, the technical problem to be solved by the present invention is to overcome the problems of high computational complexity and long processing time in the existing interferometric image denoising methods using the BM3D algorithm, which cannot meet the requirements of high real-time performance.

[0006] To address the aforementioned technical problems, this invention provides a SAR interferometric phase denoising method based on BM3D, comprising: Construct pixel-level two-dimensional threaded grids with dimensions of image width and image height, block-level two-dimensional threaded grids with dimensions of image width and search window width, block-level transformation two-dimensional threaded grids with dimensions of image width and number of pixels in image blocks, first-level layered two-dimensional threaded grids with dimensions of image width and number of similar blocks, and three-dimensional threaded grids with dimensions of image width, number of similar blocks and number of pixels in similar blocks. In a pixel-level two-dimensional threaded grid, each thread simultaneously traverses the candidate blocks within the search window of its corresponding pixel, calculates the similarity between the candidate block and the reference block centered on its corresponding pixel, filters to obtain a set of similar block position indices, and performs maximum likelihood estimation and decorrelation transformation on its corresponding pixels to obtain the real and imaginary components of its corresponding pixels. In a block-level two-dimensional thread grid, each thread simultaneously performs a two-dimensional bioorthogonal wavelet transform on its corresponding image block to obtain the transform coefficient set of its corresponding image block. In the 3D thread mesh, each thread stacks the transformation coefficients of each pixel according to the set of similar block position indices to obtain the real part 3D block group and the imaginary part 3D block group; in the block-level transformation 2D thread mesh, each thread simultaneously performs a one-dimensional fast Hadamard transformation on the third dimension of the 3D block group where its corresponding pixel is located; in the 3D thread mesh, each thread performs hard thresholding filtering on the transformation coefficients of its corresponding pixel in the transformed 3D block group. In the block-level two-dimensional thread grid, each thread simultaneously performs an inverse transform on its corresponding image block after hard thresholding to obtain the spatial domain group of its corresponding image block; in the first-level two-dimensional thread grid, each thread simultaneously performs weighted aggregation on the spatial domain pixels of its corresponding pixel and returns it to the complex interference image. In a pixel-level two-dimensional thread grid, each thread simultaneously performs an inverse decorrelation transformation on the spatial domain pixels after weighted aggregation of its corresponding pixels to obtain a preliminary denoised complex interference image. The target complex interferometric image is obtained by performing secondary filtering on the initially denoised complex interferometric image.

[0007] Preferably, the target complex interferometric image is obtained by performing a second filtering on the initially denoised complex interferometric image, including: In the pixel-level two-dimensional thread grid, each thread simultaneously calculates the similarity between the candidate block and the reference block centered on its corresponding pixel in the search window of the candidate block in the complex interferometric image after preliminary denoising. A new set of similar block position indices is obtained by filtering, and maximum likelihood estimation and decorrelation transformation are performed on the corresponding pixels to obtain the real and imaginary components of the corresponding pixels after preliminary denoising. A second-layer two-dimensional threaded grid with dimensions of image width and search window width is constructed. Each thread in each thread block of the second-layer two-dimensional threaded grid simultaneously performs two-dimensional discrete cosine transform on the real and imaginary components of each pixel in the corresponding image block of each thread block, as well as the real and imaginary components after preliminary denoising, to obtain the transform coefficient set of each image block. In the 3D thread mesh, each thread stacks the transformation coefficients of each pixel according to the new set of similar block position indices based on the transformation coefficient group of the image block, to obtain a new real part 3D block group and an imaginary part 3D block group, and performs Wiener filtering on the corresponding transformation coefficients in the new real part 3D block group and imaginary part 3D block group. In the second-layer thread grid, each thread block simultaneously performs an inverse transform on its corresponding image block after Wiener filtering to obtain the spatial domain group of its corresponding image block after Wiener filtering; in the first-layer two-dimensional thread grid, each thread simultaneously performs weighted aggregation on the spatial domain pixels of its corresponding pixels after Wiener filtering and returns to the complex interference image after preliminary denoising. In a pixel-level two-dimensional threaded grid, each thread simultaneously performs an inverse decorrelation transformation on the spatial domain pixels of its corresponding pixels in the pre-denoised complex interferometric image after weighted aggregation, to obtain the target complex interferometric image.

[0008] Preferably, the real part of the Wiener shrinkage coefficient of the pixel during Wiener filtering is... The calculation formula is: , in, Represents the transformation result matrix; This represents the estimated coherence coefficient of a pixel. Wiener shrinkage coefficient of the imaginary part of a pixel during Wiener filtering The calculation formula is: .

[0009] Preferably, when weighting and aggregating the spatial domain pixels corresponding to the Wiener-filtered pixels, the real part weight of the pixel is... The calculation formula is: , The imaginary part weight of the pixel when performing weighted aggregation of the spatial domain pixels corresponding to the Wiener filter. The calculation formula is: .

[0010] Preferably, in a pixel-level two-dimensional threaded grid, each thread simultaneously traverses the candidate blocks within the search window of its corresponding pixel, calculates the similarity between the candidate block and the reference block centered on its corresponding pixel, and filters to obtain a set of similar block position indices, including: In a pixel-level two-dimensional thread grid, each thread calculates the cosine distance between each pixel of the reference block where its corresponding pixel is located and each pixel of each candidate block in its search window. The cosine distances between each pixel are summed to obtain the similarity between the reference block and each candidate block in its search window. In a pixel-level two-dimensional thread grid, each thread compares the similarity between the reference block containing its corresponding pixel and each candidate block in its search window with a preset threshold, and takes the candidate block with a similarity greater than the preset threshold as the potential similar block of the reference block. In a pixel-level two-dimensional threaded grid, each thread sorts the potential similar blocks of the reference block containing its corresponding pixel in descending order of similarity, and obtains the set of similar block position indices of the reference block based on the position indices of the first preset number of potential similar blocks.

[0011] Preferably, in a pixel-level two-dimensional threaded mesh, each thread performs maximum likelihood estimation and decorrelation transformation on its corresponding pixel to obtain the real and imaginary components of the corresponding pixel, including: In a pixel-level two-dimensional threaded mesh, each thread calculates the mean intensity of all pixels in the similar blocks of the reference block where its corresponding pixel is located, and obtains the amplitude squared estimate of the similar blocks of the reference block. Based on the amplitude squared estimate of all similar blocks of the reference block, the complex coherence coefficient of the set of similar blocks of the reference block is calculated, thereby obtaining the coherence coefficient estimate and phase estimate of the corresponding pixel. In a pixel-level two-dimensional threaded mesh, each thread uses the decorrelation transformation model and the phase estimate of its corresponding pixel to perform orthogonal component decomposition on the complex interference signal of its corresponding pixel, thereby obtaining the real and imaginary components of the pixel.

[0012] Preferably, each thread in the 3D thread mesh performs hard thresholding filtering on the transformation coefficients of its corresponding pixels in the transformed 3D block group, including: In the 3D thread mesh, each thread calculates the real part filtering threshold and the imaginary part filtering threshold of its corresponding pixel based on the estimated coherence coefficient of that pixel. In the 3D thread mesh, each thread compares the real part 3D transformation coefficient of the corresponding pixel in the transformed real part 3D block group with the real part filtering threshold of the pixel. If the real part 3D transformation coefficient is greater than the real part filtering threshold of the pixel, the real part 3D transformation coefficient of the pixel remains unchanged. If the real part 3D transformation coefficient is less than or equal to the real part filtering threshold of the pixel, the real part 3D transformation coefficient of the pixel is set to zero. In the 3D thread mesh, each thread compares the virtual part 3D transformation coefficient of the corresponding pixel in the transformed virtual part 3D block group with the virtual part filtering threshold of the pixel. If the virtual part 3D transformation coefficient is greater than the virtual part filtering threshold of the pixel, the virtual part 3D transformation coefficient of the pixel remains unchanged. If the virtual part 3D transformation coefficient is less than or equal to the virtual part filtering threshold of the pixel, the virtual part 3D transformation coefficient of the pixel is set to zero. In the 3D thread mesh, each thread calculates the real part weight of the pixel based on the number of pixels in the similar block to which its corresponding pixel has a non-zero real part 3D transform coefficient after hard thresholding and the estimated coherence coefficient of that pixel. In the 3D thread mesh, each thread calculates the imaginary part weight of a pixel based on the number of pixels in the similar block to which its corresponding pixel has a non-zero imaginary part 3D transform coefficient after hard thresholding and the estimated coherence coefficient of that pixel.

[0013] Preferably, in the first layered two-dimensional thread grid, each thread simultaneously performs weighted aggregation of the spatial domain pixels of its corresponding pixel and returns it to the complex interferometric image, including: In the first layer of the two-dimensional threaded mesh, each thread calculates the product of the real part spatial domain and the real part weight of its corresponding pixel to obtain the weighted real part pixel value of that pixel. In the first layer of the two-dimensional threaded grid, each thread calculates the product of the imaginary part spatial domain and the imaginary part weight of its corresponding pixel to obtain the weighted imaginary part pixel value of that pixel. In the first layer of the two-dimensional threaded grid, each thread uses atomic addition to accumulate the weighted real part pixel value and the weighted imaginary part pixel value of its corresponding pixel in the three-dimensional block group to the corresponding position in the complex interference image, thereby obtaining the real part pixel value and the imaginary part pixel value of that pixel.

[0014] Preferably, the real part filtering threshold of the pixel The calculation formula is: , in, It is a constant; This represents the estimated coherence coefficient of a pixel. Imaginary part filtering threshold of a pixel The calculation formula is: , Real part weight of a pixel The calculation formula is: , in, This represents the number of pixels in the similar block to which the pixel belongs that have non-zero real-part three-dimensional transformation coefficients after hard thresholding filtering; Imaginary part weight of a pixel The calculation formula is: , in, This represents the number of pixels in the similar block to which the pixel belongs whose imaginary 3D transformation coefficients are non-zero after hard thresholding.

[0015] The present invention also provides a SAR interferometric phase denoising device based on BM3D, comprising: The thread mesh construction module is used to construct pixel-level two-dimensional thread meshes with dimensions of image width and image height, block-level two-dimensional thread meshes with dimensions of image width and search window width, block-level transformation two-dimensional thread meshes with dimensions of image width and number of pixels in image blocks, first-level layered two-dimensional thread meshes with dimensions of image width and number of similar blocks, and three-dimensional thread meshes with dimensions of image width, number of similar blocks and number of pixels in similar blocks. The similar block filtering and related transformation module is used to utilize each thread in the pixel-level two-dimensional thread grid to simultaneously traverse the candidate blocks in the search window where the corresponding pixel is located, calculate the similarity between the candidate block and the reference block centered on the corresponding pixel, filter to obtain the set of similar block position indices, and perform maximum likelihood estimation and decorrelation transformation on the corresponding pixels to obtain the real and imaginary components of the corresponding pixels. The transform coefficient acquisition module is used to simultaneously perform two-dimensional bioorthogonal wavelet transform on the corresponding image blocks by each thread in the block-level two-dimensional thread grid, and obtain the transform coefficient group of the corresponding image blocks. The collaborative filtering module is used to stack the transformation coefficients of each pixel according to the set of similar block position indices in the three-dimensional thread mesh to obtain the real part three-dimensional block group and the imaginary part three-dimensional block group; the block-level transformation module is used to simultaneously perform one-dimensional fast Hadamard transformation on the third dimension of the three-dimensional block group where the corresponding pixel is located in the two-dimensional thread mesh; and the three-dimensional thread mesh is used to perform hard thresholding filtering on the transformation coefficients of the corresponding pixels in the transformed three-dimensional block group. The inverse transformation and weighted aggregation module is used to simultaneously perform inverse transformation on the corresponding image blocks after hard thresholding by each thread in the block-level two-dimensional thread grid to obtain the spatial domain group of the corresponding image blocks; each thread in the first layer two-dimensional thread grid simultaneously performs weighted aggregation on the spatial domain pixels of its corresponding pixels and returns it to the complex interference image. The preliminary denoised image acquisition module is used to perform inverse decorrelation transformation on the spatial domain pixels after weighted aggregation of corresponding pixels by each thread in the pixel-level two-dimensional thread grid, so as to obtain the complex interference image after preliminary denoising. The secondary filtering module is used to perform secondary filtering on the initially denoised complex interferometric image to obtain the target complex interferometric image.

[0016] The SAR interferometric phase denoising method based on BM3D provided in this application has the following beneficial effects: The BM3D algorithm was modified in stages using a GPU heterogeneous parallel computing architecture. Pixel coordinates were mapped to threads in a pixel-level 2D thread grid. All threads simultaneously traversed different reference blocks and candidate blocks. Since each thread was bound to a single pixel, block-by-block searching was unnecessary during similarity block selection. Each thread then independently processed the maximum likelihood estimation and decorrelation transformation of a single pixel, enabling parallel implementation of the transformation operations. Image blocks were mapped to threads in a block-level 2D thread grid. All threads simultaneously performed 2D bioorthogonal wavelet transform on their corresponding image blocks, obtaining the transform coefficient set for each block. In a 3D thread grid, the transform coefficients of each pixel were stacked according to the similarity block position index set, resulting in a real part 3D block set and an imaginary part 3D block set. In the block-level transformation 2D thread grid, each thread simultaneously performed a 1D fast Hadamard transform on the third dimension of the 3D block set containing its corresponding pixel. Finally, in the 3D thread grid, each thread performed hard thresholding filtering on the transform coefficients of the corresponding pixels in the transformed 3D block set. After wavelet transformation, each thread in the block-level two-line threaded grid simultaneously performs an inverse transformation on its corresponding image block to obtain the spatial domain group of the corresponding image block. Then, each thread in the first-level two-dimensional threaded grid simultaneously performs weighted aggregation on the spatial domain pixels of its corresponding pixels and returns it to the complex interferometric image. Then, each thread in the pixel-level two-dimensional threaded grid simultaneously performs an inverse decorrelation transformation on the spatial domain pixels of the weighted aggregation of its corresponding pixels to obtain the preliminary denoised complex interferometric image. Thus, a second filtering is performed on the preliminary denoised complex interferometric image to obtain the target complex interferometric image. This application adopts a thread-level and block-level collaborative mechanism under the CUDA framework to implement computationally complex modules such as decorrelation transformation, block matching, three-dimensional collaborative filtering, and weighted fusion in parallel. This significantly improves the denoising efficiency of complex interferometric images while still effectively suppressing noise in low-coherence regions and maintaining the structural continuity and detail integrity of interference fringes, achieving a balance between denoising performance and accuracy. Therefore, it can meet the high real-time requirements of interferometric image denoising while ensuring denoising accuracy. Attached Figure Description

[0017] To make the content of this invention easier to understand, the invention will be further described in detail below with reference to specific embodiments and accompanying drawings, wherein: Figure 1 Flowchart of the BM3D-based SAR interferometric phase denoising method provided for this application; Figure 2 A schematic diagram of the pixel-level parallel processing structure in the BM3D-based SAR interferometric phase denoising method provided in this application; Figure 3 A schematic diagram of block-level parallel two-dimensional bioorthogonal wavelet transform in the BM3D-based SAR interferometric phase denoising method provided in this application; Figure 4 A schematic diagram of pixel-level parallel stacking of three-dimensional block groups of three-dimensional thread mesh and thread blocks in the BM3D-based SAR interferometric phase denoising method provided in this application; Figure 5 This application provides a schematic diagram of the parallel one-dimensional Hadamard transform in the BM3D-based SAR interferometric phase denoising method; Figure 6 A schematic diagram of pixel-level parallel hard thresholding filtering of three-dimensional thread mesh and thread block in the BM3D-based SAR interferometric phase denoising method provided in this application; Figure 7 A schematic diagram of parallel weighted aggregation in the BM3D-based SAR interferometric phase denoising method provided in this application; Figure 8 A schematic diagram of atomic addition in the BM3D-based SAR interferometric phase denoising method provided in this application; Figure 9 A schematic diagram of pixel-level parallel Wiener filtering of three-dimensional thread mesh and thread blocks in the BM3D-based SAR interferometric phase denoising method provided in this application; Figure 10 This is a schematic diagram comparing the interference phase recovery results obtained by the CPU and GPU for simulated interferograms, provided in an embodiment of this application; wherein, Figure 10 (a) in the diagram is a noise-free simulated interferogram. Figure 10 (b) in the figure is a simulated interferogram with noise. Figure 10 (c) shows the denoising results of the InSAR-BM3D algorithm implemented in parallel using OpenMP. Figure 10 (d) in the figure represents the denoising result of the InSAR-BM3D algorithm implemented in parallel using CUDA; Figure 11 This is a time comparison diagram of the interference phase grayscale images obtained using OpenMP and CUDA for simulated interferograms of different sizes provided in the embodiments of this application. Detailed Implementation

[0018] The present invention will be further described below with reference to the accompanying drawings and specific embodiments, so that those skilled in the art can better understand and implement the present invention. However, the embodiments described are not intended to limit the present invention.

[0019] Please see Figure 1 , Figure 1 The flowchart shown is of the SAR interferometric denoising method based on BM3D provided in this application. The method specifically includes steps S10 to S70: S10: Construct a pixel-level two-dimensional threaded grid with dimensions of image width and image height, a block-level two-dimensional threaded grid with dimensions of image width and search window width, a block-level transformation two-dimensional threaded grid with dimensions of image width and number of pixels in image blocks, a first-level two-dimensional threaded grid with dimensions of image width and number of similar blocks, and a three-dimensional threaded grid with dimensions of image width, number of similar blocks, and number of pixels in similar blocks.

[0020] S20: In a pixel-level two-dimensional threaded grid, each thread simultaneously traverses the candidate blocks within the search window of its corresponding pixel, calculates the similarity between the candidate block and the reference block centered on its corresponding pixel, filters to obtain a set of similar block position indices, and performs maximum likelihood estimation and decorrelation transformation on its corresponding pixels to obtain the real and imaginary components of its corresponding pixels.

[0021] For example, such as Figure 2 As shown, this application transmits the complex interferometric image to the global memory of the GPU and adopts a pixel-level parallel strategy based on image space on the GPU side to achieve parallel processing of the entire complex interferometric image. Specifically, the computation is organized using a two-dimensional grid and two-dimensional thread block structure. The mapping from pixel-level spatial coordinates to two-dimensional thread indices is completed through blockldx and threadldx, so that a large number of reference pixels in the image can be executed in parallel on the GPU at the same time.

[0022] Specifically, step S20 includes S200~S202: S200: In a pixel-level two-dimensional threaded grid, each thread calculates the cosine distance between each pixel of the reference block where its corresponding pixel is located and each pixel of each candidate block in its search window. The cosine distances between each pixel are accumulated to obtain the similarity between the reference block and each candidate block in its search window.

[0023] S201: In a pixel-level two-dimensional thread grid, each thread compares the similarity between its corresponding pixel's reference block and each candidate block in its search window with a preset threshold, and takes the candidate block with a similarity greater than the preset threshold as the potential similar block of the reference block.

[0024] S202: In the pixel-level two-dimensional thread grid, each thread sorts the potential similar blocks of the reference block where its corresponding pixel is located in descending order of similarity, and obtains the set of similar block position indices of the reference block based on the position indices of the first preset number of potential similar blocks.

[0025] In one specific implementation, after sorting the potential similar blocks, the top 64 potential similar blocks are obtained as similar blocks of the reference block, thereby obtaining the set of similar blocks of the reference block, and the set of similar block position indices of the reference block is written into the GPU global memory.

[0026] like Figure 2 As shown, this application adopts a pixel-level parallel processing architecture when filtering the set of similar blocks for each reference block. Each thread in the pixel-level two-dimensional thread grid independently calculates the cosine distance between the corresponding pixel in the reference block and all candidate blocks in its search window, and filters out the set of 64 most similar block indices.

[0027] Furthermore, step S20 also includes S203~S204: S203: In a pixel-level two-dimensional threaded mesh, each thread calculates the mean intensity of all pixels in the similar blocks of the reference block where its corresponding pixel is located, and obtains the amplitude squared estimate of the similar blocks of the reference block. Based on the amplitude squared estimates of all similar blocks of the reference block, the complex coherence coefficient of the set of similar blocks of the reference block is calculated, thereby obtaining the coherence coefficient estimate and phase estimate of the corresponding pixel.

[0028] Specifically, the magnitude squared estimate of similar blocks The calculation formula is: , in, This indicates the number of pixels in the similar block; , This represents the complex pixel value of the i-th pixel in the similar block of the two registered SAR images corresponding to the complex interferometric image.

[0029] Complex coherence coefficients of similar block sets The calculation is as follows: , in, This indicates the number of similar blocks in the set of similar blocks.

[0030] Pixel coherence coefficient estimate The calculation formula is: .

[0031] Pixel phase estimate The calculation formula is: .

[0032] The coherence coefficient estimate is used to calculate the noise variance and the filtering threshold, while the phase estimate is used for the decorrelation transformation of the noise.

[0033] S204: In a pixel-level two-dimensional threaded mesh, each thread uses the decorrelation transformation model and the phase estimate of its corresponding pixel to perform orthogonal component decomposition on the complex interference signal of its corresponding pixel, thereby obtaining the real and imaginary components of the pixel.

[0034] Specifically, the real and imaginary components of a pixel are represented as follows: , , in, , Let these represent the real and imaginary components after decorrelation transformation, respectively; , These represent the real and imaginary components of the complex interference signal before decorrelation transformation, respectively.

[0035] S30: In a block-level two-dimensional thread grid, each thread simultaneously performs a two-dimensional biorthogonal wavelet transform on its corresponding image block to obtain the transform coefficient set of its corresponding image block.

[0036] like Figure 3 As shown, each thread in the block-level two-dimensional thread grid is responsible for calculating the two-dimensional bioorthogonal wavelet transform of an image block, enabling all image blocks within the search window containing a row of pixels to perform the two-dimensional transform in parallel.

[0037] S40: In the 3D threaded mesh, each thread stacks the transformation coefficients of each pixel according to the set of similar block position indices to obtain the real part 3D block group and the imaginary part 3D block group; in the block-level transformation 2D threaded mesh, each thread simultaneously performs a one-dimensional fast Hadamard transformation on the third dimension of the 3D block group where its corresponding pixel is located; in the 3D threaded mesh, each thread performs hard thresholding filtering on the transformation coefficients of its corresponding pixel in the transformed 3D block group.

[0038] Specifically, when constructing 3D block groups in parallel, this application adopts a GPU-based 3D threaded mesh parallel architecture, the parallel processing principle of which is as follows: Figure 4 As shown, the three-dimensional thread mesh contains multiple threads. Each thread corresponds to a pixel in a similar block in the similar block set of each reference block in each reference block column. That is, each thread is only responsible for processing a single pixel position. Based on its own thread index and the similar block set of each reference block, it reads the pixel value corresponding to a single pixel from the two-dimensional transformation coefficients and writes it into a three-dimensional array, thereby constructing a stacked three-dimensional block group.

[0039] Furthermore, after stacking, a one-dimensional Hadamard transformation needs to be performed on the third dimension of the three-dimensional block group, such as... Figure 5 As shown, this application designs a two-dimensional thread mesh for block group transformation. There are 64 data points at each two-dimensional position of the three-dimensional block group. Each 64 data points form a one-dimensional vector. First, this one-dimensional vector is transferred to the shared memory of the corresponding GPU thread block. Each thread is responsible for the butterfly operation of one transformation coefficient in this one-dimensional vector and stores the operation result in the shared memory. After all threads in a thread block have completed the calculation, the transformation coefficient in the shared memory is transferred to the global memory.

[0040] Furthermore, such as Figure 6 As shown, each thread in the 3D thread mesh performs hard thresholding filtering on the transformation coefficients of its corresponding pixels in the transformed 3D block group, including S400~S404: S400: For each thread in the 3D threaded mesh, calculate the real part filtering threshold and the imaginary part filtering threshold for that pixel based on the estimated coherence coefficient of that pixel.

[0041] Specifically, the real part filtering threshold of a pixel The calculation formula is: , in, It is a constant, usually taken as 4.2; Represents the estimated coherence coefficient of a pixel; Imaginary part filtering threshold of a pixel The calculation formula is: .

[0042] S401: In the 3D thread mesh, each thread compares the real part 3D transformation coefficient of the corresponding pixel in the transformed real part 3D block group with the real part filtering threshold of the pixel. If the real part 3D transformation coefficient is greater than the real part filtering threshold of the pixel, the real part 3D transformation coefficient of the pixel remains unchanged. If the real part 3D transformation coefficient is less than or equal to the real part filtering threshold of the pixel, the real part 3D transformation coefficient of the pixel is set to zero.

[0043] S402: In the 3D thread mesh, each thread compares the virtual part 3D transformation coefficient of the corresponding pixel in the transformed virtual part 3D block group with the virtual part filtering threshold of the pixel. If the virtual part 3D transformation coefficient is greater than the virtual part filtering threshold of the pixel, the virtual part 3D transformation coefficient of the pixel remains unchanged. If the virtual part 3D transformation coefficient is less than or equal to the virtual part filtering threshold of the pixel, the virtual part 3D transformation coefficient of the pixel is set to zero.

[0044] S403: In the 3D thread mesh, each thread calculates the real part weight of the pixel based on the number of pixels in the similar block to which its corresponding pixel has a non-zero real part 3D transform coefficient after hard thresholding and the estimated coherence coefficient of that pixel. S404: In a 3D threaded mesh, each thread calculates the imaginary part weight of a pixel based on the number of pixels in the similar block to which its corresponding pixel has a non-zero imaginary 3D transform coefficient after hard thresholding and the estimated coherence coefficient of that pixel.

[0045] Specifically, the real part weight of the pixel The calculation formula is: , in, This represents the number of pixels in the similar block to which the pixel belongs that have non-zero real-part three-dimensional transformation coefficients after hard thresholding filtering; Imaginary part weight of a pixel The calculation formula is: , in, This represents the number of pixels in the similar block to which the pixel belongs whose imaginary 3D transformation coefficients are non-zero after hard thresholding.

[0046] S50: In the block-level two-dimensional thread grid, each thread simultaneously performs an inverse transformation on its corresponding image block after hard thresholding to obtain the spatial domain group of its corresponding image block; In the first layer two-dimensional thread grid, each thread simultaneously performs weighted aggregation on the spatial domain pixels of its corresponding pixel and returns it to the complex interference image.

[0047] Furthermore, such as Figure 7 As shown, in the first layered two-dimensional thread grid, each thread simultaneously performs weighted aggregation on the spatial domain pixels of its corresponding pixel and returns it to the complex interferometric image, including S500~S502: S500: In the first layer of two-dimensional threaded mesh, each thread calculates the product of the real part spatial domain and the real part weight of its corresponding pixel to obtain the weighted real part pixel value of that pixel.

[0048] S501: In the first layer of two-dimensional threaded grid, each thread calculates the product of the imaginary spatial domain and the imaginary weight of its corresponding pixel to obtain the weighted imaginary pixel value of that pixel.

[0049] S502: In the first layer of two-dimensional thread grid, each thread uses atomic addition to accumulate the weighted real part pixel value and weighted imaginary part pixel value of its corresponding pixel in the three-dimensional block group to the corresponding position in the complex interference image, so as to obtain the real part pixel value and imaginary part pixel value of the pixel.

[0050] Specifically, the implementation of atomic addition is as follows: Figure 8 As shown, atomic addition is a technique used in parallel computing to resolve race conditions that arise when multiple threads simultaneously write data to the same memory address. Its principle is to use hardware-provided atomic instructions to ensure that multiple threads can correctly accumulate the same global variable simultaneously, avoiding data races and unpredictable results. Specifically, when a thread performs an atomic addition operation, the memory address is locked by hardware, allowing only one thread to access it at a time, and ensuring that the operation order is atomic. This guarantees that no race conditions occur when multiple threads accumulate data at the same memory address, improving the program's parallel performance.

[0051] S60: In a pixel-level two-dimensional threaded grid, each thread simultaneously performs an inverse decorrelation transformation on the spatial domain pixels after weighted aggregation of their corresponding pixels, resulting in a preliminary denoised complex interference image.

[0052] S70: Perform secondary filtering on the initially denoised complex interferometric image to obtain the target complex interferometric image.

[0053] Furthermore, step S70 specifically includes S700~S704: S700: In a pixel-level two-dimensional threaded grid, each thread simultaneously calculates the similarity between the candidate block and the reference block centered on its corresponding pixel in the search window of the candidate block in the complex interferometric image after preliminary denoising. A new set of similar block position indices is obtained by filtering, and maximum likelihood estimation and decorrelation transformation are performed on the corresponding pixels to obtain the real and imaginary components of the corresponding pixels after preliminary denoising.

[0054] S701: Construct a second-layer two-dimensional threaded grid with dimensions of image width and search window width. Each thread in each thread block of the second-layer two-dimensional threaded grid simultaneously performs two-dimensional discrete cosine transform on the real and imaginary components of each pixel in the corresponding image block of each thread block, as well as the real and imaginary components after preliminary denoising, to obtain the transform coefficient set of each image block.

[0055] S702: In the 3D thread mesh, each thread stacks the transform coefficients of each pixel according to the new set of similar block position indices based on the transform coefficient group of the image block, resulting in a new real part 3D block group and an imaginary part 3D block group. Wiener filtering is then applied to the corresponding transform coefficients in the new real part 3D block group and the imaginary part 3D block group. For example... Figure 9 The diagram shows a pixel-level parallel Wiener filter for a three-dimensional thread grid and thread blocks.

[0056] Specifically, the real part of the Wiener shrinkage coefficient of a pixel during Wiener filtering. The calculation formula is: , in, Represents the transformation result matrix; Represents the estimated coherence coefficient of a pixel; Wiener shrinkage coefficient of the imaginary part of a pixel during Wiener filtering The calculation formula is: .

[0057] S703: In the second-layer thread grid, each thread block simultaneously performs an inverse transformation on its corresponding image block after Wiener filtering to obtain the spatial domain group of its corresponding image block after Wiener filtering; In the first-layer two-dimensional thread grid, each thread simultaneously performs weighted aggregation on the spatial domain pixels of its corresponding pixels after Wiener filtering and returns to the complex interference image after preliminary denoising.

[0058] Specifically, the real part weight of the pixel when performing weighted aggregation on the spatial domain pixels corresponding to the Wiener filter. The calculation formula is: , The imaginary part weight of the pixel when performing weighted aggregation of the spatial domain pixels corresponding to the Wiener filter. The calculation formula is: .

[0059] S704: In a pixel-level two-dimensional threaded grid, each thread simultaneously performs an inverse decorrelation transformation on the spatial domain pixels of the corresponding pixels in the pre-denoised complex interferometric image after weighted aggregation, to obtain the target complex interferometric image.

[0060] It is worth noting that the two-dimensional discrete cosine transform of the second-stage filtering in this application uses image blocks as parallel processing units to achieve efficient parallel transformation of image blocks within the two-dimensional pixel region of the search window. That is, firstly, some threads in the thread block complete the one-dimensional discrete cosine transform in the row direction, and then some threads complete the one-dimensional discrete cosine transform in the column direction, thereby realizing the separate calculation of the two-dimensional discrete cosine transform.

[0061] For example, the mathematical definition and parallel implementation of the two-dimensional discrete cosine transform (DCT8x8) for 8×8 image blocks are closely centered around the core property of separability.

[0062] The transformation formula for two-dimensional DCT is shown below. However, due to the high computational complexity of direct calculation, in practical applications, its separability is utilized to decompose it into two one-dimensional DCTs. In matrix representation, this process can be expressed as: , , in, Represents the spatial domain pixel value; x and y represent the row and column indices of the pixels in the image patch; N represents the number of rows and columns of pixels in the image patch; Represents the frequency domain coefficients; u and v represent the frequency indices in the horizontal and vertical directions, respectively; Represents the transformation matrix; , These are the normalization coefficients; A matrix representing spatial pixels; Matrix A is in the form of: , in, , , , , , ; Furthermore, matrix A simplifies to: , Through this decomposition, the one-dimensional discrete cosine transform can be divided into independent operations of even and odd components, thus significantly reducing the number of multiplications. This simplified operational structure allows each thread to execute only a fixed sequence of additions and multiplications. The even component is represented as: , Odd components are represented as: .

[0063] The above-mentioned BM3D-based SAR interferometric phase denoising method will be explained and illustrated through a specific embodiment: Step 1: Obtain experimental data: 1. Obtain four 256×256 simulated interferometric phase maps. The four simulated interferometric phase maps correspond to different ground types. The coherence increases linearly from left to right, with the leftmost coherence value being 0.3 and the rightmost coherence value being 0.9.

[0064] 2. Obtain four simulated interferometric phase diagrams of different sizes, namely 256×256, 512×512, 1024×1024 and 5000×5000, to verify the applicability and stability of the above method under different data scales.

[0065] Step 2, Experiment Setup: In this embodiment, both the CPU and GPU use the same InSAR-BM3D algorithm flow and settings. The CPU uses OpenMP for multi-threaded parallel computation, while the GPU uses CUDA for large-scale parallel computation, ensuring the fairness and consistency of the comparison results. During implementation, the CPU and GPU use the same algorithm parameter configurations, including search window size, block matching parameters, and filtering parameters, to eliminate the influence of parameter differences on the experimental results. The peak signal-to-noise ratio (PSNR) is used as a quantitative evaluation index for the interferometric phase retrieval effect, measuring the degree of difference between the recovered phase and the original simulated phase.

[0066] The CPU-side parallel implementation is based on the Xeon Platinum 8375C processor and uses OpenMP serial technology; the GPU-side parallel implementation is based on the NVIDIA RTX 3090 graphics processor and uses the CUDA parallel computing framework.

[0067] like Figure 10 The figures shown are the interferometric phase retrieval results obtained by running the InSAR-BM3D algorithm on the CPU using OpenMP parallel mode and on the GPU using CUDA parallel mode for four simulated interferograms in this embodiment. Figure 10 (a) in the diagram is a noise-free simulated interferogram. Figure 10 (b) in the figure is a simulated interferogram with noise. Figure 10 (c) shows the denoising results of the InSAR-BM3D algorithm implemented in parallel using OpenMP. Figure 10 (d) in the figure shows the denoising results of the InSAR-BM3D algorithm implemented in parallel using CUDA. As can be seen from the figure, the phase recovery map obtained by the parallel implementation on the GPU side is highly consistent with the results of the parallel implementation on the CPU side in terms of overall structure, phase distribution and detail preservation. There are no obvious differences between the two, indicating that the SAR interferometric denoising method based on BM3D proposed in this application can accurately realize the algorithm function on the CPU side while ensuring the recovery accuracy.

[0068] Table 1 presents a comparison of the overall PSNR indices obtained from four simulated interferograms under different parallelization methods: Table 1

[0069] As can be seen from the data in Table 1, regardless of whether the CPU-side OpenMP parallel implementation or the GPU-side CUDA parallel implementation is used, the PSNR values ​​in each simulation scenario are relatively small. This further verifies that the GPU parallel implementation does not introduce additional accuracy loss and can achieve efficient parallel acceleration of the algorithm while maintaining almost the same denoising effect as the CPU implementation, thereby improving denoising efficiency.

[0070] like Figure 11The figure shows a time comparison of interferometric phase retrieval maps obtained using OpenMP and CUDA for simulated interferograms of different sizes provided in this embodiment. As can be seen from the figure, OpenMP's runtime is longer than CUDA's for all interferogram sizes. However, as the image size increases, the performance advantage of CUDA becomes increasingly apparent, exhibiting a significant acceleration effect under large-size interferogram conditions. Therefore, the high-performance SAR interferometric phase denoising method based on BM3D proposed in this invention has superior computational efficiency and better scalability when processing large-scale interferometric data, effectively meeting the rapid processing requirements of high-resolution, large-size SAR interferometric data.

[0071] Based on the BM3D-based SAR interferometric phase denoising method provided in the above embodiments, this application also provides a BM3D-based SAR interferometric phase denoising device, which specifically includes: The thread mesh construction module is used to construct pixel-level two-dimensional thread meshes with dimensions of image width and image height, block-level two-dimensional thread meshes with dimensions of image width and search window width, block-level transformation two-dimensional thread meshes with dimensions of image width and number of pixels in image blocks, first-level layered two-dimensional thread meshes with dimensions of image width and number of similar blocks, and three-dimensional thread meshes with dimensions of image width, number of similar blocks, and number of pixels in similar blocks.

[0072] The similar block filtering and related transformation module is used to simultaneously traverse the candidate blocks in the search window where each thread in the pixel-level two-dimensional thread grid is located, calculate the similarity between the candidate block and the reference block centered on its corresponding pixel, filter to obtain a set of similar block position indices, and perform maximum likelihood estimation and decorrelation transformation on its corresponding pixels to obtain the real and imaginary components of its corresponding pixels.

[0073] The transform coefficient acquisition module is used to simultaneously perform two-dimensional bioorthogonal wavelet transform on the corresponding image blocks by each thread in the block-level two-dimensional thread grid, so as to obtain the transform coefficient group of the corresponding image blocks.

[0074] The collaborative filtering module is used to stack the transformation coefficients of each pixel according to the set of similar block position indices in the three-dimensional thread grid to obtain the real part three-dimensional block group and the imaginary part three-dimensional block group; in the block-level transformation, each thread in the two-dimensional thread grid simultaneously performs a one-dimensional fast Hadamard transform on the third dimension of the three-dimensional block group where the corresponding pixel is located; and each thread in the three-dimensional thread grid performs hard thresholding filtering on the transformation coefficients of the corresponding pixels in the transformed three-dimensional block group.

[0075] The inverse transform and weighted aggregation module is used to simultaneously perform inverse transform on the corresponding image blocks after hard thresholding by each thread in the block-level two-dimensional thread grid to obtain the spatial domain group of the corresponding image blocks; each thread in the first layer two-dimensional thread grid simultaneously performs weighted aggregation on the spatial domain pixels of its corresponding pixels and returns it to the complex interference image.

[0076] The preliminary denoised image acquisition module is used to perform inverse decorrelation transformation on the spatial domain pixels after the corresponding pixels are weighted and aggregated by each thread in the pixel-level two-dimensional thread grid, so as to obtain the preliminary denoised complex interference image.

[0077] The secondary filtering module is used to perform secondary filtering on the initially denoised complex interferometric image to obtain the target complex interferometric image.

[0078] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0079] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart... Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0080] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0081] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0082] Obviously, the above embodiments are merely illustrative examples for clear explanation and are not intended to limit the implementation. Those skilled in the art will recognize that other variations or modifications can be made based on the above description. It is neither necessary nor possible to exhaustively list all possible implementations here. However, obvious variations or modifications derived therefrom are still within the scope of protection of this invention.

Claims

1. A SAR interferometric phase denoising method based on BM3D, characterized in that, include: Construct pixel-level two-dimensional threaded grids with dimensions of image width and image height, block-level two-dimensional threaded grids with dimensions of image width and search window width, block-level transformation two-dimensional threaded grids with dimensions of image width and number of pixels in image blocks, first-level layered two-dimensional threaded grids with dimensions of image width and number of similar blocks, and three-dimensional threaded grids with dimensions of image width, number of similar blocks and number of pixels in similar blocks. In a pixel-level two-dimensional threaded grid, each thread simultaneously traverses the candidate blocks within the search window of its corresponding pixel, calculates the similarity between the candidate block and the reference block centered on its corresponding pixel, filters to obtain a set of similar block position indices, and performs maximum likelihood estimation and decorrelation transformation on its corresponding pixels to obtain the real and imaginary components of its corresponding pixels. In a block-level two-dimensional thread grid, each thread simultaneously performs a two-dimensional bioorthogonal wavelet transform on its corresponding image block to obtain the transform coefficient set of its corresponding image block. In the 3D thread mesh, each thread stacks the transformation coefficients of each pixel according to the set of similar block position indices to obtain the real part 3D block group and the imaginary part 3D block group; in the block-level transformation 2D thread mesh, each thread simultaneously performs a one-dimensional fast Hadamard transformation on the third dimension of the 3D block group where its corresponding pixel is located; in the 3D thread mesh, each thread performs hard thresholding filtering on the transformation coefficients of its corresponding pixel in the transformed 3D block group. In the block-level two-dimensional thread grid, each thread simultaneously performs an inverse transform on its corresponding image block after hard thresholding to obtain the spatial domain group of its corresponding image block; in the first-level two-dimensional thread grid, each thread simultaneously performs weighted aggregation on the spatial domain pixels of its corresponding pixel and returns it to the complex interference image. In a pixel-level two-dimensional thread grid, each thread simultaneously performs an inverse decorrelation transformation on the spatial domain pixels after weighted aggregation of its corresponding pixels, to obtain a preliminary denoised complex interference image. The target complex interferometric image is obtained by performing secondary filtering on the initially denoised complex interferometric image.

2. The SAR interferometric phase denoising method based on BM3D according to claim 1, characterized in that, The target complex interferometric image is obtained by performing a second filtering on the initially denoised complex interferometric image, including: In the pixel-level two-dimensional thread grid, each thread simultaneously calculates the similarity between the candidate block and the reference block centered on its corresponding pixel in the search window of the candidate block in the complex interferometric image after preliminary denoising. A new set of similar block position indices is obtained by filtering, and maximum likelihood estimation and decorrelation transformation are performed on the corresponding pixels to obtain the real and imaginary components of the corresponding pixels after preliminary denoising. A second-layer two-dimensional threaded grid with dimensions of image width and search window width is constructed. Each thread in each thread block of the second-layer two-dimensional threaded grid simultaneously performs two-dimensional discrete cosine transform on the real and imaginary components of each pixel in the corresponding image block of each thread block, as well as the real and imaginary components after preliminary denoising, to obtain the transform coefficient set of each image block. In the 3D thread mesh, each thread stacks the transformation coefficients of each pixel according to the new set of similar block position indices based on the transformation coefficient group of the image block, to obtain a new real part 3D block group and an imaginary part 3D block group, and performs Wiener filtering on the corresponding transformation coefficients in the new real part 3D block group and imaginary part 3D block group. In the second-layer thread grid, each thread block simultaneously performs an inverse transform on its corresponding image block after Wiener filtering to obtain the spatial domain group of its corresponding image block after Wiener filtering; in the first-layer two-dimensional thread grid, each thread simultaneously performs weighted aggregation on the spatial domain pixels of its corresponding pixels after Wiener filtering and returns to the complex interference image after preliminary denoising. In a pixel-level two-dimensional threaded grid, each thread simultaneously performs an inverse decorrelation transformation on the spatial domain pixels of its corresponding pixels in the pre-denoised complex interferometric image after weighted aggregation, to obtain the target complex interferometric image.

3. The SAR interferometric phase denoising method based on BM3D according to claim 2, characterized in that, Wiener shrinkage coefficient of the real part of a pixel during Wiener filtering The calculation formula is: ， in, Represents the transformation result matrix; Represents the estimated coherence coefficient of a pixel; Wiener shrinkage coefficient of the imaginary part of a pixel during Wiener filtering The calculation formula is: 。 4. The SAR interferometric phase denoising method based on BM3D according to claim 3, characterized in that, The real part weight of the pixel when performing weighted aggregation of the spatial domain pixels corresponding to the Wiener filter. The calculation formula is: ， The imaginary part weight of the pixel when performing weighted aggregation of the spatial domain pixels corresponding to the Wiener filter. The calculation formula is: 。 5. The SAR interferometric phase denoising method based on BM3D according to claim 1, characterized in that, In a pixel-level two-dimensional threaded grid, each thread simultaneously traverses the candidate blocks within its corresponding pixel's search window, calculates the similarity between the candidate block and its corresponding pixel-centered reference block, and filters to obtain a set of similar block position indices, including: In a pixel-level two-dimensional thread grid, each thread calculates the cosine distance between each pixel of the reference block where its corresponding pixel is located and each pixel of each candidate block in its search window. The cosine distances between each pixel are summed to obtain the similarity between the reference block and each candidate block in its search window. In a pixel-level two-dimensional thread grid, each thread compares the similarity between the reference block containing its corresponding pixel and each candidate block in its search window with a preset threshold, and takes the candidate block with a similarity greater than the preset threshold as the potential similar block of the reference block. In a pixel-level two-dimensional threaded grid, each thread sorts the potential similar blocks of the reference block containing its corresponding pixel in descending order of similarity, and obtains the set of similar block position indices of the reference block based on the position indices of the first preset number of potential similar blocks.

6. The SAR interferometric phase denoising method based on BM3D according to claim 1, characterized in that, In a pixel-level two-dimensional threaded mesh, each thread performs maximum likelihood estimation and decorrelation transformation on its corresponding pixel to obtain the real and imaginary components of that pixel, including: In a pixel-level two-dimensional threaded mesh, each thread calculates the mean intensity of all pixels in the similar blocks of the reference block where its corresponding pixel is located, and obtains the amplitude squared estimate of the similar blocks of the reference block. Based on the amplitude squared estimate of all similar blocks of the reference block, the complex coherence coefficient of the set of similar blocks of the reference block is calculated, thereby obtaining the coherence coefficient estimate and phase estimate of the corresponding pixel. In a pixel-level two-dimensional threaded mesh, each thread uses the decorrelation transformation model and the phase estimate of its corresponding pixel to perform orthogonal component decomposition on the complex interference signal of its corresponding pixel, thereby obtaining the real and imaginary components of the pixel.

7. The SAR interferometric phase denoising method based on BM3D according to claim 1, characterized in that, In the 3D threaded mesh, each thread performs hard thresholding filtering on the transformation coefficients of its corresponding pixels in the transformed 3D block group, including: In the 3D thread mesh, each thread calculates the real part filtering threshold and the imaginary part filtering threshold of its corresponding pixel based on the estimated coherence coefficient of that pixel. In the 3D thread mesh, each thread compares the real part 3D transformation coefficient of the corresponding pixel in the transformed real part 3D block group with the real part filtering threshold of the pixel. If the real part 3D transformation coefficient is greater than the real part filtering threshold of the pixel, the real part 3D transformation coefficient of the pixel remains unchanged. If the real part 3D transformation coefficient is less than or equal to the real part filtering threshold of the pixel, the real part 3D transformation coefficient of the pixel is set to zero. In the 3D thread mesh, each thread compares the virtual part 3D transformation coefficient of the corresponding pixel in the transformed virtual part 3D block group with the virtual part filtering threshold of the pixel. If the virtual part 3D transformation coefficient is greater than the virtual part filtering threshold of the pixel, the virtual part 3D transformation coefficient of the pixel remains unchanged. If the virtual part 3D transformation coefficient is less than or equal to the virtual part filtering threshold of the pixel, the virtual part 3D transformation coefficient of the pixel is set to zero. In the 3D thread mesh, each thread calculates the real part weight of the pixel based on the number of pixels in the similar block to which its corresponding pixel has a non-zero real part 3D transform coefficient after hard thresholding and the estimated coherence coefficient of that pixel. In the 3D thread mesh, each thread calculates the imaginary part weight of a pixel based on the number of pixels in the similar block to which its corresponding pixel has a non-zero imaginary part 3D transform coefficient after hard thresholding and the estimated coherence coefficient of that pixel.

8. The SAR interferometric phase denoising method based on BM3D according to claim 7, characterized in that, In the first layer of the two-dimensional threaded grid, each thread simultaneously performs weighted aggregation of the spatial domain pixels of its corresponding pixels and returns the result to the complex interference image, including: In the first layer of the two-dimensional threaded mesh, each thread calculates the product of the real part spatial domain and the real part weight of its corresponding pixel to obtain the weighted real part pixel value of that pixel. In the first layer of the two-dimensional threaded grid, each thread calculates the product of the imaginary part spatial domain and the imaginary part weight of its corresponding pixel to obtain the weighted imaginary part pixel value of that pixel. In the first layer of the two-dimensional threaded grid, each thread uses atomic addition to accumulate the weighted real part pixel value and the weighted imaginary part pixel value of its corresponding pixel in the three-dimensional block group to the corresponding position in the complex interference image, thereby obtaining the real part pixel value and the imaginary part pixel value of that pixel.

9. The SAR interferometric phase denoising method based on BM3D according to claim 7, characterized in that, Real part filtering threshold of a pixel The calculation formula is: ， in, It is a constant; Represents the estimated coherence coefficient of a pixel; Imaginary part filtering threshold of a pixel The calculation formula is: ， Real part weight of a pixel The calculation formula is: ， in, This represents the number of pixels in the similar block to which the pixel belongs that have non-zero real-part 3D transform coefficients after hard thresholding filtering; Imaginary part weight of a pixel The calculation formula is: ， in, This represents the number of pixels in the similar block to which the pixel belongs whose imaginary 3D transformation coefficients are non-zero after hard thresholding.

10. A SAR interferometric phase denoising device based on BM3D, characterized in that, include: The thread mesh construction module is used to construct pixel-level two-dimensional thread meshes with dimensions of image width and image height, block-level two-dimensional thread meshes with dimensions of image width and search window width, block-level transformation two-dimensional thread meshes with dimensions of image width and number of pixels in image blocks, first-level layered two-dimensional thread meshes with dimensions of image width and number of similar blocks, and three-dimensional thread meshes with dimensions of image width, number of similar blocks and number of pixels in similar blocks. The similar block filtering and related transformation module is used to utilize each thread in the pixel-level two-dimensional thread grid to simultaneously traverse the candidate blocks in the search window where the corresponding pixel is located, calculate the similarity between the candidate block and the reference block centered on the corresponding pixel, filter to obtain the set of similar block position indices, and perform maximum likelihood estimation and decorrelation transformation on the corresponding pixels to obtain the real and imaginary components of the corresponding pixels. The transform coefficient acquisition module is used to simultaneously perform two-dimensional bioorthogonal wavelet transform on the corresponding image blocks by each thread in the block-level two-dimensional thread grid, and obtain the transform coefficient group of the corresponding image blocks. The collaborative filtering module is used to stack the transformation coefficients of each pixel according to the set of similar block position indices in the three-dimensional thread mesh to obtain the real part three-dimensional block group and the imaginary part three-dimensional block group; the block-level transformation module is used to simultaneously perform one-dimensional fast Hadamard transformation on the third dimension of the three-dimensional block group where the corresponding pixel is located in the two-dimensional thread mesh; and the three-dimensional thread mesh is used to perform hard thresholding filtering on the transformation coefficients of the corresponding pixels in the transformed three-dimensional block group. The inverse transformation and weighted aggregation module is used to simultaneously perform inverse transformation on the corresponding image blocks after hard thresholding by each thread in the block-level two-dimensional thread grid to obtain the spatial domain group of the corresponding image blocks; each thread in the first layer two-dimensional thread grid simultaneously performs weighted aggregation on the spatial domain pixels of its corresponding pixels and returns it to the complex interference image. The preliminary denoised image acquisition module is used to perform inverse decorrelation transformation on the spatial domain pixels after weighted aggregation of corresponding pixels by each thread in the pixel-level two-dimensional thread grid, so as to obtain the complex interference image after preliminary denoising. The secondary filtering module is used to perform secondary filtering on the initially denoised complex interferometric image to obtain the target complex interferometric image.

Images