Remote sensing image self-adaptive processing method and system combining optimization algorithm and probability model

Abstract

The application provides a kind of remote sensing image adaptive processing method and system combining optimization algorithm and probability model, it is related to remote sensing image processing technical field, first, the original projection data set of original remote sensing image is carried out projection domain data reconstruction processing based on Landweber iteration framework, generates initial density distribution image;Then initial density distribution image is carried out adaptive segmentation processing based on Gaussian mixture model and neighborhood correlation constraint, obtains superpixel segmentation result image;Then superpixel segmentation result image is carried out Gaussian mixture model parameter estimation processing based on expectation maximization algorithm, constructs multidimensional Gaussian probability density function set;Again, the posterior probability distribution matrix of each pixel point belonging to each Gaussian component is calculated according to the set;Finally, edge preserving iterative optimization processing is carried out according to posterior probability distribution matrix and image gradient information, generates final processing result image.The application can effectively improve remote sensing image quality, realize edge enhancement and noise suppression.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image processing technology, and more specifically, to a method and system for adaptive processing of remote sensing images that combines optimization algorithms and probabilistic models. Background Technology

[0002] Traditional remote sensing image processing methods often have many limitations. For example, in the reconstruction of projection domain data, some conventional methods may not be able to fully consider the characteristics of the data, resulting in poor quality of the initial density distribution image, with problems such as blurring and noise interference, which affects the accuracy of subsequent processing.

[0003] In the image segmentation process, existing segmentation methods may struggle to effectively handle complex scenes and diverse land cover types in remote sensing images. The segmentation results often suffer from unclear boundaries and inaccurate region division, failing to accurately reflect the true structure of the image.

[0004] For Gaussian mixture model parameter estimation, traditional methods may only be based on simple statistical information and cannot fully explore the deep features of the image. The constructed set of multidimensional Gaussian probability density functions is difficult to accurately characterize the overall pixel distribution characteristics of the image.

[0005] Furthermore, in terms of edge preservation and noise suppression, existing technologies struggle to accurately preserve image edge information while effectively removing noise, resulting in the final processed image quality failing to meet the needs of practical applications. Summary of the Invention

[0006] In view of the aforementioned problems, and in conjunction with the first aspect of the present invention, embodiments of the present invention provide a remote sensing image adaptive processing method combining optimization algorithms and probabilistic models, the method comprising: The original projection data set corresponding to the acquired original remote sensing image is subjected to projection domain data reconstruction processing based on the Landweber iterative framework. The original projection data set is successively approximated and solved according to the preset relaxation factor parameters and the preset convergence judgment threshold to generate an initial density distribution image. An adaptive segmentation process based on Gaussian mixture model and neighborhood correlation constraint is performed on the initial density distribution image. Markov random field theory is introduced to construct a neighborhood dependency weight for each pixel in the initial density distribution image to measure its consistency with the neighboring pixel labels. In the iterative segmentation framework based on Gaussian mixture model and Markov random field coupling, the Gibbs energy function is constructed using the neighborhood dependency weight to perform pixel category classification processing on the initial density distribution image and generate a superpixel segmentation result image containing smooth region boundaries. The superpixel segmentation result image is subjected to Gaussian mixture model parameter estimation processing based on the expectation-maximization algorithm. The expectation parameter, variance parameter and prior probability parameter of each Gaussian component in the current Gaussian mixture model are calculated using the pixel statistics of each superpixel region in the superpixel segmentation result image. A set of multidimensional Gaussian probability density functions is constructed based on the expectation parameter, the variance parameter and the prior probability parameter to characterize the overall pixel distribution characteristics of the superpixel segmentation result image. Based on the set of multidimensional Gaussian probability density functions, Bayesian posterior probability calculation is performed on each pixel in the superpixel segmentation result image to generate a posterior probability distribution matrix for each pixel belonging to each Gaussian component. Based on the posterior probability distribution matrix and the superpixel segmentation result image, perform edge-preserving iterative optimization processing based on image gradient information. Construct a comprehensive objective function that includes a data fidelity term, a neighborhood space constraint term, and an image gradient regularization term. Determine the optimization step size parameter that minimizes the comprehensive objective function through a linear search method. Iteratively update the pixel values ​​of the superpixel segmentation result image according to the optimization step size parameter to generate the final processed result image after edge enhancement and noise suppression.

[0007] Furthermore, embodiments of the present invention also provide a remote sensing image adaptive processing system that combines optimization algorithms and probabilistic models, comprising: A processor; a machine-readable storage medium for storing machine-executable instructions of the processor; wherein the processor is configured to execute the aforementioned remote sensing image adaptive processing method combining optimization algorithms and probabilistic models by executing the machine-executable instructions.

[0008] Based on the above, the projection domain data reconstruction processing based on the Landweber iterative framework, using preset relaxation factor parameters and convergence thresholds for successive approximation solutions, can generate a relatively accurate initial density distribution image, effectively overcoming the blurring and noise problems that may occur in traditional reconstruction methods. Secondly, Markov random field theory is introduced to construct neighborhood-dependent weights, and an iterative segmentation framework coupled with Gaussian mixture models and Markov random fields is used for adaptive segmentation processing, generating superpixel segmentation result images containing smooth region boundaries. This greatly improves the accuracy of image segmentation, enabling more precise classification of pixel categories and adapting to complex terrain scenes in remote sensing images. Furthermore, the Gaussian mixture model parameters are estimated using the expectation-maximization algorithm to construct a multidimensional Gaussian probability density function set, accurately representing the overall pixel distribution characteristics of the superpixel segmentation result image. Finally, a comprehensive objective function is constructed based on the posterior probability distribution matrix and image gradient information, and edge-preserving iterative optimization processing is performed to generate a final processed image with edge enhancement and noise suppression. This effectively removes noise while accurately preserving image edge information, significantly improving the quality of remote sensing images and meeting the practical application requirements for high-quality remote sensing images. Attached Figure Description

[0009] Figure 1 This is a schematic diagram of the execution flow of the remote sensing image adaptive processing method that combines optimization algorithms and probability models provided in an embodiment of the present invention.

[0010] Figure 2 This is a schematic diagram of exemplary hardware and software components of a remote sensing image adaptive processing system that combines optimization algorithms and probability models, as provided in an embodiment of the present invention. Detailed Implementation

[0011] Figure 1 This is a flowchart illustrating an embodiment of the adaptive remote sensing image processing method combining optimization algorithms and probability models, which will be described in detail below.

[0012] Step S110: Perform projection domain data reconstruction processing based on the Landweber iterative framework on the original projection data set corresponding to the acquired original remote sensing image. According to the preset relaxation factor parameters and the preset convergence judgment threshold, the original projection data set is successively approximated to generate an initial density distribution image.

[0013] In this embodiment, taking the reconstruction of remote sensing images of urban areas as an example, the original projection data set corresponding to the original remote sensing image comes from the imaging system of a high-resolution remote sensing satellite. This imaging system collects electromagnetic radiation signals of ground objects at different orbital positions and attitudes through a detector array, forming original projection data containing the correlation between spatial location and radiation intensity. The Landweber iterative framework, as a classic iterative reconstruction algorithm, approximates the true density distribution by progressively correcting the current solution vector, and is suitable for handling noise and incomplete projection problems in remote sensing images. In the entire processing flow, the relaxation factor parameter determines the correction step size of each iteration, while the convergence criterion threshold is used to control the timing of iteration termination; both together affect the quality of the reconstructed image and the computational efficiency.

[0014] Step S111: Obtain the original projection data set of the original remote sensing image collected by the detector array during the imaging process. The original projection data set contains a sequence of ray attenuation intensity values ​​recorded by multiple detector channels at different rotation angles.

[0015] In urban remote sensing imaging scenarios, detector arrays typically consist of K detector units, each corresponding to an independent detector channel. As the satellite orbits the Earth, the imaging system rotates at a certain angular velocity, forming L different rotation angles. Each detector channel continuously collects ray attenuation intensity values ​​at T sampling points at each rotation angle. Therefore, the original projection dataset is a three-dimensional array with dimensions [K, L, T], where K represents the number of detector channels, L represents the number of rotation angles, and T represents the number of sampling points at each angle. For example, when K=128, L=360, and T=512, the original projection dataset contains 128×360×512 ray attenuation intensity values, each value representing dimensionless relative radiance (the dimension effect has been eliminated through system calibration). The data source is Level 1B data transmitted from the satellite after preliminary radiometric correction, containing no privacy-sensitive information and conforming to industry standards for remote sensing data acquisition.

[0016] Step S112: Construct a system response matrix based on the geometric parameters of the imaging system to describe the linear mapping relationship between the original projection data set and the image to be reconstructed. The row index of the system response matrix corresponds to the ray path of the detector channel, and the column index of the system response matrix corresponds to the pixel unit position in the image space to be reconstructed.

[0017] The geometric parameters of the imaging system include the physical dimensions of the detector array, center distance, rotation axis position, and focal length. The image to be reconstructed is divided into M×N pixel units, each with a spatial dimension of Δx×Δy (in meters). The system response matrix H has dimensions [K×L×T, M×N], where the total number of rows equals the total number of samples in the original projection data set (K×L×T), and the total number of columns equals the total number of pixels in the image to be reconstructed (M×N). The matrix element H(i,j) represents the contribution weight of the j-th pixel unit to the i-th ray path, and its value is calculated using a ray tracing algorithm: when the i-th ray passes through the j-th pixel unit, H(i,j) equals the ratio of the ray's path length within the pixel to the pixel's physical size; otherwise, H(i,j) is 0. For example, for a pixel cell centered at image coordinates (x0, y0), if a ray enters from the top left corner of the pixel (x0-Δx / 2, y0+Δy / 2) and exits from the bottom right corner (x0+Δx / 2, y0-Δy / 2), the path length d is obtained by calculating the distance between the two points. Then H(i, j) = d / (Δx√2) (when Δx = Δy), ensuring that the weight values ​​are normalized between 0 and 1.

[0018] Step S113: Set the initial iterative solution vector of the Landweber iterative framework. Each element in the initial iterative solution vector corresponds to the initial density estimate of each pixel unit in the image to be reconstructed. Assign all elements in the initial iterative solution vector to a preset zero initial state and perform normalization preprocessing on the initial iterative solution vector.

[0019] The initial iterative solution vector x0 is an M×N dimensional column vector, where each element x0(j) corresponds to the initial density estimate of the j-th pixel unit in the image to be reconstructed. To avoid numerical fluctuations in the early stages of iteration, all elements of x0 are initialized to 0. Normalization preprocessing is implemented as follows: the maximum value Imax and the minimum value Imin of all ray attenuation intensity values ​​in the original projection data set are calculated, and the numerical range of the initial solution vector x0 is mapped to the interval [Imin / Imax, 1]. The specific processing formula is x0_norm(j) = x0(j) / Imax + Imin / Imax, ensuring that the preprocessed data has the same dimensions and magnitude as the original projection data.

[0020] Step S114: Set the relaxation factor parameter and convergence threshold of the Landweber iteration framework. The relaxation factor parameter is used to control the magnitude of the correction of the current solution vector along the negative gradient direction in each iteration. The convergence threshold is used to determine whether the iteration process has reached the termination condition.

[0021] The relaxation factor parameter λ ranges from (0, 2 / λmax), where λmax is the largest eigenvalue of the system response matrix H. After calculating the largest eigenvalue λmax of H^TH through eigenvalue decomposition, λ is set to 0.5 / λmax to ensure the convergence of the iteration process. The convergence threshold ε is set to 0.1% of the largest element value in the original projected dataset, i.e., ε = 0.001 × Imax. The iteration is considered convergent when the norm of the iterative residual is less than ε. For example, if Imax = 1000, then ε = 1.0, and the iteration stops when the L2 norm of the residual vector is less than 1.0.

[0022] Step S115: Calculate the residual vector of the current iteration based on the system response matrix, the original projection data set, the relaxation factor parameter, and the current iteration solution vector. The residual vector is obtained by multiplying the system response matrix with the current iteration solution vector to obtain the current estimated projection data, and then subtracting the current estimated projection data from the corresponding projection value in the original projection data set element by element.

[0023] The current iterative solution vector is xk (k is the iteration number; initially, when k=0, xk=x0), and the original projected data set is vectorized into b (dimension K×L×T). The current estimated projected data is calculated using matrix multiplication: yk=H×xk, where yk is a K×L×T column vector. The residual vector rk=b-yk, where each element rk(i)=b(i)-yk(i) represents the difference between the estimated value and the true value at the i-th sampling point. For example, when k=0, y0=H×x0=0 (because x0 is a zero vector), and at this time r0=b-0=b, meaning the initial residual equals the original projected data.

