IC packaging fault positioning method based on fixed interval multilayer current distribution reconstruction

Abstract

The invention discloses an IC packaging fault positioning method and system based on fixed interval multilayer current distribution reconstruction, and relates to the technical field of integrated circuit testing. The method comprises the following steps: acquiring magnetic field measurement data of an IC package surface; establishing a transfer function model in a geometric progression form according to the fixed spacing of the conductive layers in the package; performing inversion on magnetic field data by using the model, and reconstructing three-dimensional current distribution of each layer; and realizing accurate fault positioning based on abnormal current distribution. According to the method, the prior physical constraint of the fixed spacing is converted into a mathematical structured optimization problem, the calculation complexity is remarkably reduced, the reconstruction time is shortened to several seconds from several minutes while millimeter-level high-precision positioning is realized, the method has excellent anti-noise capability, and the method is suitable for large-scale popularization and application. And an effective solution is provided for rapid and nondestructive fault diagnosis of integrated circuit packaging.

Description

Technical Field

[0001] This invention relates to integrated circuit testing technology, and in particular to integrated circuit (IC) package fault location technology based on multilayer current distribution reconstruction. Background Technology

[0002] With the rapid development of integrated circuit technology, chip integration is constantly increasing, and multi-layer packaging structures are becoming increasingly common. During IC packaging manufacturing and use, faults such as short circuits and open circuits may occur in the internal conductive layers, seriously affecting circuit performance and reliability. Traditional fault detection methods, such as X-ray inspection and acoustic microscopy, have problems such as limited resolution, long inspection time, and high cost.

[0003] Magnetic field imaging technology, as a non-contact and non-destructive detection method, has shown great potential in the field of IC fault detection. However, existing magnetic field imaging technologies mainly target single-layer or two-dimensional current distributions, and face the following technical challenges in three-dimensional current reconstruction of multi-layer package structures:

[0004] 1. Ill-conditioned inverse problem: Reconstructing the internal multilayer current distribution from the surface magnetic field is a typical ill-conditioned inverse problem, and the uniqueness and stability of the solution are difficult to guarantee.

[0005] 2. Interlayer interference: The magnetic fields generated by the currents in multiple layers are superimposed, making it difficult to accurately separate the contributions of each layer.

[0006] 3. Depth attenuation: The magnetic field decays exponentially with distance. The magnetic field signal generated by deep current is weak and has a low signal-to-noise ratio.

[0007] 4. Computational complexity: Traditional full 3D reconstruction methods have high computational complexity, making it difficult to achieve rapid detection.

[0008] In existing technologies, the Fourier method proposed by Roth et al. is only applicable to two-dimensional current reconstruction; the Bayesian method can handle noise, but requires a large amount of prior information; machine learning methods require a large amount of training data and have limited generalization ability. These methods are all difficult to effectively utilize the structured features of IC packaging.

[0009] While existing research has attempted to reconstruct the current distribution of packaged structures in full three dimensions to address multilayer problems, these methods typically rely on computationally intensive full-wave electromagnetic simulations or require solving large-scale inverse problems. Their extremely high computational complexity and massive resource requirements result in extremely long reconstruction times, failing to meet the demands for rapid location in practical fault analysis. This leads to: 1) difficulty in completing diagnosis within minutes or even seconds, thus losing its application value for online or rapid sampling inspections; 2) extreme sensitivity to noise, with reconstruction results easily becoming unstable under limited time and data, resulting in a significant decrease in location accuracy (especially in the depth direction); and 3) poor method universality, with models and computational processes difficult to adapt to different package specifications. Summary of the Invention

[0010] This invention addresses the problem of excessive computational complexity in current three-dimensional current reconstruction by providing a method that fully utilizes the deterministic prior structural knowledge that conductive layers in IC packages are stacked at fixed intervals during reconstruction. It also provides a method that can integrate the strong physical constraint of fixed intervals, thereby significantly reducing computational complexity and achieving fast, stable, and accurate three-dimensional fault location.

[0011] The technical solution adopted by this invention to solve the above-mentioned technical problems is an IC package fault location method based on fixed-gap multilayer current distribution reconstruction, comprising the following steps:

[0012] By utilizing the geometric feature of the fixed spacing of conductive layers in IC packages, the internal current distribution of multiple layers can be reconstructed from surface magnetic field measurement data to achieve fault location. This includes the following steps:

[0013] Steps for acquiring magnetic field data: Acquire magnetic field measurement data at multiple probe points on the surface of the IC package under test;

