Improved total variation based forward-looking scanning radar super-resolution imaging method

Abstract

The application discloses a forward-looking scanning radar super-resolution imaging method based on an improved total variation, and first establishes an airborne scanning radar echo model; on the basis of ensuring that the Euclidean distance of real echo data and ideal echo data is minimized, a regularization term is added and a regularization model is constructed; then the total variation method is improved in the mode of improving a gradient matrix and increasing a weight coefficient, and an improved total variation regularization imaging model is proposed; finally, an alternating direction multiplier method (ADMM) is adopted to solve a target function in the improved total variation regularization method. The method improves the gradient matrix and increases the weight coefficient, extends adjacent gradient calculation to multiple unit cells, makes the difference between the echo of one pulse and the echoes of multiple adjacent pulses feedback to the position of the original pulse echo in a certain proportion, solves the problems of false boundaries and artifacts existing in the total variation-based super-resolution imaging method, and has better contour preservation performance and imaging quality compared with the total variation method.

Description

Technical Field

[0001] This invention belongs to the field of radar imaging technology, specifically relating to a super-resolution imaging method for forward-looking scanning radar based on improved total variation. Background Technology

[0002] Forward-looking scanning radar (FLSR) uses servo antennas to scan and detect the forward-looking region, but the aperture size of real-aperture radar (PASAR) limits the azimuth resolution of the image. In recent years, an increasing number of algorithms have been proposed to achieve super-resolution imaging with FLSR. Regularization methods, as the most common approach, obtain different imaging characteristics by adding different regularization terms, achieving excellent imaging results in the super-resolution imaging problem of FLSR. The paper "Mihailo Stojnic, 'l2 / l1-optimization in block-sparse compressed sensing and its strongthresholds,' IEEE Journal of Selected Topics in Signal Processing, vol.4, no.2, pp.350-357, 2010" introduces the L1 norm for achieving super-resolution imaging of sparse targets, proposes a block combining L2 and L1 norms to achieve super-resolution imaging of sparse targets, and calculates the lower bound of sparsity for recoverable block sparse signals. The paper “EJ Candes, J. Romberg, and T. Tao, 'Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,' IEEE Transactions on Information Theory, vol. 52, no. 2, pp. 489-509, 2006” introduces a total variational operator that measures the gradient between image pixels. In the image smoothing problem, the adjacent gradient property of this operator can be used to effectively remove noise. The paper “Qiping Zhang, Yin Zhang, Yulin Huang, Yongchao Zhang, Jifang Pei, Qingying Yi, Wenchao Li, and Jianyu Yang, ‘Total variational-sparse super-resolution method for radar forward-looking imaging,’ IEEE Transactions on Geoscience and Remote Sensing, vol.58, no.9, pp.65346549, 2020” introduces the total variational operator as a regularization term in super-resolution image reconstruction, utilizing its smoothness and adjacent gradient properties to preserve image contour information, thus providing a good foundation for radar image target recognition.However, existing super-resolution imaging methods based on total variation regularization often suffer from false boundaries and artifacts, which greatly affect the quality of super-resolution imaging.

[0003] While L1 and L2 norms can avoid the ill-conditioned nature of direct deconvolution target reconstruction and improve the resolution of forward-looking radar super-resolution imaging, they cannot reconstruct target contour information. Although the total variational method has the ability to preserve target contours, it still suffers from false boundaries and artifacts, failing to achieve higher quality imaging. Summary of the Invention

[0004] To address the aforementioned technical problems, this invention provides a forward-looking scanning radar super-resolution imaging method based on improved total variation, which overcomes the problems of false boundaries and artifacts, and has better contour preservation performance and imaging quality.

[0005] The technical solution adopted in this invention is: a super-resolution imaging method for forward-looking scanning radar based on improved total variation, the specific steps of which are as follows:

[0006] Step 1: Airborne scanning radar echo modeling;

[0007] A motion geometry model based on an airborne scanning radar is constructed. The radar platform moves forward at a constant speed v while transmitting a linear frequency modulated signal to detect targets in the area ahead.

[0008] The radar coordinates of the carrier platform are (0,0,H), and the radar beam scans the target at an angular velocity ω.

[0009] For a target P on the ground, its polar coordinate position is represented as (R0,θ,α), where R0 represents the relative distance between the target and the radar, θ represents the azimuth angle of the target relative to the radar, and α represents the elevation angle of the target relative to the radar.