[0024] Step S116: Before executing the Landweber iterative framework, the original projected data set and the initial iterative solution vector are normalized and preprocessed. During the iteration process, the normalized current iterative solution vector is updated according to the transpose of the system response matrix, the relaxation factor parameter, and the residual vector to generate a new normalized iterative solution vector. The new normalized iterative solution vector is obtained by adding the product of the relaxation factor parameter and the transpose of the system response matrix multiplied by the residual vector to the normalized current iterative solution vector.

[0025] The normalization preprocessing was completed in step S113, and the normalized current iterative solution vector is denoted as xk_norm. The transpose of the system response matrix is ​​H^T (dimensions [M×N, K×L×T]). The iterative update formula is: xk+1_norm = xk_norm + λ×H^T×rk. Here, H^T×rk represents the backprojection of residual information from the projection domain to the image domain, and λ controls the intensity of the backprojection. For example, when k=0, x1_norm = x0_norm + λ×H^T×r0 = 0 + λ×H^T×b, yielding the normalized solution vector after the first iteration.

[0026] Step S117: Calculate the norm of the new residual vector corresponding to the new normalized iterative solution vector, and compare the norm of the new residual vector with the convergence determination threshold. If the norm of the new residual vector is less than the convergence determination threshold, it is determined that the convergence condition has been met and the iteration process of the Landweber iterative framework is stopped. If the norm of the new residual vector is greater than or equal to the convergence determination threshold, the new normalized iterative solution vector is used as the current iterative solution vector to continue the next iteration update.

[0027] The new residual vector rk+1 = bH × xk+1 has its L2 norm calculated as ||rk+1||2 = sqrt(sum(rk+1(i)^2)). Compare ||rk+1||2 with ε. If ||rk+1||2 < ε, stop the iteration; otherwise, let xk = xk+1, k = k+1, and repeat steps S115 to S117. For example, when the iteration reaches the 20th iteration, if ||r20||2 = 0.8 < ε = 1.0, then convergence is determined, and the iteration stops.

[0028] Step S118: Record the final normalized iterative solution vector when the convergence condition is met, perform inverse normalization on the final normalized iterative solution vector, and reorganize the inverse normalized final iterative solution vector according to the two-dimensional arrangement order of pixel units in the image space to generate the initial density distribution image with spatial dimension information.

[0029] The final normalized iterative solution vector is x_final_norm, and the inverse normalization formula is x_final(j) = (x_final_norm(j) - Imin / Imax) × Imax, restoring the original density dimensions of the pixel unit. x_final (an M×N dimensional column vector) is then rearranged into an M-row, N-column two-dimensional matrix in row-major order. Matrix element (i, j) corresponds to the density value of the pixel in the i-th row and j-th column of the image space, thus generating the initial density distribution image. This initial density distribution image has a spatial resolution of M×N, and the density value range of each pixel is consistent with the physical meaning of the original projection data (e.g., units of radiation intensity).

[0030] Step S119: Perform condition number pre-analysis on the system response matrix and calculate the ratio of the maximum singular value to the minimum singular value of the system response matrix as the condition number value.

[0031] The singular value decomposition of the system response matrix H is H = UΣV^T, where Σ is a diagonal matrix and the diagonal elements are singular values ​​σ1 ≥ σ2 ≥ ... ≥ σmin (K × L × T, M × N) ≥ 0. The condition number κ(H) = σ1 / σmin is used to measure the ill-conditioned nature of the matrix. After calculating σ1 and σmin using the singular value decomposition algorithm, κ(H) can be obtained. For example, when σ1 = 1000 and σmin = 0.1, κ(H) = 10000, indicating that the matrix may have ill-conditioned properties.

[0032] Step S1110: Compare the value of the condition number with a preset ill-condition threshold. If the value of the condition number is greater than the preset ill-condition threshold, then the system response matrix is ​​determined to exhibit ill-condition characteristics.

[0033] The preset ill-conditioned threshold is usually set to 1000 (determined based on empirical values ​​in the field of remote sensing image reconstruction). When κ(H) > 1000, the matrix is ​​deemed ill-conditioned, and directly using Landweber iteration in this case may lead to unstable solutions or slow convergence. For example, if κ(H) = 1500 > 1000, then subsequent diagonal loading preprocessing is performed.

[0034] Step S1111: Based on the determination result of the ill-conditioned characteristics, perform a diagonal loading preprocessing operation on the system response matrix to generate a corrected system response matrix after diagonal loading correction. The corrected system response matrix is ​​obtained by adding a preset diagonal loading coefficient to the system response matrix and multiplying it by the identity matrix.

[0035] The diagonal loading coefficient μ is dynamically set based on the condition number, calculated using the formula μ = σmin × (κ(H) / 1000 - 1), ensuring that the condition number of the matrix drops below 1000 after loading. The corrected system response matrix H' = H + μI, where I is the identity matrix. For example, when κ(H) = 1500 and σmin = 0.1, μ = 0.1 × (1500 / 1000 - 1) = 0.05, at which point the condition number κ(H') of H' is 1000.

[0036] Step S1112: Use the corrected system response matrix as the system response matrix used in subsequent calculations of the current estimated projection data.

[0037] In subsequent iterations, H' is used instead of H for matrix multiplication, i.e., yk=H'×xk, to improve the ill-conditioned properties of the matrix and enhance the stability of the iteration.

[0038] Step S1113: Calculate the current estimated projection data after diagonal loading correction based on the corrected system response matrix and the current iterative solution vector. The current estimated projection data after diagonal loading correction is obtained by multiplying the corrected system response matrix with the current iterative solution vector.

[0039] Similar to step S115, but using H' for calculation: yk'=H'×xk, where yk' is the corrected estimated projection data.

[0040] Step S1114: Perform element-wise subtraction between the current estimated projection data after diagonal loading correction and the corresponding projection value in the original projection data set to obtain the residual vector after diagonal loading correction.

[0041] The residual vector rk'=b-yk' is corrected for subsequent iterative updates.

[0042] Step S1115: Update the current iterative solution vector according to the transpose of the corrected system response matrix, the relaxation factor parameter, and the diagonally loaded corrected residual vector to generate a new iterative solution vector.