[0014] Modeling steps for geometric series transfer functions: Based on the number of conductive layers N and the fixed interlayer spacing d obtained from the IC package specifications, determine the depth of each layer. i is the layer number; establish a system based on the fixed layer spacing d and depth. The geometric series transfer function G is a parameter that characterizes the frequency domain mapping relationship between the single-layer current density and the generated magnetic field, and the ratio of the transfer functions of adjacent layers is constant.

[0015] Current distribution reconstruction steps: The magnetic field measurement data is inverted and calculated using the geometric series transfer function G to reconstruct the three-dimensional current density distribution of each conductive layer inside the IC package. ;

[0016] Fault location steps: Based on the reconstructed three-dimensional current density distribution of each layer Identify areas of abnormal current distribution and determine their depth within the layer. Three-dimensional fault location is achieved using planar coordinates.

[0017] Specifically, the geometric series transfer function G in the beam domain The expression is:

[0018] ;

[0019] in, These are the x and y components of the spatial frequency k, respectively. Where is the free permeability and radial wavenumber. , It is a natural constant.

[0020] Preferably, the current distribution reconstruction step includes:

[0021] Preliminary solution of frequency domain decomposition: The magnetic field measurement data is converted to the wavenumber domain. Using the analytical properties of the geometric series transfer function model G, the initial separation and reconstruction of the current contribution of each layer is achieved through recursive filtering, and the initial current density of each layer is obtained.

[0022] Iterative refinement solution: Starting from the initial current density, the current distribution of each layer is iteratively refined by minimizing the objective function L, which includes data fitting and regularization terms, to obtain the refined current density of each layer.

[0023] Global optimization: Starting with the refined current density of each layer, a multi-objective optimization problem is constructed and solved, which simultaneously considers the current distribution of all layers, data fitting terms, and regularization constraints, to obtain the globally optimal current density distribution.

[0024] This invention is based on the observed complex three-dimensional electromagnetic coupling effect introduced by the strict stacking of multiple conductive layers at a fixed spacing (determined by the thickness of the dielectric layer). When faults such as microbump cracking, TSV failure, interlayer short circuits, or high resistance occur inside the package, the current disturbance generated at the fault point will be transmitted to both the lateral and longitudinal directions in a specific mode through this fixed interlayer structure. By making full use of the prior knowledge that the conductive layer spacing in IC packaging is fixed, the complex three-dimensional reconstruction problem is transformed into a structured optimization problem.

[0025] The beneficial effects of this invention are that, by fully utilizing the structured characteristics of IC packaging, it achieves a breakthrough in the high-precision reconstruction of multilayer internal current distribution from surface magnetic field measurement data. This method cleverly transforms the physical constraint of fixed conductive layer spacing in IC packaging into a mathematical geometric series form, turning the originally ill-conditioned three-dimensional current reconstruction problem into a structured optimization problem, greatly improving the stability and accuracy of the solution. In practical applications, the reconstruction accuracy of this invention can reach the 0.05 mm level, which is 3-5 times higher than the traditional Fourier transform method. This high precision enables accurate identification of sub-circuit-level connection faults, which is of great significance for fault diagnosis in modern high-density packaging. Simultaneously, by employing geometric series decomposition technology, the computational complexity of traditional full three-dimensional reconstruction is reduced from O(N³) to O(N²logN), and the reconstruction time for a typical four-layer PCB is shortened from several minutes using traditional methods to 3-5 seconds, enabling real-time online detection. Another important advantage lies in its excellent robustness and adaptability. By introducing physically based regularization constraints, including current continuity constraints, exponential decay constraints, and interlayer coupling constraints, this method maintains stable reconstruction performance even under harsh conditions with signal-to-noise ratios as low as 10 dB, which is crucial for applications in practical industrial environments. The regularization parameters can adaptively adjust according to layer depth and noise level, eliminating the need for tedious manual parameter tuning and greatly improving the method's practicality. More importantly, this method does not require detailed internal circuit layout information as prior knowledge; it only needs to obtain interlayer spacing parameters from the IC package specification sheet to achieve high-quality reconstruction. This allows the method to be widely applied to various standard package forms, including BGA, QFP, SOP, and CSP, providing an efficient, accurate, and economical solution for quality control and fault analysis in the integrated circuit manufacturing industry. Attached Figure Description

[0026] Figure 1 This is a flowchart of an implementation example;

[0027] Figure 2 A comparison chart of multilayer current reconstruction results;

[0028] Figure 3 This is a diagram demonstrating fault location. Detailed Implementation

