Electromagnetic inverse scattering robust imaging method based on maximum correlation entropy and subspace optimization fusion

Abstract

The invention discloses an electromagnetic inverse scattering robust imaging method based on maximum correlation entropy and subspace optimization fusion. According to the method, a contrast source is decomposed into deterministic / ambiguous subspaces under an SOM framework, a data-state dual-channel robust cost function is constructed through MCC, and non-convex optimization is converted into a sequence weighted quadratic subproblem based on HQ and IRLS. And adaptively suppressing an outlier through Gaussian kernel weight, and updating a closed weighted solution of an ambiguity subspace coefficient and a contrast ratio by adopting an alternating strategy to realize effective suppression of pulse, heavy tail and speckle noise. The method is free in training and physically interpretable, has fast convergence and high robustness, and is suitable for high-fidelity inversion in non-Gaussian noise environments such as radar imaging, medical imaging and underground detection.

Description

Technical Field

[0001] This invention belongs to the field of electromagnetic wave imaging and signal processing, specifically relating to an imaging reconstruction method for electromagnetic inverse scattering problems (EISPs), and in particular a reconstruction method with high robustness and high accuracy in non-Gaussian noise environments. Background Technology

[0002] Electromagnetic backscattering (ESP) technology aims to deduce the spatial distribution of electromagnetic properties (such as dielectric constant and conductivity) within a target by measuring the electromagnetic field data scattered by the target. This technology has significant applications in radar imaging, medical diagnostics, geophysical exploration, and non-destructive testing.