[0010] Based on the movement of the radar platform, the target's range history R(t) is expressed as follows:

[0011]

[0012] Where R0 represents the initial target distance, t represents the distance travel time, and θ0 represents the initial spatial azimuth angle.

[0013] The radar transmits a linear frequency modulated signal and scans it. After pulse compression and range travel correction, the expression for the airborne scanning radar echo s(τ,t) is obtained as follows:

[0014]

[0015] Where t represents the range-fast time, τ represents the azimuth-slow time, σ0 represents the target scattering coefficient, h(t) represents the antenna pattern function, sinc{·} represents the pulse compression response function, B represents the signal bandwidth, c represents the speed of light, λ represents the wavelength, and R0 represents the initial target range.

[0016] The target echo is represented as the convolution of the antenna measurement matrix and the target scattering coefficients, plus noise, as shown in the following expression:

[0017] s=Hx+n (3)

[0018] Where s represents the echo, x represents the target scattering coefficient, and n represents noise. s, x, and n are all N×1 vectors, H is a low-rank convolution matrix of size N×N, and N represents the number of azimuth sampling points.

[0019] During the scanning process, when the antenna enters or leaves the target scene, only half of the antenna beam can illuminate the scene, ignoring the half-beamwidth at the scene edges. Therefore, the measurement matrix contains sample rows with truncated beams, and the expression for H is as follows:

[0020]

[0021] Among them, [h -m … h0 … h m ] indicates the sampling of the antenna pattern, and m indicates the number of sampling points for the half-beamwidth of the antenna pattern.

[0022] Step 2: Construct a regularization model;

[0023] According to equation (3), super-resolution imaging is transformed into a deconvolution problem. The scattering coefficients of the target are reconstructed by minimizing the Euclidean distance between the real echo data and the ideal echo data. The expression is as follows:

[0024]

[0025] Where ||·||2 represents the 2-norm.

[0026] By adding different constraints to equation (5), reconstruction performance with different characteristics can be obtained, as shown in the following expressions:

[0027]

[0028] Where A(x) represents the constraint function acting on x, and λ represents the regularization parameter. p This represents the p-norm.

[0029] Step 3: Construct an improved total variation regularized imaging model;

[0030] The gradient matrix D of the improved total variation method is given by adding weight coefficients. The expression for the gradient matrix D is as follows:

[0031]

[0032] The regularized improved total variation operator expression is as follows:

[0033] ITV(x)=α l D l x1 (8)

[0034] Where α (0 < α ≤ 1) represents the weighting factor, and l represents the echo of l units.

[0035] D l The improved D is expressed as follows:

[0036]

[0037] Each row contains *l* 1s. Different boundary effects are obtained by changing the value of *l*. When *l* = 1 and *α* = 1, it is consistent with the total variation constraint.

[0038] The improved total variation regularized imaging model expression is as follows:

[0039]

[0040] Transform x and s from vector form to matrix form. The optimization problem of equation (10) is transformed as follows:

[0041]

[0042] Where X represents the target scattering coefficient of the entire scene, and S represents the echo of the entire scene.

[0043] Step 4: Solve for the super-resolution imaging results;

[0044] The ADMM method is used to solve equation (11), which decomposes the problem into locally solvable subproblems, solves each subproblem in parallel, and finally coordinates the solutions of these subproblems according to the constraints to obtain the global solution.

[0045] First, we introduce a new variable Z to reformulate the problem, as shown in the following expression:

[0046]

[0047] The constrained optimization problem in equation (12) is transformed into an unconstrained optimization problem by constructing an augmented Lagrangian function, as shown in the following expression:

[0048]

[0049] Where U represents the Lagrange multiplier, p represents the augmented Lagrange parameter, and T represents the matrix transpose.

[0050] Then, the optimization problem of equation (13) is decomposed into three sub-problems, as shown in the following expressions:

[0051]

[0052] Where i represents the iteration number, X i+1 Z i+1 U i+1 This represents the result of the (i+1)th iteration of variables X, Z, U.

[0053] In each iteration, each subproblem is optimized sequentially. At each step, two other variables are fixed while one variable is updated. By alternating these steps, all variables are progressively optimized. After i+1 iterations, the target scattering coefficient X is obtained. i+1 The iterative optimization process is as follows:

[0054] A1. Fix Z and U, then L p Let L be a function that is only related to X, and let L be a function that makes L = X. p Update X using the smallest possible point;

[0055] L p If it is differentiable for X, then directly let The expression for obtaining the updated value of X is as follows:

[0056]

[0057] The expression for obtaining the updated value of X is as follows:

[0058] X=(H T H+pD l T D l ) -1 (H T S+pD l T (ZU)) (16)

[0059] A2, with U fixed and X updated in step A1, L is now... p Let L be a function that is only related to Z, and let L be a function that makes L = L. p Update Z using the smallest point;

[0060]

[0061] in, This represents the threshold shrinkage operator.

[0062] A3. Based on steps A1 and A2, update U with the updated X and Z;

[0063] U = U + D l (XZ) (18)

[0064] Repeat steps A1-A3 until convergence is achieved. The final solution expression is as follows:

[0065]

[0066] The beneficial effects of this invention are as follows: First, the method of this invention establishes an airborne scanning radar echo model. While minimizing the Euclidean distance between the real echo data and the ideal echo data, a regularization term is added and a regularization model is constructed. Then, the total variation method is improved by modifying the gradient matrix and increasing the weighting coefficients, and an improved total variation regularized imaging model is proposed. Finally, the Alternating Direction Multiplier Method (ADMM) is used to solve the objective function in the improved total variation regularized method. This invention improves the gradient matrix and increases the weighting coefficients, extending the calculation of adjacent gradients to multiple cells. This allows the difference between the echo of one pulse and multiple adjacent pulse echoes to be fed back to the position of the original pulse echo in a certain proportion, solving the problems of false boundaries and artifacts in total variation-based super-resolution imaging methods. Compared with the total variation method, it has better contour preservation performance and imaging quality. Attached Figure Description

[0067] Figure 1 This is a flowchart of a forward-looking scanning radar super-resolution imaging method based on improved total variation according to the present invention.

[0068] Figure 2 This is a schematic diagram of the spatial motion geometry model and angular echo data of an airborne scanning radar in an embodiment of the present invention.

[0069] Figure 3 This is a schematic diagram of the total variation and improved total variation principles in the embodiments of the present invention.

[0070] Figure 4 This is a comparison chart of the imaging results of the simulation experiment in the embodiment of the present invention.

[0071] Figure 5 This is a cross-sectional view of the imaging results in an embodiment of the present invention. Detailed Implementation

[0072] This invention employs simulation experiments to demonstrate the effectiveness of the proposed method. All steps and conclusions of this invention have been verified on the Matlab 2019b simulation platform. To enable those skilled in the art to understand the invention, the method is further described below with reference to the accompanying drawings and embodiments.

[0073] like Figure 1The flowchart of a forward-looking scanning radar super-resolution imaging method based on improved total variation is shown below. The specific steps are as follows:

[0074] Step 1: Airborne scanning radar echo modeling;

[0075] Construct a motion geometry model based on airborne scanning radar. The spatial motion geometry model of airborne scanning radar is as follows: Figure 2 As shown, the radar platform moves forward at a constant speed v while emitting a linear frequency modulated signal to detect targets in the area ahead.

[0076] The radar coordinates of the carrier platform are (0,0,H), and the radar beam scans the target at an angular velocity ω.

[0077] For a target P on the ground, its polar coordinate position is represented as (R0,θ,α), where R0 represents the relative distance between the target and the radar, θ represents the azimuth angle of the target relative to the radar, and α represents the elevation angle of the target relative to the radar.

[0078] Based on the movement of the radar platform, the target's range history R(t) is expressed as follows:

[0079]

[0080] Where R0 represents the initial target distance, t represents the distance travel time, and θ0 represents the initial spatial azimuth angle.

[0081] The radar transmits a linear frequency modulated signal and scans it. After pulse compression and range travel correction, the expression for the airborne scanning radar echo s(τ,t) is obtained as follows:

[0082]

[0083] Where t represents the range-fast time, τ represents the azimuth-slow time, σ0 represents the target scattering coefficient, h(t) represents the antenna pattern function, sinc{·} represents the pulse compression response function, B represents the signal bandwidth, c represents the speed of light, λ represents the wavelength, and R0 represents the initial target range.