[0029] like Figure 1 As shown, an IC package fault location method based on fixed-spacing multilayer current distribution reconstruction is presented. Utilizing the known fixed interlayer spacing, the solution in the depth direction (z-axis) is constrained, thereby quickly and accurately locating the depth (layer) and planar position of the fault. The specific steps include:

[0030] Step 1: Acquisition of magnetic field data

[0031] Local electromagnetic field simulation analysis of IC package structure was performed using the COMSOL Multiphysics and MATLAB co-simulation platform. A three-dimensional finite element model of the IC package was established in COMSOL, including key structures such as conductive layers, dielectric layers, and solder ball arrays, and material property parameters such as conductivity and permeability were set. Maxwell's equations were solved using COMSOL's AC / DC module to calculate the magnetic field distribution under a given current excitation. The magnetic field data calculated by COMSOL was imported into the MATLAB environment via the LiveLink for MATLAB interface. MATLAB scripts controlled various parameters of COMSOL, output data, and then performed data preprocessing and format conversion.

[0032] Step 2: Modeling the geometric series transfer function

[0033] Based on the standardized structural characteristics of IC packaging, a geometric series transfer function model between magnetic field and current is established, transforming the multilayer current reconstruction problem into a structured mathematical framework. First, key parameters are extracted from the IC packaging specifications, including the number of conductive layers N, interlayer spacing d, and substrate thickness, to determine the precise depth of each conductive layer. Where z0 is the distance from the sensor to the top layer, and i is the layer number. Based on the Biot-Savart law and the quasi-static approximation, the current density of the i-th layer is derived. Magnetic field generated in the frequency domain In the beam domain The expression is:

[0034] ;

[0035] in, For transfer functions, These are the x and y components of the spatial frequency k, respectively. Let be the depth of the i-th layer. Transfer function. , The permeability of free space, To take the modulus of k, It exhibits an exponential decay characteristic with depth. Utilizing the key constraint of fixed interlayer spacing, the transfer functions of multiple layers are organized into a geometric series, with the transfer functions of adjacent layers satisfying a proportional relationship. This regularity gives the transfer matrix a Toeplitz structure, which can be inverted using a fast algorithm. In the frequency domain, different spatial frequency components have different attenuation rates; low-frequency components attenuate slowly and penetrate deeply, while high-frequency components attenuate quickly and mainly reflect shallow information. By analyzing the contributions of different frequency components, layer-selective reconstruction can be achieved. The singular value decomposition of the transfer function reveals the ill-conditioned nature of the reconstruction problem; the condition number increases sharply with the number of layers and depth, requiring regularization techniques to stabilize the solution process.

[0036] Step 3: Preliminary solution of frequency domain decomposition

[0037] The magnetic field measurement data is converted to the wavenumber domain. Utilizing the analytical properties of the geometric series transfer function model G, a recursive filtering method is used to initially separate and reconstruct the current contributions of each layer, yielding initial estimates of the current density for each layer. Specifically:

[0038] The measured surface magnetic field data are subjected to a two-dimensional fast Fourier transform. In the frequency domain, the analytical properties of geometric series are used to achieve rapid separation and preliminary reconstruction of the currents in each layer. The measured magnetic field is the total magnetic field, which is represented in the frequency domain as a linear superposition of the contributions from each layer. N is the total number of IC layers. This constitutes a system of linear equations about the current in each layer, whose coefficient matrix is ​​composed of geometric series transfer functions. Using the geometric series summation formula, the multilayer coupling problem is transformed into a recursive form. By constructing a recursive filter, layer-by-layer stripping is achieved from shallow to deep layers. For the i-th layer (i=1,2,...,N), the stripping process is as follows:

[0039] 1. From the residual magnetic field Reconstructing the current density of the i-th layer Initial conditions: The initial residual magnetic field is equal to the measured magnetic field;

[0040] 2. Calculate the magnetic field generated by the i-th layer:

[0041] 3. Update the residual magnetic field: Used for rebuilding the next layer.

[0042] Set the critical wavenumber For each frequency component Different reconstruction strategies are adopted based on their penetration depth characteristics:

[0043] when Therefore, it is assumed that the low-frequency components contain deep information, and the initial current density is reconstructed by joint inversion: by solving the coupled multilayer linear system. Simultaneously, the initial current density of all layers is reconstructed, where For depth z i The transfer function at the location;

[0044] when If the high-frequency components are considered to mainly originate from shallow layers, the initial current density can be directly reconstructed by applying an inverse Fourier transform to each layer independently. Reconstruct the initial current density.