[0003] Its theoretical foundation lies in the Lippmann–Schwinger type integral equations expressed using Maxwell's equations and Green's function. In the continuous domain, this model manifests as the multiplicative coupling of medium contrast and the total field, along with the nonlinear superposition of multiple scattering terms. After discretization, it corresponds to a compact high-dimensional operator mapping, with rapid decay of the singular value spectrum. This makes the separation of the signal and noise subspaces unrobust, and any minute observational perturbation can be amplified in the inverse mapping, resulting in typical ill-posedness. Engineering systems commonly suffer from finite apertures, limited frequency bandwidth, and non-uniform array layouts, further truncating the angular spectrum and spatial frequency domains. This weakens the observability of weak and non-radiative modes, leading to problems such as "invisible edges," resolution degradation, and orientation-dependent artifacts. Simultaneously, model mismatch caused by simplifications in the forward model (e.g., boundary condition approximations, grid discretization errors, Green's function truncation) is iteratively amplified along the state constraint chain, becoming a hidden source affecting convergence stability and physical consistency.

[0004] Most existing mainstream reconstruction algorithms, such as source inversion, subspace optimization methods, and their variants, are based on the minimum mean squared error (MSE) as the optimization objective. The MSE criterion is mathematically equivalent to assuming that noise follows a Gaussian distribution, and it is sensitive to the second moment of the noise. When there are outliers with high energy in the data (generated by non-Gaussian noise), these outliers will dominate the optimization process, causing the algorithm to misinterpret noise as a real scattering phenomenon, thus generating an incorrect physical model. Ultimately, this leads to severe distortion and artifacts in the reconstruction results, and may even prevent convergence to a physically meaningful solution.

[0005] Therefore, there is an urgent need for a novel electromagnetic inverse scattering method that can effectively suppress non-Gaussian noise and improve the stability and accuracy of imaging in complex electromagnetic environments. Summary of the Invention

[0006] This invention proposes a robust inverse scattering imaging method (MCC-SOM) based on "maximum correlation entropy-subspace optimization fusion." By introducing information-theoretic MCC into the SOM structure, a correlation entropy cost is constructed for the data-state dual-channel architecture. A semi-quadratic (HQ) and iterative reweighted least squares (IRLS) framework is employed, along with alternating updates of deterministic / ambiguous subspaces, to adaptively suppress non-Gaussian noise such as impulses, heavy tails, and speckle, while maintaining efficient convergence and physical interpretability. This invention achieves high-fidelity reconstruction in various non-Gaussian noise environments without requiring training data or complex network structures. It possesses real-time / near-real-time potential and good interpretability, making it suitable for widespread application in fields with high safety and reliability requirements, such as medicine, electromagnetic non-destructive testing, and underground target detection.

[0007] The technical solution of this invention is a robust electromagnetic inverse scattering imaging method based on the fusion of maximum correlation entropy and subspace optimization. This method includes:

[0008] Step 1: Establish a discretized operator model for the electromagnetic scattering problem;

[0009] The target region D to be inverted is discretized into N grid cells. Based on this, the two fundamental integral equations describing the scattering phenomenon, the state equation and the data equation, are discretized to obtain the following matrix-form operator model:

[0010] Discrete state equations: describe the total field within the scattering region D. From the incident field and comparison source The resulting fields are superimposed, as shown in Equation 1:

[0011] (1);

[0012] Discrete data equation: This equation describes the scattered field measured at the receiving antenna outside the scattering region. It is a contrast source within region D. The radiation produced is shown in Equation 2:

[0013] (2);

[0014] The symbols are defined as follows:

[0015] : The known quantity is the vectorized representation of the incident field over N discrete grid cells;

[0016] : The known quantity is the vectorized representation of the scattered field data measured by Nr receiving antennas;

[0017] The unknown quantity is the vectorized representation of the comparison source to be solved on N discrete grid cells;

[0018] Unknowns are the vectorized representations of the target contrast function to be solved in N discrete grid cells; contrast sources. With contrast Satisfying Relationship:

[0019] (3);

[0020] in Represents the vector The diagonal matrix formed;

[0021] The discrete-domain Green's operator is a pre-computed matrix based on a discretized grid, used to describe the interactions between discrete elements within a region D, compared to the source. This is mapped to the scattered field generated within region D;

[0022] The discrete scattering Green's operator is a pre-calculated matrix based on a discretized grid and the location of the receiving antenna. It describes the field propagation from the source to the external receiving antenna within region D, and is compared with the source... This is mapped to the scattered field generated at the receiving antenna;

[0023] Step 2: Subspace decomposition and initialization of deterministic / ambiguous components based on SVD;

[0024] Green's operator for discrete scattering Perform singular value decomposition:

[0025] (4);

[0026] in and It is a unitary matrix, and its column vectors and These are the left singular vector and the right singular vector; diagonal matrix Includes real non-negative singular values ​​arranged in non-increasing order. Right singular vector set Constructing the contrast source space A complete orthogonal basis, Depend on composition, Depend on composition;

[0027] With a selected cutoff threshold L, the basis is divided into two positively complementary subspaces: a deterministic subspace and a non-deterministic subspace. Ambiguous subspaces The deterministic subspace consists of the first L right singular vectors. Zhang Cheng, these vectors correspond to singular values ​​above a certain noise correlation threshold, representing components in the contrast source with high radiative intensity that can be well captured by the receiver; the ambiguous subspace consists of the remaining NL right singular vectors. Zhang Cheng represents weak radiation or non-radiative source components;

[0028] For the incident field of the i-th transmitter The corresponding comparison source Projected onto these two subspaces, the deterministic part is obtained. Ambiguous parts ;in By applying scattering field data The truncated pseudoinverse is obtained by non-iterative solution, as shown in Equation 5; while the ambiguous part Located in the noise subspace, it cannot be determined by data equations alone. It is a linear combination of the basis vectors of the noise subspace, as shown in Equation 6.

[0029] (5);

[0030] (6);

[0031] in The vector of coefficients to be estimated;

[0032] Step 3: Construct a dual-channel cost function based on the maximum correlation entropy criterion;

[0033] The electromagnetic inverse scattering inversion problem described above can be simplified back to finding the contrast distribution. sum of coefficient vector sets Rephrasing it as an optimization problem, its cost function is... J MCC It consists of data fidelity and state fidelity, both of which are measured based on the maximum relevance entropy criterion;

[0034] Step 4: Perform an equivalent transformation on the cost function and establish an iterative weighted least squares framework;

[0035] Step 5: Iteratively optimize the solution using alternating steps;

[0036] In the (k+1)th iteration, an alternating optimization strategy is adopted, fixing one variable and solving for the optimal solution of the other variable, and updating the weights and fuzzy source coefficients sequentially. and contrast distribution ;

[0037] Step 6: Convergence judgment and result output;

[0038] Repeat the alternating iteration process in step 5 until the preset convergence condition is met; after the iteration terminates, output the final contrast function. This is the result of reconstructing the spatial distribution of the electromagnetic parameters of the target.

[0039] Furthermore, the specific method for step 3 is as follows:

[0040] Define the data residual corresponding to the i-th incident wave as given by the current comparison source. Calculated scattering field Compared with the actual measured scattering field The difference between them is shown in Equation 7; the state residual corresponding to the i-th incident wave is defined by the current comparison source. With current contrast The difference between the contrast source calculated from the incident field and the contrast source is shown in Equation 8;

[0041] (7);

[0042] (8);

[0043] For the residual vector The cost function based on the relevant entropy uses a Gaussian kernel. Build, designed to maximize This is mathematically equivalent to minimizing the following objective function:

[0044] (9);

[0045] By applying a Gaussian kernel to these residuals and summing over all emitters and spatial locations, the objective function is:

[0046] (10);

[0047] Where, N i This represents the total number of transmitting antennas. >0 is a regularization parameter used to balance data and state fidelity terms; while >0 is the kernel width, used to control the robustness of the estimator: smaller values ​​produce stronger outlier removal; larger values ​​are closer to the minimum mean squared error.

[0048] Furthermore, the specific method for step 4 is as follows:

[0049] For each residual term, a weight is defined as shown in Equations 11 and 12;

[0050] (11);

[0051] (12);

[0052] minimize This is equivalent to minimizing the following surrogate objective function, which is a quadratic function with respect to its variables, while keeping other variables constant:

[0053] (13);

[0054] in and It is a diagonal weight matrix; This represents the two weights established.

[0055] Furthermore, the specific method of step 5 is as follows: Step 5.1: Update weights The superscript k indicates the corresponding iteration number: using the previous iteration's... and Calculate residuals and As shown below;

[0056] (14);

[0057] (15);

[0058] Then calculate the weight matrix for all incident waves. and ;

[0059] Step 5.2: Update the fuzzy source coefficients :fixed and weight , minimize function ;set up , and The weighted form of the i-th incident wave is:

[0060] (16);

[0061] Therefore, we obtain The system of linear equations:

[0062] (17);

[0063] Equation (17) is solved using the conjugate gradient method;

[0064] Step 5.3: Update the contrast function Fixed fuzzy source coefficients and weight Total contrast source and the main field Calculated using the following formula:

[0065] (18);

[0066] (19);

[0067] For each spatial unit n, regarding objective function Minimization is decoupled; therefore, each element Closed-form solution:

[0068] (20).

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

[0070] 1) Excellent robustness: By introducing the maximum correlation entropy criterion (MCC), this method can adaptively identify and suppress outliers caused by non-Gaussian noise (especially impulse noise and heavy-tailed noise) in the measurement data. Its influence function is bounded and reducible, fundamentally avoiding the excessive influence of outliers on the reconstruction results, making the algorithm extremely stable in complex electromagnetic environments.