[0084] Since the received echo is often affected by system noise and environmental noise, the target echo can be represented as the convolution of the antenna measurement matrix and the target scattering coefficient, plus noise, as shown in the following expression:

[0085] s=Hx+n (3)

[0086] Where s represents the echo, x represents the target scattering coefficient, and n represents noise. s, x, and n are all N×1 vectors, H is a low-rank convolution matrix of size N×N, and N represents the number of azimuth sampling points.

[0087] During the scanning process, when the antenna enters or leaves the target scene, only half of the antenna beam can illuminate the scene, ignoring the half-beamwidth at the scene edges. Therefore, the measurement matrix contains sample rows with truncated beams, and the expression for H is as follows:

[0088]

[0089] Among them, [h -m … h0 … h m ] indicates the sampling of the antenna pattern, and m indicates the number of sampling points for the half-beamwidth of the antenna pattern.

[0090] Step 2: Construct a regularization model;

[0091] According to equation (3), super-resolution imaging is transformed into a deconvolution problem. However, due to the ill-conditioned nature of the antenna measurement matrix, the target cannot be directly reconstructed through deconvolution. Therefore, the method of minimizing the Euclidean distance between the real echo data and the ideal echo data is used to reconstruct the target's scattering coefficients, as expressed below:

[0092]

[0093] Where ·2 represents the 2-norm.

[0094] To achieve better image reconstruction results, different constraints are added to equation (5) to obtain reconstruction performance with different characteristics, as shown in the following expression:

[0095]

[0096] Where A(x) represents the constraint function acting on x, and λ represents the regularization parameter. p This represents the p-norm.

[0097] In forward-looking scanning radar super-resolution imaging, constraints are typically applied using L1 regularization terms with sparsity, L2 regularization terms with smoothness, and total variation regularization terms with contour preservation.

[0098] Step 3: Construct an improved total variation regularized imaging model;

[0099] Unlike L1 and L2 regularization operators, which only consider image sparsity and smoothness, the total variation operator helps preserve the contour information of the image, thus reproducing the shape of the target more realistically. However, the total variation method suffers from false boundaries and artifacts, limiting the quality of super-resolution imaging. To overcome these problems, the gradient matrix D of the total variation method is improved and weighting coefficients are increased.

[0100] The gradient matrix D is expressed as follows:

[0101]

[0102] The regularized improved total variation operator expression is as follows:

[0103] ITV(x)=α l D l x1 (8)

[0104] Where α (0 < α ≤ 1) represents the weighting factor, and l represents the echo of l units.

[0105] D l The improved D is expressed as follows:

[0106]

[0107] Each row contains *l* 1s. Different boundary effects are obtained by changing the value of *l*. When *l* = 1 and *α* = 1, it conforms to the total variation constraint. The schematic diagrams of the total variation and improved total variation are shown below. Figure 3 As shown.

[0108] The improved total variation regularized imaging model expression is as follows:

[0109]

[0110] To consider a more efficient batch processing method, the vector form of x and s is transformed into matrix form. The optimization problem of equation (10) is transformed as follows:

[0111]

[0112] Where X represents the target scattering coefficient of the entire scene, and S represents the echo of the entire scene.

[0113] Step 4: Solve for the super-resolution imaging results;

[0114] Equation (11) is solved using the ADMM method. This method decomposes the problem into locally solvable subproblems, solves each subproblem in parallel, and finally coordinates the solutions of these subproblems according to the constraints to obtain the global solution.

[0115] In order to use the ADMM method, a new variable Z must be introduced to reformulate the problem.

[0116]

[0117] The constrained optimization problem in equation (12) is transformed into an unconstrained optimization problem by constructing an augmented Lagrangian function, as shown in the following expression:

[0118]

[0119] Where U represents the Lagrange multiplier, p represents the augmented Lagrange parameter, and T represents the matrix transpose.

[0120] Then, the optimization problem of equation (13) is decomposed into three sub-problems, as shown in the following expressions:

[0121]

[0122] Where i represents the iteration number, X i+1 Z i+1 U i+1 This represents the result of the (i+1)th iteration of variables X, Z, U.

[0123] In each iteration, each subproblem is optimized sequentially. At each step, two other variables are fixed while one variable is updated. By alternating these steps, all variables are progressively optimized. After i+1 iterations, the target scattering coefficient X is obtained. i+1 The iterative optimization process is as follows:

[0124] A1. Fix Z and U, then L p Let L be a function that is only related to X, and let L be a function that makes L = X. p Update X using the smallest possible point;

[0125] L p If it is differentiable for X, then directly let The expression for obtaining the updated value of X is as follows:

[0126]

[0127] The expression for obtaining the updated value of X is as follows:

[0128] X=(H T H+pD l T D l ) -1 (H T S+pD l T (ZU)) (16)

[0129] A2, with U fixed and X updated in step A1, L is now... p Let L be a function that is only related to Z, and let L be a function that makes L = L. p Update Z using the smallest point;

[0130]

[0131] in, This represents the threshold shrinkage operator.

[0132] A3. Based on steps A1 and A2, update U with the updated X and Z;

[0133] U = U + D l (XZ) (18)

[0134] Repeat steps A1-A3 until convergence is achieved. The final solution expression is as follows:

[0135]

[0136] In this embodiment, simulation verification was performed on the Matlab2019b simulation platform. The radar simulation parameters are set as shown in Table 1. In order to simulate a real low signal-to-noise ratio environment, 10dB of noise was added to the simulation.

[0137] Table 1

[0138] Simulation parameters numerical values carrier frequency 10GHz bandwidth 30MHz Scan speed 60° / s Scan range ±5° beamwidth 3° Pulse repetition frequency 4000Hz Signal-to-noise ratio 10dB

[0139] Figure 4 This is a comparison chart of the imaging results from the simulation experiment in this embodiment. Figure 4 (a) is the original target scene, which contains two rectangular targets that are equidistant but have different azimuth angles. Figure 4 (b) is an echo, while Figure 4 (c) shows the imaging results after processing with the total variation method. Due to noise, the image is affected. Figure 4 In (c), many false targets appeared in the non-target area. At the same time, the shape of the processed target was poorly preserved. Figure 4 (d) shows the imaging results of the improved total variation method line-by-line processing. Although there are a few artifacts in the image, no false bright targets appear, and the shape of the target is well preserved. Figure 4 (e) To improve the batch imaging results of the total variation method, and Figure 4 The results of line-by-line processing in (d) are the same.

[0140] To accurately evaluate the imaging performance and contour preservation capability of the proposed method, this embodiment extracts contours from the imaging results, such as... Figure 5 As shown. Figure 5 (a) is the outline of the original scene. Figure 5 (b) is a distance cell profile of the echo. Figure 5 (c) shows the contour after processing with the total variation method, where the shape of the target on the right becomes stepped, and false targets with higher amplitudes appear. Figure 5 (d) To improve the cross-sectional view processed by the total variation method, the shape of the target remains intact. Although some small artifacts appear due to the influence of noise, their magnitudes are all lower than that of the target itself. Figure 5 (e) To improve the profile after batch processing using the total variation method, and Figure 5 The processing results in (d) are the same.

[0141] In this embodiment, batch processing takes only 0.4965 seconds to process an echo matrix of size 148×667, while line-by-line processing takes 15.4166 seconds. Although the imaging results obtained by line-by-line processing and batch processing are similar, using batch processing can shorten the processing time and improve imaging efficiency.

[0142] In summary, simulation experiments demonstrate that the method of this invention outperforms existing total variation regularized imaging methods, effectively eliminating false boundaries and artifacts under low signal-to-noise ratio conditions. The method of this invention improves the gradient matrix and increases the weighting coefficients, extending the calculation of adjacent gradients to multiple cells. This allows the difference between the echo of one pulse and the echoes of multiple adjacent pulses to be fed back to the position of the original pulse echo in a certain proportion, solving the problem of false boundaries and artifacts in total variation-based super-resolution imaging methods. Compared with total variation methods, it exhibits better contour preservation performance and imaging quality.

[0143] Those skilled in the art will recognize that the embodiments described herein are intended to help the reader understand the principles of the invention, and should be understood that the scope of protection of the invention is not limited to such specific statements and embodiments. Various modifications and variations can be made to the invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the invention should be included within the scope of the claims of the invention.

Claims