[0045] To improve numerical stability, Tikhonov regularization is introduced, which suppresses noise amplification through singular value truncation or filtering factors. The regularization parameter is adaptively selected based on the generalized cross-validation criterion. The frequency domain decomposition results provide preliminary estimates of the current in each layer. Although inter-layer crosstalk and edge effects may exist, they provide good initial values ​​for subsequent iterative optimization, significantly accelerating the convergence speed. The frequency domain solution is transformed back to the spatial domain by inverse Fourier transform, yielding a two-dimensional distribution image of the current density in each layer.

[0046] Step 4: Iterative refinement solution

[0047] Starting with the initial current density, the current distribution of each layer is iteratively refined by minimizing the objective function L, which includes data fitting and regularization terms, to obtain the refined current density of each layer. Specifically:

[0048] After obtaining the preliminary solution of the frequency domain decomposition, the iterative least squares method is used to refine the solution of the current distribution of each layer in order to improve the reconstruction accuracy and eliminate the errors that may be introduced by the frequency domain decomposition.

[0049] The iterative process obtains the initial current density of each layer from the frequency domain decomposition. Initially, the current distribution is updated in each iteration by calculating the gradient of the objective function. The specific update formula is as follows: Where α is the step size parameter, This indicates the nth iteration. Let L represent the current density of the i-th layer in the nth iteration, which is dynamically adjusted using a line search or adaptive algorithm to ensure convergence. The objective function L includes a data fitting term and a regularization term. Its gradient is obtained through the chain rule. ,in, According to the transfer function Construct the system matrix. This represents the three-dimensional current density distribution of all conductive layers within the IC package. To measure the magnetic field, For regularization parameters, Represents the regularization term gradient, Indicates transpose. For data fitting terms;

[0050] This involves the Jacobian matrix of the magnetic field transfer function and the current residual. To accelerate convergence, efficient optimization algorithms such as the conjugate gradient method or the quasi-Newton method are employed, and a momentum term is introduced to avoid getting trapped in local minima. The iterative process continues until the convergence condition is met. ,in, The L2 norm is used, and ε is a preset relative error threshold, typically set to 10. -6Simultaneously monitor the residual norm. The iteration process is terminated early when the residual no longer decreases significantly, avoiding overfitting. The entire iteration process is optimized through parallel computing, distributing the computation of different frequency components across multiple processor cores, significantly improving computational efficiency.

[0051] Step 5: Global Optimization

[0052] The objective optimization function that comprehensively considers data fitting, physical laws, and structural information is as follows:

[0053] ;

[0054] The optimal current distribution is solved using a global optimization algorithm. In the matrix The middle corresponds to the depth The submatrix. The data fitting term ensures that the reconstructed magnetic field is statistically consistent with the measured data, using a weighted least squares criterion, with weights adaptively set according to the noise level of each measurement point. That is, the regularization term. Regularization parameters Including weighting coefficients for physical constraint terms Weighting coefficients of structural constraint terms The regularization parameters λ1 and λ2 are automatically determined using the L-curve method or generalized cross-validation to ensure an optimal balance between noise suppression and detail preservation.

[0055] Physical constraints It contains four constraints:

[0056] ;

[0057] in, , , and These are the weighting coefficients for the four constraints. The four constraints are as follows:

[0058] Current continuity constraint ;

[0059] Exponential decay constraint ;

[0060] Interlayer coupling constraints ;

[0061] curl constraint ;

[0062] in, It is the attenuation coefficient. Related to the material's conductivity and frequency, and It refers to the current distribution on the adjacent upper and lower surfaces of two adjacent layers. for of Quantity, for of Quantity, express of Component pairs The partial derivative of .

[0063] Structural constraints Integrating IC layout information, including geometric features such as conductor width, spacing, and corners, and implementing it through total variation, sparsity, smoothness, and wavelet domain constraints:

[0064] ;

[0065] in, , , and These are the weighting coefficients for the four constraints. The four constraints are as follows:

[0066] Total variational constraints ;

[0067] sparsity constraints ;in The z-component represents the current density in the i-th layer. Although the current mainly flows in the x and y planes of the i-th layer, due to the current continuity equation... The constraint necessarily implies the existence of a z-direction component perpendicular to the plane, which can be achieved through... Derive from in-plane components;

[0068] Smoothness constraint ;

[0069] Wavelet domain constraints ;