[0071] 2) Higher reconstruction accuracy: Even in noise-free or low-noise environments, MCC can capture higher-order statistical properties of data than MSE, thus modeling the scattering process more accurately. Experiments show that this method significantly outperforms traditional SOM, TSOM, and KSOM algorithms in objective metrics such as structural similarity (SSIM) and peak signal-to-noise ratio (PSNR), and can reconstruct target images with sharper edges and finer structures.

[0072] 3) High computational efficiency: This invention retains the basic framework of subspace optimization methods, restricting the main optimization process to a low-dimensional fuzzy subspace, and transforming the non-convex problem into an efficient iterative reweighted least squares problem through a semi-quadratic framework. Each subproblem has an analytical solution or an efficient numerical solution, ensuring that the algorithm maintains high computational efficiency while improving robustness.

[0073] 4) Strong physical interpretability: Unlike deep learning methods that rely on large amounts of training data, this invention is an optimization method driven by a physical model. Its robustness stems from the introduction of information theory principles. The entire optimization process is transparent and interpretable, does not depend on the training dataset, and has better generalization ability and practical application potential. Attached Figure Description

[0074] Figure 1Schematic diagrams of irregular geometric shapes measured in the experiment (Austria figure, overlapping circle, MNIST handwritten font, and overlapping ring).