[0043] The updated formula is xk+1=xk+λ×(H')^T×rk', where (H')^T is the transpose of H'.

[0044] Step S1116: In each iteration, record the total variation value of the image corresponding to the new iterative solution vector after diagonal loading correction. The total variation value of the image is obtained by calculating the sum of the absolute differences of the new iterative solution vector in the horizontal and vertical directions.

[0045] The new iterative solution vector xk+1 is reorganized into a two-dimensional image matrix G (M×N). The elements of the horizontal difference matrix D_x are D_x(i,j) = |G(i,j+1) - G(i,j)| (j=1..N-1), and the elements of the vertical difference matrix D_y are D_y(i,j) = |G(i+1,j) - G(i,j)| (i=1..M-1). The total variational value TV = sum(D_x) + sum(D_y) is used to measure the smoothness of the image.

[0046] Step S1117: Compare the total variation value of the image with the preset total variation constraint threshold. If the total variation value of the image is greater than the preset total variation constraint threshold, perform soft threshold shrinkage processing on the new iterative solution vector to generate an iterative solution vector with total variation constraint, and use the iterative solution vector with total variation constraint as the current iterative solution vector for the next iteration.

[0047] The preset total variation constraint threshold TV_threshold is set according to the image complexity, for example, 50% of the initial total variation (when k=0). When TV>TV_threshold, soft thresholding is performed on the image matrix G: G'(i,j)=sign(G(i,j))×max(|G(i,j)|-τ,0), where τ is the threshold parameter, τ=TV_threshold / (M×N). The shrunk matrix G' is revectorized into xk+1', which serves as the current solution vector for the next iteration to suppress noise and artifacts in the image.

[0048] Step S120: Perform adaptive segmentation processing based on Gaussian mixture model and neighborhood correlation constraint on the initial density distribution image. Introduce Markov random field theory to construct a neighborhood dependency weight for each pixel in the initial density distribution image to measure its consistency with the neighboring pixel labels. In the iterative segmentation framework based on Gaussian mixture model and Markov random field coupling, construct the Gibbs energy function using the neighborhood dependency weight to perform pixel category classification processing on the initial density distribution image and generate a superpixel segmentation result image containing smooth region boundaries.

[0049] In urban remote sensing image scenarios, the initial density distribution image contains different land cover types such as buildings, roads, vegetation, and water bodies, each with different density distribution characteristics. Gaussian Mixture Models (GMMs) can fit complex pixel grayscale distributions using multiple Gaussian components, while Markov Random Fields (MRFs) can constrain label assignment using spatial correlations between pixels. The coupling of these two methods enables adaptive segmentation that considers both spectral features and spatial structure. Neighborhood dependency weights quantify the impact of label consistency between the center pixel and its neighbors on the segmentation result. The Gibbs energy function integrates the data term (fit to the GMM) and the spatial term (neighborhood label consistency) as the segmentation optimization objective. By minimizing the energy function, the optimal partitioning of pixel categories is achieved, ultimately generating superpixel regions with smooth boundaries.

[0050] Step S121: Read all pixels of the initial density distribution image and their corresponding gray values, and convert the initial density distribution image into a two-dimensional gray matrix. The row index of the two-dimensional gray matrix corresponds to the vertical coordinate of the pixel in the image space, and the column index of the two-dimensional gray matrix corresponds to the horizontal coordinate of the pixel in the image space.

[0051] The initial density distribution image is an M×N two-dimensional image, where each pixel (i, j) (i=1..M, j=1..N) has a corresponding gray value g(i, j) (normalized to the [0, 1] interval). The image data is stored as a two-dimensional gray matrix G, where G(i, j) = g(i, j), the row index i corresponds to the vertical coordinate of the image (e.g., the latitude direction of a remote sensing image), and the column index j corresponds to the horizontal coordinate (longitude direction). For example, for a 512×512 image, G is a 512-row, 512-column matrix, where each element represents the density normalized value of the pixel at the corresponding location.

[0052] Step S122: Perform initial cluster center point generation processing on the two-dimensional grayscale matrix according to the preset superpixel quantity parameter. Generate multiple initial cluster center point coordinates with the same value as the superpixel quantity parameter within the range of pixel coordinate values ​​of the two-dimensional grayscale matrix using a random number generation function. Each initial cluster center point coordinate corresponds to the initial center position of a superpixel region to be generated.

[0053] The preset superpixel count parameter S is set according to the image resolution and terrain complexity. For example, for a 512×512 image, S=200. A random number generation function generates S unique initial cluster center point coordinates (c1x, c1y), (c2x, c2y), ..., (csx, csy) within the coordinate range [1, M]×[1, N]. To avoid initial center point clustering, a uniform grid sampling combined with random perturbation can be used: the image is uniformly divided into a √S×√S grid, and a random perturbation within the range [-d, d] (d is half the grid side length) is added to the center of each grid to ensure uniform initial center point distribution.

[0054] Step S123: Using the coordinates of the initial cluster center points as a reference, perform a comprehensive distance metric calculation on each pixel in the two-dimensional grayscale matrix, calculating the normalized spatial Euclidean distance and the normalized grayscale difference distance between each pixel and each initial cluster center point; the normalized spatial Euclidean distance is calculated based on the difference in row coordinates and column coordinates of the pixel in spatial position; the normalized grayscale difference distance is calculated based on the difference in grayscale values ​​of the pixel; the normalized spatial Euclidean distance and the normalized grayscale difference distance are linearly combined according to a preset weight to obtain the comprehensive distance.

[0055] For pixels (i, j) and cluster centers (ckx, cky): the spatial Euclidean distance d_s = sqrt((i-ckx)^2 + (j-cky)^2), and the normalized spatial distance d_s_norm = d_s / D_max, where D_max is the length of the image diagonal sqrt(M^2 + N^2); the grayscale difference distance d_g = |G(i, j) - G(ckx, cky)|, and the normalized grayscale distance d_g_norm = d_g / G_range (G_range = 1, since the grayscale has been normalized). Preset weights α (spatial weight) and 1-α (grayscale weight), α is usually taken as 0.5, and the combined distance d = α × d_s_norm + (1-α) × d_g_norm. For example, when α = 0.5, d_s_norm = 0.3, and d_g_norm = 0.4, d = 0.35.

[0056] Step S124: Based on the comprehensive distance calculation results between each pixel and each initial cluster center point, assign each pixel to the superpixel region category corresponding to the initial cluster center point with the smallest Euclidean distance value, forming the initial superpixel region set.

[0057] For each pixel (i, j), calculate its combined distances d1, d2, ..., dS with all S cluster centers. Find the cluster center index k* = argmin(dk) corresponding to the minimum distance, and assign this pixel to the k*th superpixel region. After all pixels are assigned, S superpixel regions are obtained, each containing several pixels, forming the initial partitioning result.

[0058] Step S125: Perform cluster center point update processing on each superpixel region in the initially divided superpixel region set, calculate the average spatial coordinates and average grayscale value of all pixels in each superpixel region, and use the average spatial coordinates and average grayscale value as the new cluster center point coordinates.

[0059] For the k-th superpixel region, which contains a set of pixels Pk={(i1,j1),(i2,j2),...,(in,jn)}, calculate the average spatial coordinates (ckx_new, cky_new)=(mean(ik),mean(jk)), and the average grayscale value g_avg=mean(G(ik,jk)). Use (ckx_new, cky_new) as the new cluster center coordinates, and g_avg as the representative grayscale value of the region.

[0060] Step S126: Repeat the comprehensive distance metric calculation and combination processing, the pixel point allocation processing, and the cluster center point update processing until the change in the coordinates of the cluster center points in two adjacent iterations is less than the preset clustering convergence threshold, and obtain the initial superpixel segmentation result based on the K-Means clustering algorithm.

[0061] The clustering convergence threshold δ is set to 0.5% of the image spatial resolution. For example, for a 512×512 image, δ = 512 × 0.005 = 2.5 pixels. After each iteration, the sum of the coordinate changes of all cluster center points is calculated as Δ = sum(sqrt((ckx_new - ckx_old)^2 + (cky_new - cky_old)^2)). Iteration stops when Δ < δ. Convergence is usually achieved after 5-10 iterations, yielding the initial superpixel segmentation result. Each superpixel region contains consecutive pixels, and the region boundaries may have jagged edges.

[0062] Step S127: Based on the superpixel region category to which each pixel belongs in the initial superpixel segmentation result, construct an initial label vector for each pixel. The dimension of the initial label vector is equal to the total number of superpixel regions. Only the element corresponding to the superpixel region category to which the pixel belongs has a value of 1 in the initial label vector, and the remaining elements have a value of 0.

[0063] The initial label vector for each pixel (i, j) is an S-dimensional one-hot vector L(i, j) = [l1, l2, ..., lS], where lk = 1 (if the pixel belongs to the k-th superpixel region), otherwise lk = 0. For example, the label vector of a pixel belonging to the 3rd superpixel region is [0, 0, 1, 0, ..., 0].

[0064] Step S128: Construct a Markov random field neighborhood system based on the spatial adjacency relationship between each pixel and its surrounding neighboring pixels in the initial density distribution image. The Markov random field neighborhood system defines a set of neighboring pixels for each pixel. The set of neighboring pixels includes the pixels directly adjacent to the pixel in spatial position above, below, to the left, and to the right, as well as four adjacent pixels in the diagonal direction.

[0065] An 8-neighborhood system is used to define the neighborhood set N(i,j) = {(i-1,j-1), (i-1,j), (i-1,j+1), (i,j-1), (i,j+1), (i+1,j-1), (i+1,j), (i+1,j+1)} for each pixel (i,j), ensuring that the coordinates of all neighboring pixels are within the image range (i.e., i-1≥1, i+1≤M, j-1≥1, j+1≤N). For example, for the pixel (1,1) at the image edge, its neighborhood set only contains three valid pixels: (1,2), (2,1), and (2,2).

[0066] Step S129: Perform weight allocation processing on the set of neighboring pixels of each pixel in the Markov random field neighborhood system. Normalize the spatial distance and gray value difference between all pixels in the initial density distribution image. Calculate the corresponding spatial distance weight coefficient based on the normalized spatial distance between each neighboring pixel and the center pixel. The spatial distance weight coefficient is inversely proportional to the square of the normalized spatial distance. Calculate the corresponding gray value difference weight coefficient based on the normalized gray value difference between each neighboring pixel and the center pixel. The gray value difference weight coefficient is obtained by performing a negative exponential transformation on the absolute value of the normalized gray value difference. Perform a dot product operation between the spatial distance weight coefficient and the gray value difference weight coefficient to obtain the dimensionless comprehensive influence weight of each neighboring pixel relative to the center pixel.

[0067] For the center pixel (i, j) and its neighboring pixels (i', j'): Spatial distance d = sqrt((i-i')^2 + (j-j')^2), normalized spatial distance d_norm = d / d_max (d_max is the maximum distance √2 in the 8-neighborhood); spatial distance weight w_s = 1 / (d_norm^2 + ε), where ε = 0.01 to avoid zero in the denominator. Gray-level difference Δg = |G(i, j) - G(i', j')|, normalized gray-level difference Δg_norm = Δg / G_range (G_range = 1); gray-level difference weight w_g = exp(-β × Δg_norm), where β is a parameter controlling sensitivity (usually taken as 5). The overall influence weight is w = w_s × w_g. For example, when d_norm = √2 / √2 = 1 and Δg_norm = 0.2, w_s = 1 / (1+0.01) = 0.99, w_g = exp(-5×0.2) = exp(-1) = 0.368, and w = 0.99×0.368 ≈ 0.364.

[0068] Step S1210: Based on the comprehensive influence weight, the set of neighboring pixels of each pixel in the Markov random field neighborhood system is weighted to construct a weighted Markov random field neighborhood system.

[0069] The weight w of each neighboring pixel (i', j') is stored in the weight matrix W(i, j, i', j') to form a weighted neighborhood system, so that in the subsequent energy calculation, the influence of different neighboring pixels on the central pixel is different.

[0070] Step S1211: Construct a weighted neighborhood correlation constraint energy function based on the weighted Markov random field neighborhood system and Gibbs distribution. The weighted neighborhood correlation constraint energy function is obtained by multiplying the difference between the label vector of each neighboring pixel and the label vector of the center pixel by the corresponding comprehensive influence weight and then summing them up.

[0071] The weighted neighborhood correlation constraint energy function E_s(L(i,j))=sum_{(i',j')∈N(i,j)}w(i,j,i',j')×||L(i,j)-L(i',j')||_2^2, where ||·||_2^2 is the square of the L2 norm. For a one-hot label vector, ||L(i,j)-L(i',j')||_2^2=2(1-δ(k,k')), where δ(k,k') is the Kronecker function (1 when k=k', 0 otherwise). Therefore, E_s(L(i,j))=sumw×2(1-δ(k,k')). When the neighboring pixel has the same label as the center pixel, this term contributes 0; otherwise, it contributes 2w, ensuring that spatially adjacent pixels with different labels contribute significantly to the energy.

[0072] Step S1212: Generate a corresponding weighted neighborhood dependency weight for each pixel according to the weighted neighborhood correlation constraint energy function. The weighted neighborhood dependency weight is obtained by performing a negative exponential transformation on the weighted neighborhood correlation constraint energy function and then normalizing it.

[0073] The weighted neighborhood dependency weight w_dep(i,j) = exp(-E_s(L(i,j))) / Z, where Z is a normalization constant (the sum of exp(-E_s) for all pixels), ensuring that w_dep(i,j) is in the interval [0,1] and the sum is 1. This weight reflects the degree to which the pixel label is constrained by the neighborhood; the lower the energy (the more consistent the label), the higher the dependency weight.

[0074] Step S1213: Construct a neighborhood correlation constraint energy function based on Markov random field theory and Gibbs distribution. The neighborhood correlation constraint energy function is used to measure the consistency between the current label vector of each pixel and the label vectors of all pixels in its neighborhood pixel set.

[0075] In the unweighted case, the neighborhood correlation constraint energy function is E_s(L(i,j))=sum_{(i',j')∈N(i,j)}||L(i,j)-L(i',j')||_2^2, where all neighboring pixels have equal weights (all 1), which is used for the basic spatial consistency measure.

[0076] Step S1214: Generate a corresponding neighborhood dependency weight for each pixel according to the neighborhood correlation constraint energy function. The neighborhood dependency weight is obtained by performing a negative exponential transformation on the neighborhood correlation constraint energy function and then normalizing it.

[0077] The neighborhood dependency weight w_dep(i,j) = exp(-E_s(L(i,j))) / Z is similar to step S1212, but uses unweighted E_s for calculation, which is suitable for scenarios where spatial correlation requirements are not high. In this embodiment, the weighted neighborhood dependency weight calculated in step S1212 is used.

[0078] Step S1215: Input the gray value of each pixel and the initial label vector into an iterative segmentation framework based on the coupling of Gaussian mixture model and Markov random field, and use the weighted neighborhood dependency weights to construct the Gibbs energy function as the segmentation constraint. Iteratively estimate the parameters of the Gaussian mixture model and the posterior probability of each pixel by alternately executing the expectation step and the maximization step through the expectation-maximization algorithm.

[0079] The Gaussian Mixture Model (GMM) contains K Gaussian components (K=S, the same as the number of superpixels), each with parameters (μk, σk^2, πk), where μk is the mean, σk^2 is the variance, and πk is the prior probability. The Gibbs energy function is E=E_d+γ×E_s, where E_d is the data term (negative log-likelihood of the GMM), E_s is the spatial term (weighted neighborhood correlation constraint energy), and γ is the balance parameter (set to 1.0). The expectation step (E-step) of the Expectation-Maximization (EM) algorithm calculates the posterior probability of a pixel belonging to each Gaussian component; the maximization step (M-step) updates the GMM parameters and simultaneously adjusts the label vector based on the spatial constraints of the MRF, achieving coupled optimization.

[0080] Step S1216: Repeat the expectation step and the maximization step until the parameter change of the Gaussian mixture model is less than the preset model convergence threshold. Based on the posterior probability of each pixel belonging to each Gaussian component obtained in the final iteration, select the category corresponding to the Gaussian component with the largest posterior probability value as the final superpixel region category of the pixel.

[0081] The model convergence threshold is set to the sum of the changes in all Gaussian component parameters being less than 1e-3, i.e., sum(|μk_new-μk_old|+|σk_new^2-σk_old^2|+|πk_new-πk_old|)<1e-3. After iterative convergence, for each pixel, the index of the Gaussian component with the highest posterior probability is selected as its final class label, thus achieving pixel class classification.

[0082] Step S1217: Perform region labeling processing on the two-dimensional grayscale matrix according to the final superpixel region category of each pixel, merge adjacent pixels with the same final superpixel region category into the same superpixel region, and smooth the boundary of the merged superpixel region to generate the superpixel segmentation result image containing the smoothed region boundary.

[0083] Region labeling is achieved through connected component analysis, merging spatially adjacent pixels with the same category label into superpixel regions. Boundary smoothing employs the closing operation (dilation followed by erosion) in morphological operations, with a 3×3 circular kernel as the structuring element, eliminating fine jagged edges on the boundaries and making the superpixel region boundaries smoother. In the final generated superpixel segmentation image, pixels within each superpixel region have similar grayscale characteristics and spatial continuity.

[0084] Step S130: Perform Gaussian mixture model parameter estimation processing based on the expectation-maximization algorithm on the superpixel segmentation result image. Calculate the expectation parameter, variance parameter, and prior probability parameter of each Gaussian component in the current Gaussian mixture model using the pixel statistics of each superpixel region in the superpixel segmentation result image. Construct a set of multidimensional Gaussian probability density functions to characterize the overall pixel distribution characteristics of the superpixel segmentation result image based on the expectation parameter, the variance parameter, and the prior probability parameter.

[0085] In the superpixel segmentation result image, each superpixel region represents a land cover region with similar gray-level characteristics, such as building areas, vegetation areas, etc. The parameter estimation of the Gaussian Mixture Model (GMM) aims to optimize the model parameters using the pixel statistics (mean, variance, etc.) of each superpixel region, enabling the model to accurately describe the overall pixel distribution of the image. The Expectation-Maximization (EM) algorithm achieves maximum likelihood estimation of the parameters by iteratively estimating the probability of a pixel's association with each Gaussian component (E-step) and updating the model parameters (M-step). The set of multidimensional Gaussian probability density functions is the specific mathematical expression of each Gaussian component, used for subsequent pixel category probability calculations. In urban remote sensing images, this step can decompose the image pixel distribution into multiple Gaussian distributions corresponding to different land cover types.

[0086] Step S131: Analyze the superpixel segmentation result image to obtain the total number of superpixel regions contained in the superpixel segmentation result image and the grayscale value sequence of all pixels contained in each superpixel region.

[0087] The superpixel segmentation result image contains S superpixel regions (consistent with the superpixel count parameter in step S122). Using the region labeling information, the category label of each pixel is traversed to obtain the set of pixels Pk={(i1, j1), (i2, j2), ..., (in, jn)} contained in each superpixel region k (k=1..S), and the grayscale value sequence Gk=[G(i1, j1), G(i2, j2), ..., G(in, jn)] of these pixels is extracted, where n is the number of pixels in the region. For example, the 5th superpixel region contains 100 pixels, and its grayscale value sequence is 100 normalized grayscale values.

[0088] Step S132: Based on the gray value of each pixel in the superpixel segmentation result image, calculate the global mean vector and global covariance matrix of the gray values ​​of all pixels. Use the global mean vector as the initial value of the expected parameter of each Gaussian component in the Gaussian mixture model, and use the global covariance matrix as the initial value of the variance parameter of each Gaussian component in the Gaussian mixture model.

[0089] The global mean vector μ_global = mean(G(i,j))forall(i,j) is a scalar (because grayscale is a single channel). The global covariance matrix Σ_global = var(G(i,j)) is also a scalar. During Gaussian mixture model initialization, the expected parameter μk_initial = μ_global, the variance parameter σk_initial^2 = Σ_global, and the prior probability parameter πk_initial = 1 / S (equal probability initialization) for each Gaussian component k. For example, if the global mean is 0.5 and the global variance is 0.04, then the initial mean of all Gaussian components is 0.5, the initial variance is 0.04, and the initial prior probability is 1 / S.

[0090] Step S133: Perform statistical feature extraction processing on the pixel grayscale value sequence in each superpixel region, and calculate the regional mean vector and regional covariance matrix of all pixel grayscale values ​​in each superpixel region.

[0091] For a superpixel region k, the region mean μk_region=mean(Gk) and the region covariance σk_region^2=var(Gk) are used to assist in subsequent parameter updates. For example, in the M-step, they can be used as prior information for parameter updates.

[0092] Step S134: Based on the current expectation parameter, current variance parameter, and current prior probability parameter of each Gaussian component in the Gaussian mixture model, perform Gaussian probability density calculation processing on each pixel in the superpixel segmentation result image. By substituting the gray value of each pixel into the Gaussian probability density function defined by the current expectation parameter and the current variance parameter, calculate the probability density value generated by each Gaussian component for each pixel, and obtain the multidimensional probability density vector corresponding to each pixel.

[0093] For the gray value g of pixel (i, j), the probability density function of the k-th Gaussian component is p(g|μk, σk^2)=(1 / sqrt(2πσk^2))*exp(-(g-μk)^2 / (2σk^2)). Calculate the probability density values ​​of all K Gaussian components to obtain the multidimensional probability density vector p(i, j)=[p(g|μ1, σ1^2), p(g|μ2, σ2^2), ..., p(g|μK, σK^2)]. For example, when g=0.6, μk=0.5, σk^2=0.04, p(g|μk,σk^2)=(1 / sqrt(2π×0.04))*exp(-(0.6-0.5)^2 / (2×0.04))≈1.995×exp(-0.01 / 0.08)≈1.995×0.8825≈1.76.

[0094] Step S135: Based on the multidimensional probability density vector corresponding to each pixel and the current prior probability parameter of each Gaussian component, perform the expectation step calculation of the expectation maximization algorithm. By multiplying the current prior probability parameter of each Gaussian component by the probability density value generated by that Gaussian component for the corresponding pixel, and then dividing by the sum of the current prior probability parameters of all Gaussian components multiplied by the corresponding probability density values, calculate the belonging probability of each pixel to each Gaussian component, and obtain the belonging probability vector corresponding to each pixel.

[0095] The assignment probability γk(i,j) = (πk×p(g|μk,σk^2)) / sum_{m=1toK}(πm×p(g|μm,σm^2)), where γk(i,j) represents the probability that pixel (i,j) belongs to the k-th Gaussian component. The assignment probability vector for each pixel is γ(i,j) = [γ1(i,j), γ2(i,j), ..., γK(i,j)], satisfying sum(γk(i,j)) = 1. For example, if π1=0.3, p1=1.76, π2=0.7, p2=0.5, then γ1=(0.3×1.76) / (0.3×1.76+0.7×0.5)=0.528 / (0.528+0.35)=0.528 / 0.878≈0.601, γ2=1-0.601=0.399.

[0096] Step S136: Perform spatial smoothing filtering on the attribution probability vector corresponding to each pixel. By weighting and averaging the attribution probability vector of each pixel with the attribution probability vectors of all pixels in its neighboring pixel set, a spatially smoothed attribution probability vector corresponding to each pixel is generated.

[0097] The assignment probability is smoothed using the comprehensive influence weight w(i,j,i',j') calculated in step S129. For the assignment probability vector γ(i,j) of pixel (i,j), the spatially smoothed assignment probability vector γ_smooth(i,j)=sum_{(i',j')∈N(i,j)}w(i,j,i',j')×γ(i',j') / sum(w(i,j,i',j')). The influence of noise on the assignment probability is suppressed by neighborhood weighted averaging, thereby enhancing spatial consistency.

[0098] Step S137: The weights used in the weighted average are determined based on the spatial distance between pixels and the difference in grayscale values.

[0099] The same comprehensive influence weight calculation method as in step S129 ensures that neighboring pixels that are close in distance and have similar gray levels have a greater impact on the probability of the center pixel being assigned.

[0100] Step S138: Use the spatial smoothing assignment probability vector as the assignment probability vector used in the subsequent maximization step calculation.

[0101] Replace γ(i,j) with γ_smooth(i,j) for subsequent M-step parameter updates to improve the stability of parameter estimation.

[0102] Step S139: Based on the assignment probability vector corresponding to each pixel, perform the maximization step calculation of the expectation maximization algorithm, sum the assignment probability values ​​of all pixels belonging to the Gaussian component in the assignment probability vector corresponding to each Gaussian component, and obtain the total assignment weight of each Gaussian component.

[0103] The total assignment weight Nk=sum_{(i,j)}γk(i,j) is the sum of the probabilities of all pixels belonging to the k-th Gaussian component, reflecting the overall contribution of this component to the image pixels.

[0104] Step S1310: Calculate the new expected parameters for spatial smoothing of each Gaussian component based on the total spatial smoothing weight of each Gaussian component and the gray values ​​of all pixels belonging to that Gaussian component.

[0105] The total spatial smoothing assignment weight is Nk_smooth=sum_{(i,j)}γ_smoothk(i,j), and the new spatial smoothing expectation parameter is μk_new=sum_{(i,j)}(γ_smoothk(i,j)×G(i,j)) / Nk_smooth. The mean parameter is updated by weighted averaging, and the weight is the spatial smoothing assignment probability.

[0106] Step S1311: Calculate the new spatial smoothing variance parameter for each Gaussian component based on the new spatial smoothing expectation parameter for each Gaussian component, the grayscale values ​​of all pixels belonging to that Gaussian component, and the corresponding spatial smoothing assignment probability values.

[0107] The spatial smoothing new variance parameter σk_new^2=sum_{(i,j)}(γ_smoothk(i,j)×(G(i,j)-μk_new)^2) / Nk_smooth, also uses the spatial smoothing assignment probability as the weight to calculate the weighted squared difference between the gray value and the new mean.

[0108] Step S1312: Calculate the new prior probability parameter of spatial smoothing for each Gaussian component based on the total spatial smoothing weight of each Gaussian component and the total number of pixels in the superpixel segmentation result image.

[0109] The total number of pixels is N = M × N. The new prior probability parameter for spatial smoothing is πk_new = Nk_smooth / N, ensuring that sum(πk_new) = 1.

[0110] Step S1313: Construct the spatially smoothed intermediate Gaussian mixture model based on the new spatial smoothing expectation parameter, new spatial smoothing variance parameter, and new spatial smoothing prior probability parameter for each Gaussian component.

[0111] The parameters of the intermediate Gaussian mixture model are {(μk_new, σk_new^2, πk_new)|k=1..K}, which incorporates spatial smoothness information and has better noise resistance.

[0112] Step S1314: Replace the intermediate Gaussian mixture model with the spatially smoothed intermediate Gaussian mixture model for subsequent comparison of parameter changes and convergence determination.

[0113] In subsequent iterations, the spatially smoothed model parameters are used for convergence judgment to ensure the stability of parameter estimation.

[0114] Step S1315: Calculate the new expected parameter for each Gaussian component based on the total assignment weight of each Gaussian component and the gray values ​​of all pixels belonging to that Gaussian component. The new expected parameter is obtained by weighting the gray values ​​of all pixels belonging to that Gaussian component with their corresponding assignment probabilities as weights.

[0115] Without spatial smoothing, the new expected parameter μk_new=sum_{(i,j)}(γk(i,j)×G(i,j)) / Nk is similar to step S1310, but uses the original attribution probability γk(i,j).

[0116] Step S1316: Calculate the new variance parameter of each Gaussian component based on the new expected parameter of each Gaussian component, the gray values ​​of all pixels belonging to that Gaussian component, and the corresponding belonging probability values. The new variance parameter is obtained by weighting the squared difference between the gray values ​​of all pixels belonging to that Gaussian component and the new expected parameter, with the corresponding belonging probability values ​​as weights.

[0117] The new variance parameter σk_new^2=sum_{(i,j)}(γk(i,j)×(G(i,j)-μk_new)^2) / Nk is calculated using the original attribution probability.

[0118] Step S1317: Calculate a new prior probability parameter for each Gaussian component based on the total assigned weight of each Gaussian component and the total number of pixels in the superpixel segmentation result image. The new prior probability parameter is obtained by dividing the total assigned weight of each Gaussian component by the total number of pixels in the superpixel segmentation result image.

[0119] The new prior probability parameter πk_new=Nk / N is calculated using the original total attribution weights.

[0120] Step S1318: Based on the new expectation parameter, new variance parameter, and new prior probability parameter of each Gaussian component, construct an intermediate Gaussian mixture model to characterize the overall pixel distribution characteristics of the superpixel segmentation result image in the current iteration step, and calculate the log-likelihood function value corresponding to the intermediate Gaussian mixture model.

[0121] The intermediate Gaussian mixture model parameters are {(μk_new, σk_new^2, πk_new)|k=1..K}, and the log-likelihood function value LL=sum_{(i,j)}log(sum_{k=1toK}(πk_new×p(G(i,j)|μk_new, σk_new^2))) is used to evaluate the model's fit to the data.

[0122] Step S1319: Compare the new expected parameters, new variance parameters, and new prior probability parameters of each Gaussian component in the intermediate Gaussian mixture model with the expected parameters, variance parameters, and prior probability parameters of each Gaussian component in the Gaussian mixture model in the previous iteration step, and calculate the sum of the absolute values ​​of the changes in all parameters.

[0123] The parameter change Δ = sum_{k=1toK}(|μk_new-μk_old|+|σk_new^2-σk_old^2|+|πk_new-πk_old|) measures the magnitude of the parameter update.

[0124] Step S1320: Determine whether the sum of the absolute values ​​of the changes of all parameters is less than the preset expectation-maximization algorithm convergence threshold. If it is less than the preset expectation-maximization algorithm convergence threshold, stop the iteration. If it is not less than the preset expectation-maximization algorithm convergence threshold, use the new expectation parameter, new variance parameter, and new prior probability parameter of each Gaussian component in the intermediate Gaussian mixture model as the expectation parameter, variance parameter, and prior probability parameter of each Gaussian component in the current Gaussian mixture model, and repeat the Gaussian probability density calculation, expectation step calculation, maximization step calculation, and log-likelihood function numerical calculation until the convergence condition is reached.

[0125] The convergence threshold of the expectation-maximization algorithm is set to 1e-3. Iteration stops when Δ < 1e-3. Otherwise, let μk_old = μk_new, σk_old^2 = σk_new^2, πk_old = πk_new, and return to step S134 to continue iteration.

[0126] Step S1321: Determine the expected parameter, variance parameter, and prior probability parameter of each Gaussian component in the intermediate Gaussian mixture model when the convergence condition is met as the final Gaussian mixture model parameters. Construct a corresponding one-dimensional Gaussian probability density function based on the expected parameter and variance parameter of each Gaussian component in the final Gaussian mixture model parameters. Combine the one-dimensional Gaussian probability density functions of all Gaussian components in descending order of prior probability parameters to generate the set of multi-dimensional Gaussian probability density functions used to characterize the overall pixel distribution characteristics of the superpixel segmentation result image.

[0127] The final Gaussian mixture model parameters are {(μk_final, σk_final^2, πk_final)|k=1..K}, and the one-dimensional Gaussian probability density function corresponding to each component is p_k(g)=(1 / sqrt(2πσk_final^2))*exp(-(g-μk_final)^2 / (2σk_final^2)). These functions are sorted in descending order of πk_final to form a set of multi-dimensional Gaussian probability density functions, for example [p_1(g), p_2(g), ..., p_K(g)], where π1_final≥π2_final≥...≥πK_final.

[0128] Step S140: Perform Bayesian posterior probability calculation on each pixel in the superpixel segmentation result image according to the set of multidimensional Gaussian probability density functions to generate the posterior probability distribution matrix of each pixel belonging to each Gaussian component.

[0129] Each function in the set of multidimensional Gaussian probability density functions corresponds to the grayscale distribution characteristics of a land cover category. Bayesian posterior probability calculation aims to infer the probability of a pixel belonging to each land cover category based on its grayscale value. The posterior probability distribution matrix organizes the category probability of each pixel in matrix form. In urban remote sensing images, this posterior probability distribution matrix can clearly show the probability that each pixel belongs to categories such as buildings, roads, and vegetation. For example, the posterior probability distribution of a pixel might be [0.8 (buildings), 0.15 (roads), 0.05 (vegetation)], indicating that it is likely to be a building.

[0130] Step S141: Read the expected parameter, variance parameter and prior probability parameter corresponding to each Gaussian component in the multidimensional Gaussian probability density function set, and arrange each Gaussian component in descending order of prior probability parameter to generate a sorted Gaussian component sequence.

[0131] Extract the parameters (μk, σk^2, πk) of each Gaussian component k from the set of multidimensional Gaussian probability density functions, and sort them in descending order of πk to obtain the sequence [(μ1, σ1^2, π1), (μ2, σ2^2, π2), ..., (μK, σK^2, πK)], where π1≥π2≥...≥πK. For example, the first three sorted components may correspond to the three largest land cover types in the image.

[0132] Step S142: Obtain the grayscale value of each pixel in the superpixel segmentation result image, and arrange the grayscale values ​​of each pixel in the spatial coordinate order of the superpixel segmentation result image to form a grayscale value matrix of the pixels to be processed.

[0133] The gray value matrix G of the pixels to be processed is a two-dimensional matrix of M×N, which has the same spatial size as the superpixel segmentation result image. G_to_process(i,j) = G_superpixel(i,j), which is the gray value of the pixel in the i-th row and j-th column.

[0134] Step S143: For each pixel in the grayscale value matrix of the pixels to be processed, substitute it sequentially into the Gaussian probability density function corresponding to each Gaussian component in the sorted Gaussian component sequence. Subtract the expected parameter of the Gaussian component from the grayscale value of the pixel to obtain the difference. Calculate the square of the difference and divide it by twice the variance parameter of the Gaussian component to obtain the exponent part. Then multiply it by the normalization coefficient to calculate the conditional probability density value generated by the Gaussian component for the pixel, thus obtaining the conditional probability density value matrix of each pixel corresponding to each Gaussian component.

[0135] For the gray value g = G(i,j) of pixel (i,j), the conditional probability density p_k(g) of the k-th Gaussian component is calculated as follows: difference d = g - μk; exponential part e = d^2 / (2×σk^2); normalization coefficient c = 1 / sqrt(2×π×σk^2); p_k(g) = c×exp(-e). The conditional probability density values ​​of all pixels constitute the conditional probability density numerical matrix P, with dimensions M×N×K, where P(i,j,k) = p_k(g).

[0136] Step S144: Based on the conditional probability density matrix of each pixel corresponding to each Gaussian component and the prior probability parameter of each Gaussian component, perform Bayesian formula calculation for each pixel. Multiply the prior probability parameter of each Gaussian component by the conditional probability density value of the pixel corresponding to that Gaussian component to obtain the numerator. Sum the products of the prior probability parameters of all Gaussian components multiplied by the conditional probability density values ​​of the pixel corresponding to their respective Gaussian components to obtain the denominator. Divide the numerator by the denominator to obtain the posterior probability value of the pixel belonging to that Gaussian component.

[0137] Bayes' theorem states that the posterior probability P(k|g) = (πk × p_k(g)) / sum_{m=1 to K}(πm × p_m(g)). For a pixel (i, j), calculate K posterior probability values ​​to obtain the posterior probability vector P(i, j, :) = [P(1|g), P(2|g), ..., P(K|g)].

[0138] Step S145: Repeat the Bayesian formula calculation process for each pixel until the posterior probability value of each pixel belonging to all Gaussian components is calculated, and the posterior probability vector corresponding to each pixel is obtained.

[0139] Iterate through all pixels (i, j) (i=1..M, j=1..N) in the grayscale value matrix of the pixels to be processed, and perform the calculation of step S144 for each pixel to obtain M×N posterior probability vectors.

[0140] Step S146: Arrange the posterior probability vectors corresponding to all pixels in the spatial coordinate order of the superpixel segmentation result image to construct a three-dimensional posterior probability distribution matrix. The first dimension of the posterior probability distribution matrix corresponds to the vertical coordinate of the pixel, the second dimension of the posterior probability distribution matrix corresponds to the horizontal coordinate of the pixel, and the third dimension of the posterior probability distribution matrix corresponds to the index of the Gaussian component.

[0141] The posterior probability distribution matrix P_dist has dimensions M×N×K, where P_dist(i,j,k)=P(k|g(i,j)), representing the posterior probability that the pixel at vertical coordinate i and horizontal coordinate j belongs to the k-th Gaussian component.

[0142] For example, after arranging the posterior probability vectors corresponding to all pixels according to the spatial coordinate order of the superpixel segmentation result image to construct a three-dimensional posterior probability distribution matrix, the method further includes: Step S147: Perform principal component analysis on the three-dimensional posterior probability distribution matrix along the third dimension, calculate the covariance matrix of the posterior probability distribution matrix, and solve for the eigenvalues ​​and corresponding eigenvectors of the covariance matrix.

[0143] The posterior probability distribution matrix P_dist is reshaped into a two-dimensional matrix X (with dimensions of (M×N)×K), where each row corresponds to the posterior probability vector of a pixel. The covariance matrix C=cov(X) (K×K dimension) is obtained through eigenvalue decomposition C=VΛV^T, yielding the eigenvalue matrix Λ (a diagonal matrix with diagonal elements of λ1≥λ2≥...≥λK) and the eigenvector matrix V (column vectors are eigenvectors v1, v2, ..., vK).

[0144] Step S148: Sort the corresponding feature vectors according to the feature values ​​from largest to smallest, and select the feature vectors corresponding to the two largest feature values ​​as the principal component projection directions.

[0145] v1 and v2 (corresponding to λ1 and λ2) are selected as the principal component projection directions, which contain the main information of the posterior probability distribution.

[0146] Step S149: Project the posterior probability vector corresponding to each pixel in the three-dimensional posterior probability distribution matrix onto the principal component projection direction to generate a two-dimensional principal component score vector corresponding to each pixel.

[0147] For a pixel's posterior probability vector x (K-dimensional), the two-dimensional principal component score vector s = [v1^Tx, v2^Tx], where v1^Tx and v2^Tx are the projection values ​​of x in the v1 and v2 directions, respectively.

[0148] Step S1410: Construct a two-dimensional principal component score distribution matrix based on the two-dimensional principal component score vector corresponding to each pixel. The first dimension of the two-dimensional principal component score distribution matrix corresponds to the vertical coordinate of the pixel, the second dimension corresponds to the horizontal coordinate of the pixel, and the third dimension corresponds to the two principal component scores.

[0149] The two-dimensional principal component score distribution matrix S_dist has a dimension of M×N×2, where S_dist(i,j,1)=v1^TP_dist(i,j,:) and S_dist(i,j,2)=v2^TP_dist(i,j,:).

[0150] Step S1411: Perform visualization mapping processing on the two-dimensional principal component score distribution matrix, mapping the two principal component scores of each pixel in the two-dimensional principal component score distribution matrix to the red channel value and green channel value of the color image respectively, and setting the blue channel value to zero, to generate a pseudo-color posterior probability distribution image.

[0151] The principal component scores are normalized to the interval [0, 255]: s1_norm = 255 × (s1 - s1_min) / (s1_max - s1_min), s2_norm = 255 × (s2 - s2_min) / (s2_max - s2_min). The red channel R = s1_norm, the green channel G = s2_norm, and the blue channel B = 0 of the pseudo-color image, forming a pseudo-color posterior probability distribution image that visually displays the spatial distribution of the posterior probability.

[0152] Step S1412: The pseudo-color posterior probability distribution image and the superpixel segmentation result image are superimposed and fused. The pseudo-color posterior probability distribution image and the superpixel segmentation result image are weighted and summed according to a preset transparency parameter to generate a fused display image containing posterior probability visualization information.

[0153] The transparency parameter α (value 0.5) is used to control the superposition weight. The fused image F = α × pseudo-color image + (1-α) × superpixel segmentation result image (grayscale image converted to RGB format) makes the superpixel region boundary and posterior probability distribution features visualized simultaneously.

[0154] Step S1413: Perform maximum value extraction processing on the posterior probability vector corresponding to each pixel in the posterior probability distribution matrix, find the element with the largest value in each posterior probability vector, and record the Gaussian component index corresponding to the largest element to generate the maximum posterior probability class index matrix for each pixel.

[0155] The maximum a posteriori probability class index matrix Idx is an M×N two-dimensional matrix, Idx(i,j)=argmax(P_dist(i,j,:)), which is the Gaussian component index with the largest a posteriori probability, representing the most likely class of the pixel.

[0156] Step S1414: Based on the maximum posterior probability category index matrix and the posterior probability distribution matrix, normalize and verify all elements in the posterior probability vector corresponding to each pixel, calculate whether the sum of all elements in each posterior probability vector is equal to one, and if there are pixels whose sum deviates from one, then normalize the posterior probability vector of that pixel again.

[0157] For each pixel (i, j), calculate sum(P_dist(i, j, :)). If |sum-1|>1e-6 (numerical error threshold), then perform renormalization: P_dist(i, j, k)=P_dist(i, j, k) / sum, ensuring that the posterior probability vector satisfies the probability normalization condition.

[0158] Step S1415: Determine the posterior probability distribution matrix after normalization verification as the final posterior probability distribution matrix.

[0159] After normalization verification, each posterior probability vector in P_dist satisfies sum(P_dist(i,j,:))=1, which can be used for subsequent edge-preserving iterative optimization.

[0160] Step S150: Perform edge-preserving iterative optimization processing based on image gradient information according to the posterior probability distribution matrix and the superpixel segmentation result image. Construct a comprehensive objective function that includes a data fidelity term, a neighborhood space constraint term, and an image gradient regularization term. Determine the optimization step size parameter that minimizes the comprehensive objective function through a linear search method. Iteratively update the pixel values ​​of the superpixel segmentation result image according to the optimization step size parameter to generate the final processed result image after edge enhancement and noise suppression.

[0161] Superpixel segmentation results may suffer from intra-region noise and blurred region boundaries. Edge-preserving iterative optimization constructs a comprehensive objective function to suppress noise while preserving image edge information. A data fidelity term ensures consistency between the optimized result and the original superpixel image, a neighborhood space constraint term promotes smoothness within pixels, and an image gradient regularization term enhances gradient values ​​in edge regions to highlight boundaries. A linear search method is used to determine the optimal iteration step size, enabling the objective function to converge quickly to its minimum. In urban remote sensing images, this step can make building edges clearer, road areas smoother, and remove speckle noise from vegetation areas, improving image visual quality and the reliability of subsequent analysis.

[0162] Step S151: Read the current pixel value matrix and the posterior probability distribution matrix of the superpixel segmentation result image, and normalize the current pixel value matrix.

[0163] The current pixel value matrix I_current is an M×N two-dimensional matrix, with elements representing the gray values ​​of the superpixel segmentation result image (in the [0, 255] interval). Normalization transforms it to the [0, 1] interval: I_current_norm(i, j) = I_current(i, j) / 255. The posterior probability distribution matrix P_dist is the M×N×K dimensional matrix obtained in step S1415.

[0164] Step S152: Construct a data fidelity term in the comprehensive objective function based on the current pixel value matrix and the posterior probability distribution matrix. The data fidelity term is calculated by weighting and combining the gray value of each pixel in the current pixel value matrix with the posterior probability vector of the corresponding pixel in the posterior probability distribution matrix, and calculating the sum of squared differences between the estimated value after weighting and the original observation value.

[0165] The data fidelity term E_data = sum_{i,j}(I_current_norm(i,j)-sum_{k=1toK}(P_dist(i,j,k)×μk_final))^2, where sum_{k}(P_dist(i,j,k)×μk_final) is the weighted average gray-level estimate based on the posterior probability, and μk_final is the final expectation parameter of the Gaussian component. This data term measures the difference between the current image and the expected image based on the probabilistic model.

[0166] Step S153: Construct a neighborhood spatial constraint term based on the current pixel value matrix and the spatial adjacency relationship between pixels in the superpixel segmentation result image. The neighborhood spatial constraint term is obtained by calculating the sum of the squares of the gray value differences between each pixel and all pixels in its neighborhood pixel set, and then accumulating the calculation results for all pixels.

[0167] An 8-neighborhood system is adopted, with the neighborhood space constraint term E_space=sum_{i,j}sum_{(i',j')∈N(i,j)}(I_current_norm(i,j)-I_current_norm(i',j'))^2. By penalizing the grayscale difference between adjacent pixels, smoothness within the region is promoted.

[0168] Step S154: Perform gradient operator convolution processing on the current pixel value matrix, calculate the horizontal gradient matrix of the current pixel value matrix in the horizontal direction and the vertical gradient matrix in the vertical direction, construct an image gradient regularization term based on the horizontal gradient matrix and the vertical gradient matrix. The image gradient regularization term is obtained by calculating the sum of the squares of each element in the horizontal gradient matrix and the squares of each element in the vertical gradient matrix, and accumulating the calculation results of all elements.

[0169] The horizontal gradient operator uses [-1, 0, 1], and the horizontal gradient matrix Gx(i, j) = I_current_norm(i, j+1) - I_current_norm(i, j-1) (boundary pixels are zero-padded); the vertical gradient operator uses [-1; 0; 1], and the vertical gradient matrix Gy(i, j) = I_current_norm(i+1, j) - I_current_norm(i-1, j). The image gradient regularization term E_grad = sum_{i, j}(Gx(i, j)^2 + Gy(i, j)^2) enhances the gradient values ​​to highlight edge features.

[0170] Step S155: Before constructing the comprehensive objective function, the current pixel value matrix is ​​normalized. Based on the normalized pixel value matrix, the data fidelity term, the neighborhood space constraint term, and the image gradient regularization term are calculated. The comprehensive objective function is constructed based on the data fidelity term, the neighborhood space constraint term, and the image gradient regularization term calculated after normalization. The comprehensive objective function is the sum of the product of the data fidelity term, the neighborhood space constraint term, and the first weight coefficient, and the product of the image gradient regularization term and the second weight coefficient.

[0171] The overall objective function is E = E_data + α × E_space + β × E_grad, where α is the neighborhood space constraint weight (set to 0.1) and β is the image gradient regularization weight (set to 0.05). The weight parameters are adjusted according to the image noise level and edge intensity; increasing α enhances the smoothing effect, and increasing β enhances the edge enhancement effect.

[0172] Step S156: Set the initial iteration variable of the comprehensive objective function to the current pixel value matrix, and set an initial step size parameter candidate set for linear search, wherein the initial step size parameter candidate set contains multiple candidate step size values ​​distributed according to a preset interval.

[0173] The initial iteration variable is I = I_current_norm, and the initial step size parameter candidate set is τ_candidates = [0.001, 0.005, 0.01, 0.05, 0.1, 0.2] with logarithmic distribution, covering step sizes of different magnitudes.

[0174] Step S157: For each candidate step size value in the candidate set of initial step size parameters, subtract the candidate step size value multiplied by the gradient matrix of the comprehensive objective function with respect to the current pixel value matrix from the current pixel value matrix to obtain the corresponding candidate updated pixel value matrix.

[0175] The gradient of the overall objective function, ∇E = ∂E / ∂I, is obtained by differentiating E (the specific derivative calculation involves the partial derivatives of each energy term). The candidate updated pixel value matrix is ​​I_candidate = I - τ × ∇E, where τ is the candidate step size.

[0176] Step S158: Calculate the candidate comprehensive objective function value corresponding to each candidate updated pixel value matrix by substituting each candidate updated pixel value matrix into the expression of the comprehensive objective function.

[0177] For each I_candidate, calculate E_candidate=E_data(I_candidate)+α×E_space(I_candidate)+β×E_grad(I_candidate).

[0178] Step S159: Compare the candidate synthesis objective function values ​​corresponding to all candidate step size values ​​in the initial step size parameter candidate set, find the candidate step size value that minimizes the candidate synthesis objective function value, and determine the candidate step size value as the optimized step size parameter for the current iteration step.

[0179] Optimize the step size parameter τ_opt=argmin(E_candidate(τ)), and select the τ value that minimizes E_candidate.

[0180] Step S1510: Update the current pixel value matrix according to the optimization step size parameter and the gradient matrix of the comprehensive objective function with respect to the current pixel value matrix to generate a new pixel value matrix. The new pixel value matrix is ​​obtained by subtracting the optimization step size parameter from the current pixel value matrix and multiplying it by the gradient matrix.

[0181] The new pixel value matrix I_new = I - τ_opt × ∇E.

[0182] Step S1511: Calculate the difference metric between the new pixel value matrix and the current pixel value matrix, and compare the difference metric with a preset iteration termination threshold. If the difference metric is less than the preset iteration termination threshold, it is determined that the convergence condition has been met and the iteration stops. If the difference metric is not less than the preset iteration termination threshold, the new pixel value matrix is ​​used as the current pixel value matrix, and the process of constructing the comprehensive objective function, setting the initial step size parameter candidate set, calculating the candidate updated pixel value matrix, determining the optimized step size parameter, and updating is repeated until the convergence condition is met.

[0183] The difference metric D = sum_{i,j}|I_new(i,j)-I(i,j)| / (M×N) is set, and the iteration termination threshold is set to 1e-4. If D < 1e-4, the iteration stops; otherwise, I = I_new, and the iteration returns to step S152.

[0184] Step S1512: Perform inverse normalization on the new normalized pixel value matrix that has reached the convergence condition, and determine the inverse normalized pixel value matrix as the intermediate result image after edge enhancement and noise suppression.

[0185] Inverse normalization: I_mid(i,j)=I_new(i,j)×255, convert back to the gray range of [0,255] to obtain the intermediate result image.

[0186] Step S1513: Perform edge sharpening enhancement processing on the intermediate result image. Calculate the horizontal and vertical second derivative matrices of the intermediate result image, construct a Laplacian sharpening operator based on the horizontal and vertical second derivative matrices, multiply the Laplacian sharpening operator by a preset sharpening intensity coefficient, and then superimpose it onto the intermediate result image to generate the final processed result image.

[0187] The horizontal second derivative matrix is ​​Lx(i,j) = I_mid(i,j+1) - 2 × I_mid(i,j) + I_mid(i,j-1); the vertical second derivative matrix is ​​Ly(i,j) = I_mid(i+1,j) - 2 × I_mid(i,j) + I_mid(i-1,j); the Laplacian operator is L = Lx + Ly. The sharpening intensity coefficient k = 0.5, and the final processed image is I_final(i,j) = I_mid(i,j) - k × L(i,j). By subtracting the Laplacian operator to enhance edge contrast, the final adaptive processing result of the remote sensing image is generated.

[0188] In one exemplary embodiment, a remote sensing image adaptive processing system combining optimization algorithms and probabilistic models is provided. This system can be a terminal, server, etc., and its internal structure diagram can be as follows: Figure 2 As shown, the remote sensing image adaptive processing system combining optimization algorithms and probabilistic models includes a processor, memory, input / output interface, communication interface, display unit, and input device. The processor, memory, and input / output interface are connected via a system bus, and the communication interface, display unit, and input device are also connected to the system bus via the input / output interface. The processor provides computational and control capabilities. The memory includes a non-volatile storage medium and internal memory. The non-volatile storage medium stores the operating system and computer programs. The internal memory provides the environment for the operation of the operating system and computer programs in the non-volatile storage medium. The input / output interface is used for exchanging information between the processor and external devices. The communication interface is used for wired or wireless communication with external terminals; wireless communication can be achieved through Wi-Fi, mobile cellular networks, near-field communication, or other technologies. When the computer program is executed by the processor, it implements a remote sensing image adaptive processing method combining optimization algorithms and probabilistic models. The display unit is used to form a visually visible image and can be a display screen, projection device, or virtual reality imaging device. The display screen can be an LCD screen or an e-ink screen. The input device can be a touch layer covering the display screen, or a button, trackball, or touchpad set on the shell of a remote sensing image adaptive processing system that combines optimization algorithms and probability models. It can also be an external keyboard, touchpad, or mouse, etc.

[0189] It should be noted that, in order to simplify the description of the present invention and thus help to understand one or more embodiments of the invention, multiple features may sometimes be grouped into one embodiment, drawing or description thereof in the foregoing description of the embodiments of the present invention.

Claims

1. A remote sensing image adaptive processing method combining optimization algorithms and probabilistic models, characterized in that, The method includes: The original projection data set corresponding to the acquired original remote sensing image is subjected to projection domain data reconstruction processing based on the Landweber iterative framework. The original projection data set is successively approximated and solved according to the preset relaxation factor parameters and the preset convergence judgment threshold to generate an initial density distribution image. An adaptive segmentation process based on Gaussian mixture model and neighborhood correlation constraint is performed on the initial density distribution image. Markov random field theory is introduced to construct a neighborhood dependency weight for each pixel in the initial density distribution image to measure its consistency with the neighboring pixel labels. In the iterative segmentation framework based on Gaussian mixture model and Markov random field coupling, the Gibbs energy function is constructed using the neighborhood dependency weight to perform pixel category classification processing on the initial density distribution image and generate a superpixel segmentation result image containing smooth region boundaries. The superpixel segmentation result image is subjected to Gaussian mixture model parameter estimation processing based on the expectation-maximization algorithm. The expectation parameter, variance parameter and prior probability parameter of each Gaussian component in the current Gaussian mixture model are calculated using the pixel statistics of each superpixel region in the superpixel segmentation result image. A set of multidimensional Gaussian probability density functions is constructed based on the expectation parameter, the variance parameter and the prior probability parameter to characterize the overall pixel distribution characteristics of the superpixel segmentation result image. Based on the set of multidimensional Gaussian probability density functions, Bayesian posterior probability calculation is performed on each pixel in the superpixel segmentation result image to generate a posterior probability distribution matrix for each pixel belonging to each Gaussian component. Based on the posterior probability distribution matrix and the superpixel segmentation result image, perform edge-preserving iterative optimization processing based on image gradient information. Construct a comprehensive objective function that includes a data fidelity term, a neighborhood space constraint term, and an image gradient regularization term. Determine the optimization step size parameter that minimizes the comprehensive objective function through a linear search method. Iteratively update the pixel values ​​of the superpixel segmentation result image according to the optimization step size parameter to generate the final processed result image after edge enhancement and noise suppression.

2. The remote sensing image adaptive processing method combining optimization algorithms and probabilistic models according to claim 1, characterized in that, The process of performing projection domain data reconstruction based on the Landweber iterative framework on the original projection data set corresponding to the acquired original remote sensing image includes: successive approximation solutions are performed on the original projection data set according to preset relaxation factor parameters and preset convergence thresholds to generate an initial density distribution image. The original projection data set of the original remote sensing image is acquired by the detector array during the imaging process. The original projection data set contains a sequence of ray attenuation intensity values ​​recorded by multiple detector channels at different rotation angles. A system response matrix is ​​constructed based on the geometric parameters of the imaging system to describe the linear mapping relationship between the original projection data set and the image to be reconstructed. The row index of the system response matrix corresponds to the ray path of the detector channel, and the column index of the system response matrix corresponds to the pixel unit position in the image space to be reconstructed. The initial iterative solution vector of the Landweber iterative framework is set, each element of the initial iterative solution vector corresponds to the initial density estimate of each pixel unit in the image to be reconstructed, and all elements of the initial iterative solution vector are assigned to a preset zero initial state, and the initial iterative solution vector is normalized preprocessed. The Landweber iterative framework is configured with a relaxation factor parameter and a convergence threshold. The relaxation factor parameter controls the magnitude of the correction to the current solution vector along the negative gradient direction during each iteration. The convergence threshold is used to determine whether the iteration process has reached the termination condition. The residual vector for the current iteration is calculated based on the system response matrix, the original projection data set, the relaxation factor parameter, and the current iteration solution vector. The residual vector is obtained by multiplying the system response matrix with the current iteration solution vector to obtain the current estimated projection data, and then subtracting the current estimated projection data from the corresponding projection value in the original projection data set element by element. Before executing the Landweber iterative framework, the original projected data set and the initial iterative solution vector are normalized and preprocessed. During the iteration process, the normalized current iterative solution vector is updated according to the transpose of the system response matrix, the relaxation factor parameter, and the residual vector to generate a new normalized iterative solution vector. The new normalized iterative solution vector is obtained by adding the product of the relaxation factor parameter and the transpose of the system response matrix multiplied by the residual vector to the normalized current iterative solution vector. Calculate the norm of the new residual vector corresponding to the new normalized iterative solution vector, and compare the norm of the new residual vector with the convergence determination threshold. If the norm of the new residual vector is less than the convergence determination threshold, it is determined that the convergence condition has been met and the iteration process of the Landweber iterative framework is stopped. If the norm of the new residual vector is greater than or equal to the convergence determination threshold, the new normalized iterative solution vector is used as the current iterative solution vector to continue the next iteration update. Record the final normalized iterative solution vector when the convergence condition is met, perform inverse normalization on the final normalized iterative solution vector, and reorganize the inverse normalized final iterative solution vector according to the two-dimensional arrangement order of pixel units in the image space to generate the initial density distribution image with spatial dimension information.

3. The remote sensing image adaptive processing method combining optimization algorithms and probabilistic models according to claim 1, characterized in that, The process involves performing adaptive segmentation on the initial density distribution image based on a Gaussian mixture model and neighborhood correlation constraints. Markov random field theory is introduced to construct a neighborhood dependency weight for each pixel in the initial density distribution image, measuring its consistency with neighboring pixel labels. Within an iterative segmentation framework coupled with the Gaussian mixture model and Markov random field, the neighborhood dependency weight is used to construct a Gibbs energy function to classify the initial density distribution image into pixel categories, generating a superpixel segmentation result image containing smooth region boundaries. This includes: Read all pixels of the initial density distribution image and their corresponding gray values, and convert the initial density distribution image into a two-dimensional gray matrix. The row index of the two-dimensional gray matrix corresponds to the vertical coordinate of the pixel in the image space, and the column index of the two-dimensional gray matrix corresponds to the horizontal coordinate of the pixel in the image space. The initial cluster center point generation process is performed on the two-dimensional grayscale matrix according to the preset superpixel number parameter. Multiple initial cluster center point coordinates with the same value as the superpixel number parameter are generated within the range of pixel coordinate values ​​of the two-dimensional grayscale matrix by a random number generation function. Each initial cluster center point coordinate corresponds to the initial center position of a superpixel region to be generated. Using the initial cluster center coordinates as a reference, a comprehensive distance metric calculation is performed on each pixel in the two-dimensional grayscale matrix. The normalized spatial Euclidean distance and the normalized grayscale difference distance between each pixel and each initial cluster center are calculated. The normalized spatial Euclidean distance is calculated based on the difference in row and column coordinates of the pixel in its spatial location. The normalized grayscale difference distance is calculated based on the difference in grayscale values ​​of the pixel. The normalized spatial Euclidean distance and the normalized grayscale difference distance are linearly combined according to preset weights to obtain the comprehensive distance. Based on the comprehensive distance calculation results between each pixel and each initial cluster center point, each pixel is assigned to the superpixel region category corresponding to the initial cluster center point with the smallest Euclidean distance value, forming the initial superpixel region set. For each superpixel region in the initially divided superpixel region set, perform cluster center point update processing, calculate the average spatial coordinates and average gray value of all pixels in each superpixel region, and use the average spatial coordinates and average gray value as the new cluster center point coordinates; Repeat the comprehensive distance metric calculation and combination processing, the pixel point allocation processing, and the cluster center point update processing until the change in the coordinates of the cluster center points in two adjacent iterations is less than the preset clustering convergence threshold, and obtain the initial superpixel segmentation result based on the K-Means clustering algorithm. Based on the superpixel region category to which each pixel belongs in the initial superpixel segmentation result, an initial label vector is constructed for each pixel. The dimension of the initial label vector is equal to the total number of superpixel regions. Only the element corresponding to the superpixel region category to which the pixel belongs has a value of 1 in the initial label vector, and the other elements have a value of 0. A Markov random field neighborhood system is constructed based on the spatial adjacency relationship between each pixel and its surrounding neighboring pixels in the initial density distribution image. The Markov random field neighborhood system defines a set of neighboring pixels for each pixel, which includes the pixels directly adjacent to the pixel in spatial position above, below, to the left, and to the right, as well as four adjacent pixels in the diagonal direction. Based on Markov random field theory and Gibbs distribution, a neighborhood correlation constraint energy function is constructed. The neighborhood correlation constraint energy function is used to measure the degree of consistency between the current label vector of each pixel and the label vectors of all pixels in its neighborhood pixel set. The neighborhood dependency weight is generated for each pixel according to the neighborhood correlation constraint energy function. The neighborhood dependency weight is obtained by performing a negative exponential transformation on the neighborhood correlation constraint energy function and then normalizing it. The gray value of each pixel and the initial label vector are input into an iterative segmentation framework based on the coupling of Gaussian mixture model and Markov random field. The Gibbs energy function is constructed using the neighborhood dependency weights as the constraint condition for segmentation. The expectation step and the maximization step are executed alternately through the expectation-maximization algorithm to iteratively estimate the parameters of the Gaussian mixture model and the posterior probability of each pixel. Repeat the expectation step and the maximization step until the parameter change of the Gaussian mixture model is less than the preset model convergence threshold. Based on the posterior probability of each pixel belonging to each Gaussian component obtained in the final iteration, select the category corresponding to the Gaussian component with the largest posterior probability value as the final superpixel region category of the pixel. The two-dimensional grayscale matrix is ​​labeled according to the final superpixel region category of each pixel. Adjacent pixels with the same final superpixel region category are merged into the same superpixel region. The boundary of the merged superpixel region is smoothed to generate the superpixel segmentation result image containing the smoothed region boundary.

4. The remote sensing image adaptive processing method combining optimization algorithms and probabilistic models according to claim 1, characterized in that, The step involves performing Gaussian mixture model parameter estimation processing based on the expectation-maximization algorithm on the superpixel segmentation result image. This process calculates the expectation parameter, variance parameter, and prior probability parameter of each Gaussian component in the current Gaussian mixture model using the pixel statistics of each superpixel region in the superpixel segmentation result image. Based on the expectation parameter, variance parameter, and prior probability parameter, a set of multidimensional Gaussian probability density functions characterizing the overall pixel distribution of the superpixel segmentation result image is constructed, including: The superpixel segmentation result image is analyzed to obtain the total number of superpixel regions contained in the superpixel segmentation result image and the gray value sequence of all pixels contained in each superpixel region; Based on the grayscale value of each pixel in the superpixel segmentation result image, calculate the global mean vector and global covariance matrix of the grayscale values ​​of all pixels. Use the global mean vector as the initial value of the expected parameter of each Gaussian component in the Gaussian mixture model, and use the global covariance matrix as the initial value of the variance parameter of each Gaussian component in the Gaussian mixture model. Perform statistical feature extraction processing on the pixel grayscale value sequence within each superpixel region, and calculate the regional mean vector and regional covariance matrix of all pixel grayscale values ​​within each superpixel region; Based on the current expectation parameter, current variance parameter, and current prior probability parameter of each Gaussian component in the Gaussian mixture model, Gaussian probability density calculation is performed on each pixel in the superpixel segmentation result image. By substituting the gray value of each pixel into the Gaussian probability density function defined by the current expectation parameter and the current variance parameter, the probability density value generated by each Gaussian component for each pixel is calculated, and the multidimensional probability density vector corresponding to each pixel is obtained. Based on the multidimensional probability density vector corresponding to each pixel and the current prior probability parameter of each Gaussian component, the expectation step of the expectation maximization algorithm is performed. By multiplying the current prior probability parameter of each Gaussian component by the probability density value generated by that Gaussian component for the corresponding pixel, and then dividing by the sum of the current prior probability parameters of all Gaussian components multiplied by the corresponding probability density values, the probability of each pixel belonging to each Gaussian component is calculated, and the probability vector corresponding to each pixel is obtained. Based on the assignment probability vector corresponding to each pixel, the maximization step of the expectation maximization algorithm is performed to sum the assignment probability values ​​of all pixels belonging to that Gaussian component in the assignment probability vector corresponding to each Gaussian component, and obtain the total assignment weight of each Gaussian component. Based on the total assignment weight of each Gaussian component and the gray values ​​of all pixels belonging to that Gaussian component, a new expected parameter for each Gaussian component is calculated. The new expected parameter is obtained by weighting the gray values ​​of all pixels belonging to that Gaussian component with their corresponding assignment probabilities as weights. Based on the new expectation parameter of each Gaussian component, the gray values ​​of all pixels belonging to that Gaussian component, and the corresponding assignment probability values, the new variance parameter of each Gaussian component is calculated. The new variance parameter is obtained by weighting the squared difference between the gray values ​​of all pixels belonging to that Gaussian component and the new expectation parameter, with the corresponding assignment probability values ​​as weights. Based on the total assigned weight of each Gaussian component and the total number of pixels in the superpixel segmentation result image, a new prior probability parameter for each Gaussian component is calculated. The new prior probability parameter is obtained by dividing the total assigned weight of each Gaussian component by the total number of pixels in the superpixel segmentation result image. Based on the new expectation parameter, new variance parameter, and new prior probability parameter of each Gaussian component, an intermediate Gaussian mixture model is constructed to characterize the overall pixel distribution characteristics of the superpixel segmentation result image in the current iteration step, and the log-likelihood function value corresponding to the intermediate Gaussian mixture model is calculated. The new expected parameters, new variance parameters, and new prior probability parameters of each Gaussian component in the intermediate Gaussian mixture model are compared with the expected parameters, variance parameters, and prior probability parameters of each Gaussian component in the Gaussian mixture model in the previous iteration step, and the sum of the absolute values ​​of the changes of all parameters is calculated. Determine whether the sum of the absolute values ​​of the changes in all parameters is less than a preset convergence threshold for the expectation-maximization algorithm. If it is less than the preset convergence threshold, stop the iteration. If it is not less than the preset convergence threshold, use the new expectation parameter, new variance parameter, and new prior probability parameter of each Gaussian component in the intermediate Gaussian mixture model as the expectation parameter, variance parameter, and prior probability parameter of each Gaussian component in the current Gaussian mixture model, and repeat the Gaussian probability density calculation, expectation step calculation, maximization step calculation, and log-likelihood function numerical calculation until the convergence condition is met. The expected parameter, variance parameter, and prior probability parameter of each Gaussian component in the intermediate Gaussian mixture model when the convergence condition is met are determined as the parameters of the final Gaussian mixture model. A corresponding one-dimensional Gaussian probability density function is constructed based on the expected parameter and variance parameter of each Gaussian component in the final Gaussian mixture model parameters. The one-dimensional Gaussian probability density functions of all Gaussian components are combined in descending order of prior probability parameters to generate the set of multi-dimensional Gaussian probability density functions used to characterize the overall pixel distribution characteristics of the superpixel segmentation result image.

5. The remote sensing image adaptive processing method combining optimization algorithms and probabilistic models according to claim 1, characterized in that, The step of performing Bayesian posterior probability calculation on each pixel in the superpixel segmentation result image based on the multidimensional Gaussian probability density function set to generate a posterior probability distribution matrix for each pixel belonging to each Gaussian component includes: Read the expected parameter, variance parameter and prior probability parameter corresponding to each Gaussian component in the set of multidimensional Gaussian probability density functions, and arrange each Gaussian component in descending order of prior probability parameter to generate a sorted Gaussian component sequence. Obtain the grayscale value of each pixel in the superpixel segmentation result image, and arrange the grayscale values ​​of each pixel according to the spatial coordinate order of the superpixel segmentation result image to form a grayscale value matrix of the pixels to be processed. For each pixel in the grayscale value matrix of the pixels to be processed, the pixel is substituted into the Gaussian probability density function corresponding to each Gaussian component in the sorted Gaussian component sequence. The difference is obtained by subtracting the expected parameter of the Gaussian component from the grayscale value of the pixel. The square of the difference is calculated and divided by twice the variance parameter of the Gaussian component to obtain the exponent. Then, the exponent is multiplied by the normalization coefficient to calculate the conditional probability density value generated by the Gaussian component for the pixel. This yields the conditional probability density value matrix for each pixel corresponding to each Gaussian component. Based on the conditional probability density matrix of each Gaussian component corresponding to each pixel and the prior probability parameter of each Gaussian component, Bayes' theorem is performed on each pixel. The prior probability parameter of each Gaussian component is multiplied by the conditional probability density value of the pixel corresponding to that Gaussian component to obtain the numerator. The products of the prior probability parameters of all Gaussian components multiplied by the conditional probability density values ​​of the pixel corresponding to their respective Gaussian components are summed to obtain the denominator. The numerator is divided by the denominator to obtain the posterior probability value of the pixel belonging to that Gaussian component. The Bayesian formula calculation process is repeated for each pixel until the posterior probability value of each pixel belonging to all Gaussian components is calculated, thus obtaining the posterior probability vector corresponding to each pixel. Arrange the posterior probability vectors corresponding to all pixels according to the spatial coordinate order of the superpixel segmentation result image to construct a three-dimensional posterior probability distribution matrix. The first dimension of the posterior probability distribution matrix corresponds to the vertical coordinate of the pixel, the second dimension of the posterior probability distribution matrix corresponds to the horizontal coordinate of the pixel, and the third dimension of the posterior probability distribution matrix corresponds to the index of the Gaussian component. For each pixel in the posterior probability distribution matrix, perform maximum value extraction processing to find the element with the largest value in each posterior probability vector and record the Gaussian component index corresponding to the largest element to generate the maximum posterior probability class index matrix for each pixel. Based on the maximum posterior probability category index matrix and the posterior probability distribution matrix, normalization verification is performed on all elements in the posterior probability vector corresponding to each pixel. The sum of all elements in each posterior probability vector is calculated to see if it is equal to one. If there are pixels whose sum deviates from one, the posterior probability vector of that pixel is normalized again. The posterior probability distribution matrix after normalization verification is determined as the final posterior probability distribution matrix.

6. The remote sensing image adaptive processing method combining optimization algorithms and probabilistic models according to claim 1, characterized in that, The process involves performing edge-preserving iterative optimization based on image gradient information according to the posterior probability distribution matrix and the superpixel segmentation result image. This constructs a comprehensive objective function including a data fidelity term, a neighborhood space constraint term, and an image gradient regularization term. An optimization step size parameter that minimizes the comprehensive objective function is determined using a linear search method. The pixel values ​​of the superpixel segmentation result image are iteratively updated according to the optimization step size parameter to generate a final processed result image with edge enhancement and noise suppression. This includes: Read the current pixel value matrix and the posterior probability distribution matrix of the superpixel segmentation result image, and normalize the current pixel value matrix. Based on the current pixel value matrix and the posterior probability distribution matrix, a data fidelity term is constructed in the comprehensive objective function. The data fidelity term is calculated by weighting and combining the gray value of each pixel in the current pixel value matrix with the posterior probability vector of the corresponding pixel in the posterior probability distribution matrix, and calculating the sum of squared differences between the weighted combined estimated value and the original observed value. Based on the current pixel value matrix and the spatial adjacency relationship between pixels in the superpixel segmentation result image, a neighborhood spatial constraint term is constructed. The neighborhood spatial constraint term is obtained by calculating the sum of the squares of the gray value differences between each pixel and all pixels in its neighborhood pixel set, and then summing the calculation results for all pixels. Gradient operator convolution processing is performed on the current pixel value matrix to calculate the horizontal gradient matrix in the horizontal direction and the vertical gradient matrix in the vertical direction. An image gradient regularization term is constructed based on the horizontal gradient matrix and the vertical gradient matrix. The image gradient regularization term is obtained by calculating the sum of the squares of each element in the horizontal gradient matrix and the squares of each element in the vertical gradient matrix, and then accumulating the calculation results of all elements. Before constructing the comprehensive objective function, the current pixel value matrix is ​​normalized. Based on the normalized pixel value matrix, the data fidelity term, the neighborhood space constraint term, and the image gradient regularization term are calculated. The comprehensive objective function is constructed based on the data fidelity term, the neighborhood space constraint term, and the image gradient regularization term calculated after normalization. The comprehensive objective function is the sum of the product of the data fidelity term, the neighborhood space constraint term, and the first weight coefficient, and the product of the image gradient regularization term and the second weight coefficient. The initial iteration variable of the comprehensive objective function is set to the current pixel value matrix, and an initial step size parameter candidate set is set for linear search. The initial step size parameter candidate set contains multiple candidate step size values ​​distributed according to a preset interval. For each candidate step size value in the candidate set of initial step size parameters, subtract the candidate step size value from the current pixel value matrix by multiplying it by the gradient matrix of the comprehensive objective function with respect to the current pixel value matrix to obtain the corresponding candidate updated pixel value matrix; The candidate comprehensive objective function value corresponding to each candidate updated pixel value matrix is ​​calculated by substituting each candidate updated pixel value matrix into the expression of the comprehensive objective function. Compare the candidate comprehensive objective function values ​​corresponding to all candidate step size values ​​in the initial step size parameter candidate set, find the candidate step size value that minimizes the candidate comprehensive objective function value, and determine the candidate step size value as the optimized step size parameter for the current iteration step; Based on the optimization step size parameter and the gradient matrix of the comprehensive objective function with respect to the current pixel value matrix, the current pixel value matrix is ​​updated to generate a new pixel value matrix. The new pixel value matrix is ​​obtained by subtracting the optimization step size parameter from the current pixel value matrix and multiplying it by the gradient matrix. Calculate the difference metric between the new pixel value matrix and the current pixel value matrix, and compare the difference metric with a preset iteration termination threshold. If the difference metric is less than the preset iteration termination threshold, it is determined that the convergence condition has been met and the iteration stops. If the difference metric is not less than the preset iteration termination threshold, the new pixel value matrix is ​​used as the current pixel value matrix, and the process of constructing the comprehensive objective function, setting the initial step size parameter candidate set, calculating the candidate updated pixel value matrix, determining the optimized step size parameter, and updating is repeated until the convergence condition is met. The new normalized pixel value matrix that reaches the convergence condition is denormalized, and the denormalized pixel value matrix is ​​determined as the intermediate result image after edge enhancement and noise suppression. Edge sharpening enhancement processing is performed on the intermediate result image. The horizontal and vertical second derivative matrices of the intermediate result image are calculated, and a Laplacian sharpening operator is constructed based on the horizontal and vertical second derivative matrices. The Laplacian sharpening operator is multiplied by a preset sharpening intensity coefficient and then superimposed on the intermediate result image to generate the final processed result image.

7. The remote sensing image adaptive processing method combining optimization algorithm and probabilistic model according to claim 2, characterized in that, Before calculating the residual vector of the current iteration based on the system response matrix, the original projection data set, the relaxation factor parameters, and the current iteration solution vector, the method further includes: Condition number pre-analysis is performed on the system response matrix to calculate the ratio of the maximum singular value to the minimum singular value of the system response matrix as the condition number value. The condition number value is compared with a preset ill-condition threshold. If the condition number value is greater than the preset ill-condition threshold, the system response matrix is ​​determined to exhibit ill-condition characteristics. Based on the determination result of the ill-conditioned characteristics, a diagonal loading preprocessing operation is performed on the system response matrix to generate a corrected system response matrix after diagonal loading correction. The corrected system response matrix is ​​obtained by adding a preset diagonal loading coefficient to the system response matrix and multiplying it by the identity matrix. The corrected system response matrix is ​​used as the system response matrix for subsequent calculations of the current estimated projection data; The current estimated projection data after diagonal loading correction is calculated based on the corrected system response matrix and the current iterative solution vector. The current estimated projection data after diagonal loading correction is obtained by multiplying the corrected system response matrix with the current iterative solution vector. The diagonally loaded and corrected residual vector is obtained by subtracting the corresponding projection values ​​in the original projection data set element by element from the current estimated projection data after diagonal loading correction. The current iterative solution vector is updated based on the transpose of the corrected system response matrix, the relaxation factor parameter, and the diagonally loaded corrected residual vector to generate a new iterative solution vector. In each iteration, the total variation value of the image corresponding to the new iterative solution vector after diagonal loading correction is recorded. The total variation value of the image is obtained by calculating the sum of the absolute values ​​of the differences between the new iterative solution vector in the horizontal and vertical directions. The total variation value of the image is compared with a preset total variation constraint threshold. If the total variation value of the image is greater than the preset total variation constraint threshold, soft threshold shrinkage processing is performed on the new iterative solution vector to generate an iterative solution vector that has passed the total variation constraint. The iterative solution vector that has passed the total variation constraint is used as the current iterative solution vector for the next iteration.

8. The remote sensing image adaptive processing method combining optimization algorithm and probabilistic model according to claim 3, characterized in that, After constructing the Markov random field neighborhood system based on the spatial adjacency relationship between each pixel in the initial density distribution image and its surrounding neighboring pixels, the method further includes: Weighting is performed on the set of neighboring pixels of each pixel in the Markov random field neighborhood system. The spatial distance and grayscale difference between all pixels in the initial density distribution image are normalized. A spatial distance weight coefficient is calculated based on the normalized spatial distance between each neighboring pixel and the center pixel. The spatial distance weight coefficient is inversely proportional to the square of the normalized spatial distance. A grayscale difference weight coefficient is calculated based on the normalized grayscale difference between each neighboring pixel and the center pixel. The grayscale difference weight coefficient is obtained by performing a negative exponential transformation on the absolute value of the normalized grayscale difference. The spatial distance weight coefficient and the grayscale difference weight coefficient are multiplied by a dot product to obtain the dimensionless comprehensive influence weight of each neighboring pixel relative to the center pixel. The neighborhood pixel set of each pixel in the Markov random field neighborhood system is weighted according to the comprehensive influence weight to construct a weighted Markov random field neighborhood system. Based on the weighted Markov random field neighborhood system and Gibbs distribution, a weighted neighborhood correlation constraint energy function is constructed. The weighted neighborhood correlation constraint energy function is obtained by multiplying the difference between the label vector of each neighboring pixel and the label vector of the center pixel by the corresponding comprehensive influence weight and then summing them up. According to the weighted neighborhood correlation constraint energy function, a corresponding weighted neighborhood dependency weight is generated for each pixel. The weighted neighborhood dependency weight is obtained by performing a negative exponential transformation on the weighted neighborhood correlation constraint energy function and then normalizing it. The weighted neighborhood dependency weights are used to replace the neighborhood dependency weights and are input into the iterative segmentation framework based on the coupling of Gaussian mixture model and Markov random field for subsequent pixel category classification processing.

9. The remote sensing image adaptive processing method combining optimization algorithm and probabilistic model according to claim 4, characterized in that, Before performing the maximization step of the expectation-maximization algorithm based on the attribution probability vector corresponding to each pixel, the method further includes: Spatial smoothing filtering is performed on the attribution probability vector corresponding to each pixel. By weighting and averaging the attribution probability vector of each pixel with the attribution probability vectors of all pixels in its neighboring pixel set, a spatial smoothed attribution probability vector corresponding to each pixel is generated. The weights used in the weighted average are determined based on the spatial distance between pixels and the difference in grayscale values. The spatially smoothed attribution probability vector is used as the attribution probability vector in the subsequent maximization step calculation; Based on the probability value of each pixel in the spatial smoothing assignment probability vector belonging to each Gaussian component, the probability values ​​of all pixels belonging to that Gaussian component in the spatial smoothing assignment probability vector corresponding to each Gaussian component are summed to obtain the total spatial smoothing assignment weight of each Gaussian component. Based on the total spatial smoothing assignment weight of each Gaussian component and the grayscale values ​​of all pixels belonging to that Gaussian component, calculate the new spatial smoothing expectation parameter for each Gaussian component. Based on the new spatial smoothing expectation parameter of each Gaussian component, the gray values ​​of all pixels belonging to that Gaussian component, and the corresponding spatial smoothing assignment probability value, calculate the new spatial smoothing variance parameter of each Gaussian component. Calculate the new prior probability parameter of spatial smoothing for each Gaussian component based on the total spatial smoothing weight of each Gaussian component and the total number of pixels in the superpixel segmentation result image. Based on the new spatial smoothing expectation parameter, new spatial smoothing variance parameter, and new spatial smoothing prior probability parameter of each Gaussian component, construct the spatially smoothed intermediate Gaussian mixture model. The intermediate Gaussian mixture model is replaced by the spatially smoothed intermediate Gaussian mixture model for subsequent comparison of parameter changes and convergence determination.

10. A remote sensing image adaptive processing system combining optimization algorithms and probabilistic models, characterized in that, include: processor; A machine-readable storage medium for storing machine-executable instructions of the processor; The processor is configured to execute the remote sensing image adaptive processing method combining optimization algorithms and probabilistic models as described in any one of claims 1 to 9 by executing the machine-executable instructions.

Images