[0070] in, Let L represent the weighting coefficients of the high-frequency components of the l-th wavelet layer, and L represent the total number of wavelet decomposition layers. This represents the high-pass filter operator in the x-direction when performing wavelet decomposition at the l-th level. This represents the y-direction high-pass filter operator when performing wavelet decomposition at the l-th level. The operator representing the approximate components (low-frequency components) of the L-th level wavelet decomposition. It is an L1 norm. This represents the weighting coefficients of the wavelet approximation component (low-frequency component).

[0071] The optimization solution employs the Alternating Direction Multiplier Method (ADMM) or the Interior Point Method, which decomposes the complex constrained optimization problem into a series of easily solvable subproblems and achieves global optimum through iterative coordination.

[0072] Step 6: Fault Location and Analysis

[0073] Based on the optimized current distribution of each layer, a comprehensive fault analysis and precise fault location are performed, generating an intuitive diagnostic report for engineers' decision-making reference. First, the current density amplitude distribution of each layer is calculated. It then compares the data point by point with a reference pattern that is in normal working condition, and identifies anomalous areas with significant deviations through statistical hypothesis testing.

[0074] Anomaly detection employs multiple criteria: when the local current density exceeds the dynamic threshold of each layer, it is marked as a potential fault point, where the dynamic threshold of the i-th layer... , and Let be the mean and standard deviation of the current in the i-th layer, and let be the time markers as potential fault points;

[0075] Calculate the spatial gradient of current density , Indicates Fourier transform, This represents the inverse Fourier transform. It is a frequency domain high-pass filter core. , These are parameters used to control the high-pass cutoff frequency. Among them... The parameters used to control the Qualcomm cutoff frequency typically take the following values: , The Nyquist limit is the sampling frequency. Abrupt spatial gradient changes often correspond to open-circuit or short-circuit faults; abnormal high-frequency components are identified through spectral analysis, which typically indicate partial discharge or electromigration. Cluster analysis is performed on the identified abnormal regions, merging adjacent fault points to form fault areas, and classifying them according to fault characteristics, such as open-circuit faults exhibiting near-zero current density, short-circuit faults exhibiting abnormally high current density, and poor welding exhibiting uneven current distribution. The final output is a comprehensive diagnostic report including the fault layer number, precise two-dimensional coordinate location, fault type, severity score, and recommended remedial measures. It also generates pseudo-color images of the current distribution in each layer and 3D visualization results for intuitive understanding of the fault location and its impact range.

[0076] Example

[0077] This embodiment verifies the effectiveness of the method of the present invention through numerical simulation, and constructs a test model containing two conductive layers, such as... Figure 2As shown, the first layer is a loop current circuit with a radius of 20 micrometers and a width of 5 micrometers. The second layer is a straight conductor with a width of 4 micrometers and an open circuit defect introduced in the y-coordinate range of 40-60 micrometers. The depths of the two layers are set to 5 micrometers and 10 micrometers, respectively, with a fixed layer spacing Δz = 5 micrometers. The simulation area is 100×100 micrometers, using a 100×100 pixel grid. The magnetic field generated by the current in each layer is calculated using the Biot-Savart law, and the total magnetic field is obtained by superimposing the results. Added standard deviation of 1×10⁻ 9 Tesla's Gaussian noise simulates the actual measurement environment. The geometric series transfer function method of this invention is applied to construct the transfer function matrix in the frequency domain. and Verify its ratio As a constant, the conjugate gradient method is used for joint optimization, and the regularization parameter is set to λ1 = 2 × 10. -29 (Physical constraint weights), λ² = 8 × 10 -27 (Structural constraint weights, where the total variation TV sub-item weight is 1×10) -20 The parameter differences reflect the differences in the numerical ranges of different constraint terms. After 328 iterations, it converged to a residual of 9.87 × 10⁻⁶. -7 Two layers of current distribution were successfully reconstructed. The reconstruction quality indicators for the first layer, a defect-free region, were MSE = 4.3908e-02, PSNR = 23.57dB, SSIM = 0.8654, and correlation coefficient = 0.8382. The reconstruction quality indicators for the second layer, a defective region, were MSE = 9.1314×10⁻⁶. -2 PSNR = 14.10 dB, SSIM = 0.7012, correlation coefficient = 0.6523. The difference in reconstructed current amplitude between normal and defective conditions was calculated. Figure 3 As shown, the defect center was accurately located (x=0.0μm, y=50.0μm, z=10.0μm) with zero positioning error, verifying the feasibility and high accuracy of multilayer current joint reconstruction using a geometric series transfer function constructed with a fixed interlayer spacing.