[0075] Figure 2 This is a schematic diagram of the experimental results for the backscattering imaging of the Austrian pattern.

[0076] Figure 3 The graphs represent the iterative curves of the Austrian graph under various algorithms.

[0077] Figure 4 This is a schematic diagram of the experimental results of backscattering imaging of overlapping circular images.

[0078] Figure 5 The graph shows the iterative curves of the overlapping circle image under various algorithms.

[0079] Figure 6 This is a schematic diagram of the experimental results of backscattering imaging of handwritten characters in MNIST.

[0080] Figure 7 The iterative curves of MNIST handwritten fonts under various algorithms.

[0081] Figure 8 This is a schematic diagram of the experimental results of backscattering imaging of overlapping ring images.

[0082] Figure 9 The graph shows the iterative curves of the overlapping ring image under various algorithms. Detailed Implementation

[0083] To make the objectives, technical solutions, and advantages of this invention clearer, a specific embodiment will be used to describe the technical solution of this invention in detail below. It should be understood that the specific embodiment described herein is only for explaining this invention and is not intended to limit the scope of protection of this invention.

[0084] This embodiment uses a multi-station static radar system under two-dimensional transverse magnetic waves as an example to demonstrate the process of accurate target imaging by the method of the present invention in a strong impulse noise environment.

[0085] I. Experimental Environment and Parameter Settings

[0086] 1. Physical Model: A two-dimensional scattering scene is set, with free space as the background. The target to be imaged is located within a 2m × 2m square detection area D. The operating frequency of the electromagnetic waves is set to 300MHz.

[0087] 2. Detection System: A multi-station static radar configuration is employed. N radars are uniformly positioned around the detection area D. t =16 transmitting antennas and N r=32 receiving antennas. During operation, 16 transmitting antennas sequentially emit TM-polarized cylindrical waves to illuminate the detection area, and all 32 receiving antennas synchronously receive the scattered field data.

[0088] 3. Imaging Target: The grayscale image of the digits extracted from the irregular geometric image is used as the true contrast distribution of the target to be imaged. The grayscale values ​​of the image are mapped to the contrast values ​​of the medium; for example, the pixel value range [0,1] is mapped to the contrast range [0,1.5].

[0089] 4. Discretization: The 2m×2m detection area D is discretized into a uniform grid of 80×80, with a total of N=6400 discrete elements.

[0090] 5. Operator Generation: Based on the discretized grid and antenna layout described above, the discrete-domain Green's operator is calculated using the method of moments. and discrete scattering Green's operator These two matrices are fixed and are calculated and stored all at once before the inversion begins. G D The dimensions are 6400×6400, G S The dimensions are (32×4096).

[0091] 6. Simulation data generation and noise injection:

[0092] First, the scattered field data generated by the target under ideal noise-free conditions are calculated using a forward scattering model. .

[0093] Then, in order to simulate harsh measurement environments, [the following was done]... Inject strong noise. Specifically:

[0094] 1) To simulate the granular artifacts unique to coherent imaging systems, a multiplicative noise model is used. The expression is as follows:

[0095] (12)

[0096] in This represents zero-mean complex Gaussian white noise. It is a scalar quantity for controlling noise intensity.

[0097] 2) To explain the multipath propagation effect and other non-Gaussian interference sources that produce significant outliers, additive noise from the Student t-distribution is introduced.

[0098] (13)

[0099] in For complex noise, the real and imaginary parts are independently sampled with degrees of freedom of . The t-distribution. Smaller. This will result in a heavy-tailed distribution, and It then acts as a scaling factor.

[0100] 3) To represent sensor malfunctions or data transmission errors, we simulate impulse noise. Specifically, we randomly select 3% of the data points and forcibly replace their amplitudes with an abnormally large value (scaled to 10 times the original amplitude in this example), thereby generating simulated measurement data with strong outliers, as shown in the following expression.

[0101]

[0102] in , The original data was magnified 10 times.

[0103] The three noise scenarios mentioned above are recorded as Noise1, Noise2, and Noise3, respectively.