1. An improved total variation based forward-looking scanning radar super-resolution imaging method, the specific steps are as follows: Step one, modeling of airborne scanning radar echo; A kinematic model is constructed based on an airborne scanning radar, the radar platform moves at a constant speed Proceeding while transmitting a linear frequency modulated signal to detect targets in the area ahead wherein, The radar coordinate position of the carrier platform is , and the radar beam scans the target at an angular velocity . For a target P on the ground, its polar coordinate position is denoted as , denotes the relative distance between the target and the radar, denotes the azimuth angle of the target relative to the radar, denotes the elevation angle of the target relative to the radar; Distance history of a target according to motion of the radar platform The expression is as follows: ； wherein, denotes the target initial distance, denotes the distance-to-go, denotes the initial spatial azimuth angle; The radar transmits a linear frequency modulated signal and scans to obtain an airborne scanning radar echo after pulse compression and range walk correction The expression is as follows: ； wherein, denotes the range fast time, denotes the azimuth slow time, denotes the target scattering coefficient, denotes the antenna pattern function, denotes the pulse compression response function, denotes the signal bandwidth, denotes the speed of light, denotes the wavelength, denotes the target initial range; The target echo is expressed as the convolution superposition of antenna measurement matrix and target scattering coefficient plus noise, and the expression is as follows: ； wherein, denotes an echo, denotes a target scattering coefficient, denotes a noise; , , are N x 1 vectors, is a low-rank convolution matrix of size N x N, N denoting the number of azimuth samples. In the scanning process, when the antenna enters or leaves the target scene, only half of the antenna beams can irradiate the scene, ignoring the half-beam width at the scene edge; then there are beam-truncated sampling rows in the measurement matrix, The expression is as follows: ； wherein, denotes a sampling of an antenna pattern, denotes a number of antenna pattern half-beamwidth sampling points; Step two, construction of regularization model; According to formula (3), the super-resolution imaging is converted into an inverse convolution problem, and the method of minimizing the Euclidean distance between the real echo data and the ideal echo data is used to reconstruct the scattering coefficient of the target, and the expression is as follows: ； wherein denotes the 2-norm; Different constraint conditions are added in formula (5) to obtain different characteristics of the reconstruction performance, and the expression is as follows: ； wherein, represents a constraint function acting on , represents a regularization parameter, represents a p-norm; Step three, construction of improved total variation regularization imaging model; Improved gradient matrix for total variation method and increasing the weight coefficient; gradient matrix The expression is as follows: ； The expression of the improved total variation regularization operator is as follows: ； wherein 0 ≤ 1 represents a weight factor, represents echoes; For the improved , the expression is as follows: ； where each row has ones; different boundary effects are obtained by varying ; when = 1, = 1, it is consistent with the total variation constraint; The expression of the improved total variation regularization imaging model is as follows: ； Transforming the vector form of and into matrix form; transforming the optimization problem of equation (10) as follows: ； wherein, a target albedo representing the entire scene, a return representing the entire scene; Step four, solving the super-resolution imaging result; The ADMM method is used to solve formula (11), the problem is decomposed into subproblems that can be solved locally, then each subproblem is solved in parallel, and finally the solutions of these subproblems are coordinated according to the constraint conditions to obtain the global solution; First introduce a new variable The problem is restated as follows: ； The augmented Lagrangian function is constructed to convert the constraint optimization problem of formula (12) into an unconstrained optimization problem, and the expression is as follows: ； wherein, denotes a Lagrange multiplier, denotes an augmented Lagrange parameter, denotes matrix transposition; Then the optimization problem of formula (13) is decomposed into three subproblems, and the expression is as follows: ； wherein, denotes the number of iterations, denotes the variable the result of the first iteration; In each iteration, the optimization is performed sequentially for each subproblem; each step fixes the other two variables to update one variable; by repeating these steps alternately, all variables are gradually optimized; after the iterations, the target scattering coefficient is obtained; the iterative optimization process is as follows: A1, fix and at this time denotes only with respect to the function, update with the point that minimizes ; To derivable, directly let , get The update value expression of is as follows: ； The updated value expression of the above is as follows: The updated value expression of the above is as follows: ； A2, fix and step Al updated at this time denotes only with respect to the function, update with the point that ; ； wherein denotes a threshold shrink operator; A3. Based on steps Al, A2, update the and update ; ； Steps A1-A3 are repeatedly performed until convergence, and the final solution expression is as follows: 。

Images