Claims

1. A method for IC package fault location based on fixed-gap multilayer current distribution reconstruction, characterized in that, By utilizing the geometric feature of the fixed spacing of conductive layers in IC packages, the internal current distribution of multiple layers can be reconstructed from surface magnetic field measurement data to achieve fault location. This includes the following steps: Steps for acquiring magnetic field data: Acquire magnetic field measurement data at multiple probe points on the surface of the IC package under test; Modeling steps for geometric series transfer functions: Based on the number of conductive layers N and the fixed interlayer spacing d obtained from the IC package specifications, determine the depth of each layer. i is the layer number; establish a system based on the fixed layer spacing d and depth. The geometric series transfer function G is a parameter that characterizes the frequency domain mapping relationship between the single-layer current density and the generated magnetic field, and the ratio of the transfer functions of adjacent layers is constant. Current distribution reconstruction steps: The magnetic field measurement data is inverted and calculated using the geometric series transfer function G to reconstruct the three-dimensional current density distribution of each conductive layer inside the IC package. ; Fault location steps: Based on the reconstructed three-dimensional current density distribution of each layer Identify areas of abnormal current distribution and determine their depth within the layer. Three-dimensional fault location is achieved using planar coordinates.

2. The method as described in claim 1, characterized in that, The geometric series transfer function G in the beam domain The expression is: ； in, These are the x and y components of the spatial frequency k, respectively. Where is the free permeability and radial wavenumber. , It is a natural constant.

3. The method as described in claim 2, characterized in that, The current distribution reconstruction steps include: Preliminary solution of frequency domain decomposition: The magnetic field measurement data is converted to the wavenumber domain. Using the analytical properties of the geometric series transfer function model G, the initial separation and reconstruction of the current contribution of each layer is achieved through recursive filtering, and the initial current density of each layer is obtained. Iterative refinement solution: Starting from the initial current density, the current distribution of each layer is iteratively refined by minimizing the objective function L, which includes data fitting and regularization terms, to obtain the refined current density of each layer. Global optimization: Starting with the refined current density of each layer, a multi-objective optimization problem is constructed and solved, which simultaneously considers the current distribution of all layers, data fitting terms, and regularization constraints, to obtain the globally optimal current density distribution.

4. The method as described in claim 3, characterized in that, In the preliminary solution of frequency domain decomposition, in the frequency domain decomposition sub-step, for each spatial frequency component... Based on its radial wavenumber Based on the comparison results with the preset critical wavenumber, an adaptive reconstruction strategy is selected: when Then, the initial current density is reconstructed by joint inversion; when Then, the initial current density of the multilayer can be directly inverted and reconstructed.

5. The method as described in claim 3, characterized in that, In the iterative refinement solution, the objective function L is expressed as: ； gradient of objective function L Represented as: ； in, It is the L2 norm. According to the transfer function Construct the system matrix. For regularization parameters, Represents the regularization term gradient, Indicates transpose. The residual magnetic field; This is the data fitting term.

6. The method as described in claim 5, characterized in that, Regularization constraints in global optimization include physical constraints. and structural constraints ; Physical constraints Including at least one of the following: current continuity constraint, exponential decay constraint, interlayer coupling constraint, and curl constraint; Structural constraints It includes at least one of the following: total variation constraints, sparsity constraints, smoothness constraints, and wavelet domain constraints.

7. The method as described in claim 5, characterized in that, The constraints on current continuity, exponential decay, interlayer coupling, curl, total variation, sparsity, smoothness, and wavelet domain are as follows: Current continuity constraint ; Exponential decay constraint ; Interlayer coupling constraints ; curl constraint ; in, Let be the current density of the i-th layer. It is the attenuation coefficient. It is the current distribution on the lower surface of the i-th layer, and Current distribution on the upper surface of the (i+1)th layer for of Quantity, for of Quantity, express of Component pairs The partial derivative; Total variational constraints ; sparsity constraints ; for The z-component; Smoothness constraint ; Wavelet domain constraints ; in, Let L represent the weighting coefficients of the high-frequency components of the l-th wavelet layer, and L represent the total number of wavelet decomposition layers. This represents the high-pass filter operator in the x-direction when performing wavelet decomposition at the l-th level. This represents the y-direction high-pass filter operator when performing wavelet decomposition at the l-th level. This represents the low-frequency component operator of the L-th level wavelet decomposition. It is an L1 norm. This represents the weighting coefficient of the low-frequency components of the wavelet.

Images