[0104] Detailed experimental results are attached. Figures 1 to 9 The quantitative data of the experimental results are shown in Tables 1 and 2. The results show that, in different datasets, the method proposed in this invention has higher output accuracy and better stability compared with the classical SOM method, the dual subspace method (TSOM), and the kernel subspace optimization method (KSOM).

[0105] The entire pseudocode process of this invention is shown in the table below.

[0106]

[0107] The innovation of this invention lies in:

[0108] 1. This invention is the first to introduce the "Maximum Correlation Entropy Criterion" (MCC) from information theory into the field of electromagnetic inverse scattering, replacing the "Minimum Mean Square Error" (MSE) criterion relied upon by traditional methods. MSE is extremely sensitive to noise, while MCC, using the kernel function method, can fundamentally "passivate" the interference of non-Gaussian noise such as strong impulses and heavy tails. This is the theoretical foundation for the robustness of this invention.

[0109] 2. This invention designs a highly efficient iterative adaptive weighting algorithm. In each calculation, it automatically identifies "bad data" contaminated by noise and significantly reduces its weight, while amplifying the influence of "good data." This ensures that even when the data is severely contaminated, the inversion process converges in the correct direction.

[0110] 3. This invention solves the industry problem of traditional methods failing completely in non-Gaussian noise environments, achieving stable and high-precision robust imaging. More importantly, while achieving strong robustness, it does not sacrifice performance in noise-free or conventional Gaussian noise environments, achieving a "full-scene" performance improvement, and has extremely high practical value and wide applicability.

[0111] The above description of specific embodiments details the implementation process and application scenarios of the maximum correlation entropy comparison source inversion method, enabling those skilled in the art to implement the maximum correlation entropy comparison source inversion method according to the present invention. The above embodiments are only used to illustrate the technical solution of the present invention and do not limit the scope of protection of the present invention. Those skilled in the art can make various modifications and variations based on the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

[0112] Table 1. Comprehensive evaluation results of SSIM and PSNR under different noise conditions;

[0113] Table 2. Comparison of the improvement effects of SSIM and PSNR under different noise conditions;

[0114]

Claims

1. A robust electromagnetic inverse scattering imaging method based on the fusion of maximum correlation entropy and subspace optimization, the method comprising: Step 1: Establish a discretized operator model for the electromagnetic scattering problem; The target region D to be inverted is discretized into N grid cells. Based on this, the two fundamental integral equations describing the scattering phenomenon, the state equation and the data equation, are discretized to obtain the following matrix-form operator model: Discrete state equations: describe the total field within the scattering region D. From the incident field and comparison source The resulting fields are superimposed, as shown in Equation 1: (1)； Discrete data equation: This equation describes the scattered field measured at the receiving antenna outside the scattering region. It is a contrast source within region D. The radiation produced is shown in Equation 2: (2)； The symbols are defined as follows: : The known quantity is the vectorized representation of the incident field over N discrete grid cells; : The known quantity is the vectorized representation of the scattered field data measured by Nr receiving antennas; The unknown quantity is the vectorized representation of the comparison source to be solved on N discrete grid cells; Unknowns are the vectorized representations of the target contrast function to be solved in N discrete grid cells; contrast sources. With contrast Satisfying Relationship: (3)； in Represents the vector The diagonal matrix formed; The discrete-domain Green's operator is a pre-computed matrix based on a discretized grid, used to describe the interactions between discrete elements within a region D, compared to the source. This is mapped to the scattered field generated within region D; The discrete scattering Green's operator is a pre-calculated matrix based on a discretized grid and the location of the receiving antenna. It describes the field propagation from the source to the external receiving antenna within region D, and is compared with the source... This is mapped to the scattered field generated at the receiving antenna; Step 2: Subspace decomposition and initialization of deterministic / ambiguous components based on SVD; Green's operator for discrete scattering Perform singular value decomposition: (4)； in and It is a unitary matrix, and its column vectors and These are the left singular vector and the right singular vector; diagonal matrix Includes real non-negative singular values ​​arranged in non-increasing order. Right singular vector set Constructing the contrast source space A complete orthogonal basis, Depend on composition, Depend on composition; With a selected cutoff threshold L, the basis is divided into two positively complementary subspaces: a deterministic subspace and a non-deterministic subspace. Ambiguous subspaces The deterministic subspace consists of the first L right singular vectors. Zhang Cheng, these vectors correspond to singular values ​​above a certain noise correlation threshold, representing components in the contrast source with high radiative intensity that can be well captured by the receiver; the ambiguous subspace consists of the remaining NL right singular vectors. Zhang Cheng represents weak radiation or non-radiative source components; For the incident field of the i-th transmitter The corresponding comparison source Projected onto these two subspaces, the deterministic part is obtained. Ambiguous parts ;in By applying scattering field data The truncated pseudoinverse is obtained by non-iterative solution, as shown in Equation 5; while the ambiguous part Located in the noise subspace, it cannot be determined by data equations alone. It is a linear combination of the basis vectors of the noise subspace, as shown in Equation 6. (5)； (6)； in The vector of coefficients to be estimated; Step 3: Construct a dual-channel cost function based on the maximum correlation entropy criterion; The electromagnetic inverse scattering inversion problem described above can be simplified back to finding the contrast distribution. sum of coefficient vector sets Rephrasing it as an optimization problem, its cost function is... J MCC It consists of data fidelity and state fidelity, both of which are measured based on the maximum relevance entropy criterion; Step 4: Perform an equivalent transformation on the cost function and establish an iterative weighted least squares framework; Step 5: Iteratively optimize the solution using alternating steps; In the (k+1)th iteration, an alternating optimization strategy is adopted, fixing one variable and solving for the optimal solution of the other variable, and updating the weights and fuzzy source coefficients sequentially. and contrast distribution ; Step 6: Convergence judgment and result output; Repeat the alternating iteration process in step 5 until the preset convergence condition is met; after the iteration terminates, output the final contrast function. This is the result of reconstructing the spatial distribution of the electromagnetic parameters of the target.

2. The robust electromagnetic inverse scattering imaging method based on the fusion of maximum correlation entropy and subspace optimization as described in claim 1, characterized in that, The specific method for step 3 is as follows: Define the data residual corresponding to the i-th incident wave as given by the current comparison source. Calculated scattering field Compared with the actual measured scattering field The difference between them is shown in Equation 7; Define the state residual corresponding to the i-th incident wave as given by the current comparison source. With current contrast The difference between the contrast source calculated from the incident field and the contrast source is shown in Equation 8; (7)； (8)； For the residual vector The cost function based on the relevant entropy uses a Gaussian kernel. Build, designed to maximize This is mathematically equivalent to minimizing the following objective function: (9)； By applying a Gaussian kernel to these residuals and summing over all emitters and spatial locations, the objective function is: (10)； Where, N i This represents the total number of transmitting antennas. >0 is a regularization parameter used to balance data and state fidelity terms; while >0 is the kernel width, used to control the robustness of the estimator: smaller values ​​produce stronger outlier removal; larger values ​​are closer to the minimum mean squared error.

3. The robust electromagnetic inverse scattering imaging method based on the fusion of maximum correlation entropy and subspace optimization as described in claim 1, characterized in that, The specific method for step 4 is as follows: For each residual term, a weight is defined as shown in Equations 11 and 12; (11)； (12)； minimize This is equivalent to minimizing the following surrogate objective function, which is a quadratic function with respect to its variables, while keeping other variables constant: (13)； in and It is a diagonal weight matrix; This represents the two weights established.

4. The robust electromagnetic inverse scattering imaging method based on the fusion of maximum correlation entropy and subspace optimization as described in claim 1, characterized in that, The specific method for step 5 is as follows: Step 5.1: Update weights The superscript k indicates the corresponding iteration number: using the previous iteration's... and Calculate residuals and As shown below; (14)； (15)； Then calculate the weight matrix for all incident waves. and ; Step 5.2: Update the fuzzy source coefficients :fixed and weight , minimize function ;set up , and The weighted form of the i-th incident wave is: (16)； Therefore, we obtain The system of linear equations: (17)； Equation (17) is solved using the conjugate gradient method; Step 5.3: Update the contrast function Fixed fuzzy source coefficients and weight Total contrast source and the main field Calculated using the following formula: (18)； (19)； For each spatial unit n, regarding objective function Minimization is decoupled; therefore, each element Closed-form solution: (20)。

Images