A SAR control point extraction method and device, electronic equipment and storage medium

Abstract

The present application relates to the technical field of remote sensing image processing, and particularly relates to a SAR control point extraction method and device, electronic equipment and storage medium, wherein the SAR control point extraction method comprises: acquiring multi-scene SAR images and matched candidate point observation data and imaging geometric parameters; establishing a multi-view positioning observation model; performing initial three-dimensional coordinate calculation to obtain a primary calculation result; combining prior constraints of imaging geometric configuration to perform variance component estimation and dynamically optimize observation weights of each scene image; recalculating three-dimensional coordinates; iteratively performing weight optimization and coordinate calculation until a convergence condition is met; and screening and determining SAR control points according to an accuracy index of a final calculation result. The above method constructs a complete framework of robust calculation of multi-view stereo SAR control points, embeds geometric prior constraints such as intersection angles and orbit combinations into the variance component estimation process, and realizes end-to-end optimization processing from candidate homonymic points to high-reliability control points.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image processing technology, and more specifically, to a SAR control point extraction method, apparatus, electronic device, and storage medium. Background Technology

[0002] Synthetic Aperture Radar (SAR) satellite remote sensing technology can acquire high-resolution microwave imagery of the Earth's surface. In urban and suburban environments, man-made features such as metal lampposts, traffic signs, and metal water tanks appear as high-amplitude, high-signal-to-noise-ratio pixel-level bright spots in SAR imagery due to their strong backscattering characteristics, good temporal stability, and near-pointy geometry. These features are easily detected automatically and matched across scenes. These strong scattering points can be extracted as SAR control points using spaceborne SAR stereo measurement technology. In the absence of ground control points or under weak control conditions, they serve as an effective supplement to traditional field measurement control points, supporting high-precision terrain reconstruction and the generation of geographic information products.

[0003] For matched SAR candidate points with the same name from multiple scenes, existing 3D positioning methods mainly construct and solve range-Doppler observation equations. However, existing technologies suffer from a lack of physical priors in their weighting strategies: they typically employ static equal-weighting or purely data-driven variance component estimation weighting mechanisms, failing to incorporate the geometric quality and orbital combination characteristics of the observations into the weighting model. This results in the solution process failing to adaptively reflect the reliability of the observations. Although variance component estimation can dynamically adjust the weights, it relies entirely on the observation residuals. In cases where strong scattering points have small residuals but weak geometric conditions, they are easily misjudged as high-precision observations, leading to a significant decrease in solution accuracy under weak geometric conditions. Furthermore, in multi-scene configurations, the imbalance in weight allocation causes degradation of internal consistency performance, contradicting the expectation of improving accuracy through multi-scene setups. Summary of the Invention

[0004] The purpose of this application is to provide a SAR control point extraction method, apparatus, electronic device, and storage medium to solve the above-mentioned technical problems.

[0005] In a first aspect, embodiments of this application provide a SAR control point extraction method, the method comprising: acquiring multiple SAR images and observation data and imaging geometric parameters of matched candidate points in the multiple SAR images; establishing a multi-view positioning observation model based on the data and parameters; performing initial three-dimensional coordinate calculation on the observation model to obtain initial calculation results; performing variance component estimation based on the initial calculation results and prior constraints of imaging geometric configuration to dynamically optimize the observation weights of each image; recalculating the three-dimensional coordinates based on the optimized observation weights; iteratively executing the weight optimization step and coordinate calculation step until the convergence condition is met; and selecting and determining SAR control points based on the accuracy index of the final calculation results.

[0006] Furthermore, the prior constraints of the imaging geometry configuration are used to perform variance component estimation, which includes: for each image, calculating the intersection angle based on the spatial relationship between the satellite position vector of the candidate point at the imaging time of the image and the satellite position vectors at the corresponding times of the other images, and determining the geometric robustness level of the stereo image pair formed by the image and the other images based on the intersection angle, and selecting the highest level as the comprehensive geometric robustness level of the image; determining the orbit combination type between images based on the orbital affiliation relationship of each image, introducing a co-orbital system error penalty factor for co-orbital image combinations to suppress the influence of co-orbital system error correlation on the solution results; constructing a prior variance model based on the comprehensive geometric robustness level and the co-orbital system error penalty factor, and introducing the prior variance model as a regularization term into the objective function of variance component estimation.

[0007] Furthermore, the objective function for the variance component estimation includes an observation residual term and a regularization term based on the prior variance model; the dynamic optimization of the observation weights of each image scene includes: obtaining the optimized variance components of each image scene by minimizing the objective function, and generating an adaptive weight matrix based on the optimized variance components to dynamically optimize the observation weights of each image scene.

[0008] Further, the step of selecting and determining SAR control points based on the accuracy index of the final solution result includes: calculating the variance-covariance matrix at the final solution three-dimensional coordinates according to the type of three-dimensional coordinate solution model used, and extracting the mean square error of the candidate point in the plane and elevation directions from the diagonal elements of the variance-covariance matrix; defining the intra-point conformance accuracy based on the mean square error in the plane and the mean square error in the elevation direction, wherein the intra-point conformance accuracy includes the intra-plane conformance accuracy and the intra-elevation conformance accuracy; determining whether the intra-plane conformance accuracy and the intra-elevation conformance accuracy both meet the preset threshold conditions; if they meet the conditions, the candidate point is determined as a SAR control point and output; if they do not meet the conditions, the candidate point is discarded as an invalid point.

[0009] Further, the initial three-dimensional coordinate calculation of the observation model includes: based on the multi-view positioning observation model, the initial equal-weight matrix is ​​used in combination with the initial three-dimensional coordinates of the candidate points for the first calculation; the first calculation is performed using either a linearized model or a nonlinear model: if a linearized model is used, the observation model is expanded in first order Taylor at the initial three-dimensional coordinates to construct a design matrix and constant terms, and a weighted least squares problem is solved to obtain coordinate corrections to update the initial calculation results; if a nonlinear model is used, the initial three-dimensional coordinates are used as initial values, and the weighted residual objective function is minimized using an iterative optimization algorithm with damping factors. The iterative optimization process includes calculating the residual vector and partial derivative matrix, constructing the damping normal equation and solving for the coordinate corrections, and outputting the initial calculation results after multiple rounds of iterative optimization.

[0010] Furthermore, the convergence condition includes at least one of the following: the coordinate update amount between the coordinate solution result of the current iteration step and the coordinate solution result of the previous iteration step is less than a preset convergence threshold; the current iteration step number reaches a preset maximum number of iterations.

[0011] Furthermore, the establishment of the multi-view positioning observation model based on the data and parameters includes: for each image scene, constructing a distance observation equation and a Doppler observation equation based on the image coordinates of the candidate point in the image, the satellite position vector and velocity vector corresponding to the image imaging time; combining the distance observation equation and the Doppler observation equation of each image scene to form a set of nonlinear equations with the three-dimensional geographic coordinates of the candidate point as unknowns, which serves as the multi-view positioning observation model.

[0012] Secondly, embodiments of this application provide a SAR control point extraction device, comprising: The image acquisition module is used to acquire multiple SAR images and the observation data and imaging geometric parameters of the matched candidate points in the multiple SAR images; The observation model building module is used to build a multi-view positioning observation model based on the data and parameters. The initial solution module is used to perform initial three-dimensional coordinate solution on the observation model and obtain the initial solution result; The adaptive optimization module is used to perform variance component estimation based on the initial solution results and in combination with the prior constraints of the imaging geometry configuration, and dynamically optimize the observation weights of each image scene. The secondary solution module is used to recalculate the three-dimensional coordinates based on the optimized observation weights; The control point selection module is used to iteratively execute the weight optimization step of the adaptive optimization module and the coordinate calculation step of the secondary calculation module until the convergence condition is met; and to select and determine SAR control points based on the accuracy index of the final calculation result.

[0013] Thirdly, embodiments of this application provide an electronic device, the device including: a memory and a processor, wherein the memory and the processor communicate with each other through an internal connection path, the memory is used to store instructions, the processor is used to execute the instructions stored in the memory, and when the processor executes the instructions stored in the memory, the processor causes the processor to perform the method in any of the above-described embodiments.

[0014] Fourthly, embodiments of this application provide a computer-readable storage medium that stores a computer program, wherein when the computer program is run on a computer, the methods in any of the above-described embodiments are executed.

[0015] The advantages or beneficial effects of the above technical solutions include at least the following: First, by constructing a robust solution framework for multi-view stereo SAR control points, and by embedding geometric prior constraints such as intersection angle and orbit combination into the variance component estimation process, end-to-end optimization processing from candidate corresponding points to high-reliability control points is achieved. Second, by constructing a prior variance model based on geometric robustness level and error penalty factor of the same track system, and introducing it as a regularization term into the objective function of variance component estimation, the problem of weak geometric observation misjudgment caused by lack of physical prior in traditional weight setting is solved, and negative variance estimation and iterative divergence are avoided. Third, the use of an adaptive weight matrix for iterative coordinate calculation effectively suppresses the weight contribution of small intersection angles and combined observations on the same orbit in the calculation process, overcomes the problem of accuracy degradation caused by blindly introducing low-quality observations in multi-scene SAR configuration, and improves the robustness of three-dimensional positioning under weak geometric conditions. Fourth, the generated adaptive weight matrix can be seamlessly adapted to the dual-path solution framework of linearized weighted least squares and nonlinear iterative optimization, taking into account both computational efficiency and accuracy requirements. It is suitable for different hardware platforms and processing flows and has strong compatibility. Fifth, based on the variance-covariance matrix of the final solution, the mean square error in the plane and elevation directions is extracted. Through the point-to-point accuracy threshold screening mechanism, the automated quality assessment and invalid point removal of SAR control points are realized, ensuring the accuracy and reliability of the output control points.

[0016] Other features and advantages of the invention will be set forth in the following description, and will be apparent in part from the description, or may be learned by practicing the embodiments of this application. The objects and other advantages of the invention may be realized and obtained by means of the structures particularly pointed out in the written description, claims, and drawings. Attached Figure Description

[0017] To more clearly illustrate the technical solutions of the embodiments of this application, the accompanying drawings used in the embodiments of this application will be briefly introduced below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as a limitation on the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0018] Figure 1 A flowchart illustrating the SAR control point extraction method provided in this application embodiment; Figure 2 The experimental results diagram provided in this application embodiment reflects the systematic error of the geometric measurement model for different incident angle differences under different intersection angles; Figure 3 This is a schematic diagram of the SAR control point extraction device provided in the embodiments of this application; Figure 4 This is a schematic diagram of the structure of an electronic device provided in an embodiment of this application. Detailed Implementation

[0019] The embodiments of the technical solution of the present invention will be described in detail below with reference to the accompanying drawings. The following embodiments are only used to more clearly illustrate the technical solution of the present invention, and are therefore only examples, not intended to limit the scope of protection of the present invention. Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains; the terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit the invention; the terms "comprising" and "having," and any variations thereof, in the specification, claims, and the foregoing description of the accompanying drawings, are intended to cover non-exclusive inclusion. In the description of the embodiments of this application, technical terms such as "first," "second," etc., are only used to distinguish different objects and should not be construed as indicating or implying relative importance or implicitly indicating the number, specific order, or primary or secondary relationship of the indicated technical features. In the description of the embodiments of this application, "a plurality of" means two or more, unless otherwise explicitly specified. The reference to "embodiment" herein means that a particular feature, structure, or characteristic described in connection with an embodiment may be included in at least one embodiment of the present invention. The appearance of this phrase in various places in the specification does not necessarily refer to the same embodiment, nor is it an independent or alternative embodiment mutually exclusive with other embodiments. It will be explicitly and implicitly understood by those skilled in the art that the embodiments described herein can be combined with other embodiments.

[0020] Long-term stable man-made features in urban and suburban environments (such as metal lampposts, traffic signs, and metal water tanks) naturally appear as high-amplitude, high-signal-to-noise-ratio pixel-level bright spots in SAR imagery due to their strong backscattering characteristics, good temporal stability, and near-point geometric shape. These features are easily detected automatically and matched across scenes. These strong scattering points can be extracted as SAR ground control points (GCPs) using spaceborne SAR stereo measurement technology. In the absence of ground control points or under weak control conditions, these GCPs serve as an effective supplement to traditional field measurement control points, supporting high-precision terrain reconstruction and the generation of geographic information products.

[0021] For matched SAR points of the same name (i.e. candidate control points), conventional three-dimensional positioning methods construct no less than 4 sets of range-Doppler observation equations and solve for their three-dimensional coordinates. The current mainstream solution methods mainly include two categories: (1) Weighted least squares adjustment based on linearized range-Doppler model: usually adopts equal weight processing, or combines variance component estimation (VCE) for dynamic weighting; (2) Nonlinear optimization method based on spatial triangulation (such as Levenberg–Marquardt algorithm, LMA): achieves three-dimensional positioning by minimizing the distance from the ground point to the line-of-sight vector of each SAR scene, but usually assumes that all observations are equally reliable, i.e., implicitly assumes equal weight.

[0022] However, the above methods share the following common drawbacks in their solution strategies:

[0023] (1) Lack of physical priors in weight setting: Existing methods do not incorporate the geometric quality of observations (such as intersection angles) and orbital combination characteristics into weight modeling. Static equal weights or purely data-driven VCE weighting mechanisms are commonly used, which cannot adaptively reflect the reliability of observations during the solution process. Although VCE can dynamically adjust weights, it relies entirely on observation residuals. In cases where strong scattering points have small residuals but weak geometric conditions (such as small intersection angles), they are easily misjudged as high-precision observations, leading to negative variance estimation or iteration non-convergence.

[0024] (2) The accuracy of the solution is significantly reduced under weak geometry conditions: When multiple images contain small intersection angles (e.g., <20°) or combinations of the same orbit, the sensitivity of the imaging geometry to elevation is reduced, and the systematic errors (e.g., orbital deviations) are highly correlated, making the elevation solution fragile. However, traditional methods do not apply a priori suppression to such observations and still assign them high weights, thus contaminating the overall solution results.

[0025] (3) The phenomenon of “pseudo-gain” in multiple scenes is prominent: In typical multiple scene configurations (such as 2 ascending orbits + 1 descending orbit), due to the failure to effectively distinguish the observation quality, the introduction of additional low-precision observations not only fails to improve accuracy, but also leads to a decline in internal consistency performance due to the imbalance of weight distribution, which violates the expectation of “multiple scenes improve accuracy”.

[0026] In summary, currently no technology integrates the physical priors of intersection angle and orbit combination into weighted modeling and combines it with an internal coincidence accuracy screening mechanism to achieve robust solution throughout the entire process from matched SAR candidate corresponding points to high-reliability control points. This invention is a technical solution proposed to address the above-mentioned problems.

[0027] To address the issues of poor robustness and unreasonable weight allocation leading to pseudo-gain and accuracy degradation in existing spaceborne SAR stereo measurement control point extraction methods under weak geometric conditions, this application provides a SAR control point extraction method. By constructing a complete framework for robust multi-view stereo SAR control point extraction, and embedding geometric prior constraints such as intersection angle and orbit combination into the variance component estimation process, end-to-end optimization processing from candidate corresponding points to high-reliability control points is achieved.

[0028] like Figure 1 As shown, the SAR control point extraction method provided in this application embodiment may include: Step S110: Acquire multiple SAR images and the observation data and imaging geometric parameters of the matched candidate points in the multiple SAR images.

[0029] The aforementioned multi-scene SAR imagery refers to a collection of radar remote sensing images with overlapping coverage acquired by a spaceborne synthetic aperture radar satellite at different time phases, orbital directions, or incident angles for the same target area. This collection must contain at least two images. The acquisition of such images relies on the raw echo data acquisition and imaging processing workflow of the spaceborne SAR system. Specifically, range pulse compression, azimuth synthetic aperture focusing, and geometric correction are performed on the raw radar echo signals transmitted from the satellite to generate slant range or ground distance image products containing precise imaging geometric parameters. During data acquisition, the precise orbital state vector, radar wavelength, pulse repetition frequency, Doppler center frequency, and image metadata information corresponding to each image must be extracted simultaneously. The image coordinates of candidate control points in each image are also recorded. These image coordinates are composed of the range sampling time and the azimuth imaging time, forming the basic observation input for subsequent multi-view stereo positioning.

[0030] The candidate point observation data and imaging geometric parameters matched in the aforementioned multi-scene SAR images refer to the measurement information and image auxiliary parameter set derived from the corresponding strong scattering points automatically detected and established from two or more SAR images. The candidate point observation data includes the image coordinates of the candidate points in each image. These image coordinates are composed of range-time and azimuth-time components. The range-time component is obtained by converting the candidate point's image column coordinates with the near-range imaging time and range sampling rate in the image metadata. The azimuth-time component is obtained by converting the candidate point's image row coordinates with the zero-Doppler time and pulse repetition frequency in the image metadata. The imaging geometric parameters include the precise orbit parameters, radar system parameters, and image metadata corresponding to the candidate point imaging time in each image. The precise orbit parameters cover the satellite position vector and velocity vector, representing the sensor's three-dimensional spatial coordinates and operating status at the imaging time, respectively. The radar system parameters include the center wavelength and propagation speed of the radar transmitted signal. The image metadata records the range sampling frequency, pulse repetition frequency, and Doppler center frequency information. The aforementioned observation data and imaging geometric parameters together constitute the input basis of the range-Doppler positioning model.

[0031] Step S120: Establish a multi-view positioning and observation model based on data and parameters.

[0032] The aforementioned multi-view positioning and observation model is essentially a nonlinear equation system composed of the range-Doppler equations from multiple SAR images. This model uses the 3D geographic coordinates of candidate points as unknown parameters and the image coordinates of corresponding points in multiple images as observation inputs. It achieves multi-view geometric convergence positioning by constructing at least four sets of range-Doppler constraints. Specifically, for each SAR image involved in the solution, the range sampling time and azimuth imaging time are calculated based on the image row and column coordinates of the candidate points in that image. Combined with the satellite position and velocity vectors, radar wavelength, and light speed parameters at the corresponding imaging time, range observation equations characterizing the slant range from the candidate point to the sensor and Doppler observation equations characterizing the relative motion relationship are established. After simultaneously solving the 2N range-Doppler equations corresponding to N images, a nonlinear equation system is formed with the 3D geographic coordinates of the candidate points as the unknowns to be solved. This equation system establishes a mathematical relationship between the image spatial observation data and the object space 3D coordinates through multi-view geometric constraints, forming the observation basis for subsequent 3D positioning solutions.

[0033] Optionally, step S120 may include: for each image, constructing a distance observation equation and a Doppler observation equation based on the image coordinates of the candidate point in the image, the satellite position vector and velocity vector corresponding to the image imaging time; combining the distance observation equation and the Doppler observation equation of each image to form a set of nonlinear equations with the three-dimensional geographic coordinates of the candidate point as unknowns, as a multi-view positioning observation model.

[0034] The above step S120, based on the matched SAR candidate point image coordinates, precise orbit parameters for each scene, and image metadata, establishes 2N range-Doppler observation equations for N (N≥2) SAR images, forming a nonlinear observation model for multi-view three-dimensional positioning, and creating a multi-view geometric constraint system, mainly including: For each SAR image i (i=1,2,…,N), the row and column coordinates (u) of the candidate points in that image are combined. i ,v i According to its imaging time t az,i Corresponding satellite position r s,i (t) az,i ) = [X s Y s Z s ] T With velocity v s,i (t az,i )=[V x,s V y,s V z,s ] T The following two nonlinear observation equations are established: ; Where, r t =[X, Y, Z] T λ represents the three-dimensional geographic coordinates of the candidate point to be solved; c is the speed of light, and λ is the radar wavelength. and For distance and Doppler observation noise terms; t rg,i From distance to time, from image column coordinates v i Combine image metadata (such as proximity time, distance sampling rate f) s ) is obtained by conversion; t az,i For azimuth time, derived from the image row coordinate u i It is calculated by combining image metadata (such as zero Doppler time and pulse repetition frequency, PRF); The center frequency of the Doppler vector is usually set to 0 in SAR image processing (i.e., satisfying the zero Doppler condition). In this case, the second equation simplifies to the velocity vector being orthogonal to the line-of-sight vector.

[0035] By combining the 2N distance-Doppler equations of the above N images, a system is formed with r t A system of nonlinear equations with unknowns is used as the observation model for subsequent iterative optimization, enabling high-precision three-dimensional positioning of SAR control points.

[0036] Step S130: Perform initial three-dimensional coordinate calculation on the observation model to obtain the initial calculation results.

[0037] The aforementioned initial 3D coordinate solution is the first parameter estimation process for the multi-view positioning observation model using an initial equal-weight matrix, aiming to provide initial values ​​for subsequent adaptive weight optimization. This solution uses the initial 3D geographic coordinates of candidate points as the starting point and is completed through two optional paths: The first is the linearization model path, which involves performing a first-order Taylor expansion of the distance-Doppler observation equation at the initial coordinates, constructing a design matrix composed of the partial derivatives of the observation equations to be solved for the coordinates, and a constant term composed of the difference between the observed values ​​and the model-calculated values. Then, a weighted least squares problem is solved to obtain the coordinate correction vector, which is superimposed onto the initial coordinates to obtain the initial solution result. The second is the nonlinear model path, which uses the initial coordinates as the initial iteration value and employs an iterative optimization algorithm with a damping factor to minimize the weighted residual objective function. In each iteration, the residual vector between the observed values ​​and the current coordinate estimates, as well as the partial derivative matrix of the observation equations with respect to the coordinates, are calculated. Based on this, a damped normal equation is constructed and solved to obtain the coordinate correction amount. After multiple iterations until the convergence condition is met, the initial solution result is output. Neither of the two paths introduced geometric prior constraints, and both adopted an equal-weighting strategy. Their outputs served as the basic inputs for subsequent geometrically constrained variance component estimation.

[0038] Optionally, step S130 above may include: based on the multi-view positioning observation model, performing the first solution using an initial equal-weight matrix combined with the initial three-dimensional coordinates of candidate points; the first solution is performed using either a linearized model or a nonlinear model: if a linearized model is used, a first-order Taylor expansion is performed on the observation model at the initial three-dimensional coordinates to construct a design matrix and constant terms, and the weighted least squares problem is solved to obtain coordinate corrections to update the initial solution results; if a nonlinear model is used, the initial three-dimensional coordinates are used as initial values, and the weighted residual objective function is minimized using an iterative optimization algorithm with damping factors. The iterative optimization process includes calculating the residual vector and partial derivative matrix, constructing the damping normal equation and solving for the coordinate corrections, and outputting the initial solution results after multiple rounds of iterative optimization.

[0039] In the above scheme, the initial weight matrix can be set as an identity matrix (i.e., equally weighted) P. (0) =I 2N Combined with the initial three-dimensional coordinates r of the candidate points t (0) The first three-dimensional coordinate solution of the candidate control points was completed, and the initial solution result r was obtained. t (1) Specifically: If a linearized model is used, a weighted least squares solution is performed based on the linearized range-Doppler equations. The nonlinear range-Doppler observation equations established in step 1 are then applied to the initial coordinates r. t (0)A first-order Taylor expansion is performed to obtain the linearized observation equation, and then weighted least squares are used to solve for the coordinate correction vector. : ; Among them, B (0) ∈R 2N×3 For in r t (0) The design matrix calculated at the location is the observation equation with respect to r. t Partial derivatives of L; (0) =yf(r t (0) ) represents a constant term (y is the observation vector converted from image coordinates); v lin This is the linearized residual vector, representing the deviation between the observed values ​​and the predicted values ​​of the linear approximation model; This is the initial weight matrix.

[0040] Update 3D coordinates: This process is a single matrix solution and does not include internal iterations.

[0041] If a nonlinear model is used, the objective function is solved by minimizing the weighted residuals using the LMA algorithm. (Using coordinate r...) t (0) Using the initial values, perform LMA inner iteration optimization: (1) Initialization: Let m=0, r t (0,0) =r t (0) Damping factor μ (0) =0.01; (2) Calculate the current residual vector and Jacobian matrix J (m) Construct the damping normal equation and solve for the coordinate correction. : ; Where I3 is a 3×3 identity matrix; This is the Jacobian matrix at the m-th iteration; The damping factor at the m-th iteration; (3) Update coordinates: μ is adaptively adjusted as the objective function decreases. (m) After convergence, the three-dimensional coordinates are output. , This represents the final number of inner iterations. The process involves multiple iterations, i.e., the inner loop.

[0042] In summary, the initial three-dimensional coordinates r t (1)Different solution mechanisms are employed depending on the selected model: in linearized models, the weighted least squares problem is solved directly based on the design matrix and constant terms, without internal iteration; in nonlinear models, the LMA algorithm is used for internal iterative optimization until the convergence condition is met. Regardless of the path used, the output is used as input to the geometrically constrained variance component estimation (GC-VCE) and participates in subsequent outer iterative optimization loops, namely: updating the weight matrix, resolving, and determining overall convergence.

[0043] Iteratively execute steps S140 and S150 above until the convergence condition is met: Step S140: Based on the initial solution results, perform variance component estimation in conjunction with the prior constraints of the imaging geometry configuration, and dynamically optimize the observation weights of each image scene.

[0044] Step S150: Recalculate the three-dimensional coordinates based on the optimized observation weights.

[0045] The prior constraints of the above imaging geometry configuration are qualitative or quantitative constraints constructed based on the spatial geometric relationships during SAR image acquisition. These constraints quantify the influence of the intersection angle and orbit combination type on positioning accuracy, transforming physical geometric characteristics into prior information for weight adjustment.

[0046] The aforementioned variance component estimation refers to the adjustment method in a mixed observation model that determines the variance components based on the residual statistical characteristics and prior constraints of various observations. This method obtains the optimized variance components of each image observation by constructing an objective function containing residual and regularization terms and performing a minimization solution, thereby generating an adaptive weight matrix to achieve differentiated weighting of observations with different geometric qualities.

[0047] Optionally, step S140 above combines prior constraints of imaging geometry configuration to perform variance component estimation, including: for each image, calculating the intersection angle based on the spatial relationship between the satellite position vector of the candidate point at the image imaging time and the satellite position vectors of the other images at the corresponding times, and determining the geometric robustness level of the stereo image pair formed by the image and the other images based on the intersection angle, and selecting the highest level as the comprehensive geometric robustness level of the image; determining the orbit combination type between images based on the orbital affiliation relationship of each image, introducing a co-orbital systematic error penalty factor for co-orbital image combinations to suppress the influence of co-orbital systematic error correlation on the solution results; constructing a prior variance model based on the comprehensive geometric robustness level and the co-orbital systematic error penalty factor, and introducing the prior variance model as a regularization term into the objective function of variance component estimation.

[0048] Optionally, the objective function for the above variance component estimation includes an observation residual term and a regularization term based on the prior variance model; The above step S140 dynamically optimizes the observation weights of each image, including: minimizing the objective function to obtain the optimized variance components of each image, and generating an adaptive weight matrix based on the optimized variance components to dynamically optimize the observation weights of each image.

[0049] The above step S140 uses the initial solution result r t (1) Calculate the current residual vector v lin (1) (Linear case) or v nl (1) (In the nonlinear case), perform geometrically constrained variance component estimation (GC-VCE) to update the variance components of each image observation. Generate an adaptive weight matrix P (1) Specifically: like Figure 2 As shown in the figure, the experimental results reflect the systematic error of the geometric measurement model for different incident angle differences under different intersection angles. Figure 2 Based on real multi-scene SAR image data (SV2-03 / 04), statistical analysis was performed on the residuals of 3D coordinate calculation for a large number of corresponding points under different intersection angles. The horizontal axis represents the intersection angle (unit: degrees), and the vertical axis represents the 3D orientation positioning error (unit: meters). To scientifically reflect the geometric reliability of different stereo image pairs, the embodiments of this application are based on the intersection angle θ. ij The image pairs are divided into three robustness levels: Type I (Strong): θ ij >40°, high geometric strength, and elevation calculation is not sensitive to errors; Class II (Medium): 20°≤θ ij <40°, with moderate positioning capability; Class III (weak): 0° < θ ij <20°, geometrically fragile, positioning errors are easily amplified.

[0050] The experimental results verify the rationality of the robustness classification based on the intersection angle in this invention, and provide an empirical basis for prior modeling in subsequent geometrically constrained variance component estimation (GC-VCE).

[0051] (1) Determine the comprehensive geometric grade of the i-th image: For the i-th image, enumerate all stereo image pairs (including those with different tracks and those with the same track) formed by it and other images, and select the one with the highest intersection angle grade as its comprehensive geometric grade. If there are multiple Class I image pairs, the intersection angle θ closest to 90° can be further selected. i best Used for accurate modeling.

[0052] (2) Calculate the overall geometric quality score of the i-th scene, and define the geometric score as: ; Where k∈[0.03,0.07] is the geometric sensitivity coefficient. This function naturally imposes an exponential penalty on weak intersection angles. The score reflects the optimal solid geometry quality that the i-th scene can participate in. The smaller the value, the stronger the geometry. It is the smallest (geometrically strong) in class I image pairs and the largest (geometrically weak) in class III image pairs.

[0053] (3) Introducing a penalty for co-orbital systematic errors: Although the geometric score already reflects the intersection angle intensity, co-orbital images are geometrically highly similar, resulting in a high correlation of their systematic errors, making them difficult to separate through adjustment. Therefore, to suppress the high correlation of systematic errors between co-orbital images, an additional penalty factor p is introduced. i In the Jing did not participate in any parallel image pairing. ; in the When a scene participates in at least one co-track image pair, Where β∈[0.6,0.8]. This mechanism ensures that even if a scene can participate in high-quality inter-track rendezvous, it will be appropriately downweighted as long as it also participates in same-track pairing, in order to reduce the risk of systematic error propagation.

[0054] (4) Construct the prior variance and perform GC-VCE: Suppose there are N SAR images involved in the calculation. The observations are grouped by image (N categories in total), and the prior variance of each image is calculated. And construct the GC-VCE objective function with prior regularization terms. : ; in, The reference variance (can be set to 1); Q i Let d be the observation structure matrix of the i-th scene, defined as a 2N×2N diagonal matrix. This matrix depends only on the observation group, with its diagonal elements being 1 at the corresponding distances to the Doppler observation positions and 0 elsewhere; i For the effective degrees of freedom, based on the geometric score g i Correcting the effective degrees of freedom: If g i >30 (intersection angle is approximately less than 20°), then let d i Multiply by 0.5, the correction factor is 0.5, preferably in the range of 0.4 to 0.6, and can be adaptively adjusted according to the geometric characteristics of the region; λ∈[0.5,1.0] is the regularization weight.

[0055] minimize To optimize the variance components And generate a 2N×2N diagonal weight matrix: ; The updated weight matrix co-integrates geometric strengths (via a continuous exponential function g). i ) and systematic error risk (through same-track penalty p) i It can achieve physically reasonable and precision-oriented adaptive weighting in various scenarios such as ascent + descent, full ascent orbit, and full descent orbit.

[0056] The above step S150 is based on the updated adaptive weight matrix P (1) The distance-Doppler model equations are solved again in three dimensions to obtain the second three-dimensional coordinate solution for the candidate control points, thus yielding the optimization result r. t (2) Specifically: Based on the solution path selected in step S130 above, perform the following operations respectively: (1) Linear model path: based on the initial solution result r t (1) As a new initial value, the nonlinear distance-Doppler model is relinearized at this point to obtain the updated design matrix B. (1) With constant term L (1) Perform weighted least squares calculation to determine the coordinate corrections: And update the three-dimensional coordinates: .

[0057] (2) Nonlinear Model Path (LMA): with r t (1) =r t (1,0) As the new initial values, the LMA algorithm is used to minimize the weighted residual objective function and solve the damping normal equation: ; The solution is optimized through m inner iterations until the convergence condition is met, and then the convergent solution is output. (The inner layer eventually converges to the first step) step).

[0058] Optionally, the above convergence condition includes at least one of the following: the coordinate update amount between the coordinate solution result of the current iteration step and the coordinate solution result of the previous iteration step is less than a preset convergence threshold; the current iteration step number reaches a preset maximum number of iterations. An example of this implementation is:

[0059] Determine if the convergence condition is met: coordinate update amount ||r t (k+1) -r t (k) ||<ε or the number of iterations reaches the preset maximum value K max If convergence fails, then the current solution result r is used. t (k+1)Based on this, repeat steps 3 and 4, recalculate the residuals and update the weights, and continue iterative optimization. Otherwise, output the final 3D coordinates. (in , (This refers to the final number of iterations). Specifically: ε>0 is the coordinate convergence threshold, typically 0.05; K max This represents the maximum number of iterations, typically 5. If the convergence condition is not met, the current solution result r is used. t (k+1) With the new initial values, return to step 3, recalculate the residual vector, perform geometrically constrained variance component estimation (GC-VCE), update the adaptive weight matrix, and complete the next round of weighted 3D coordinate calculation in step 4, continuing the outer iterative optimization. If any convergence condition is met, terminate the iteration and output the final 3D coordinates: ,in Let be the actual number of iterations at termination, and satisfy . .

[0060] Step S160: Select and determine SAR control points based on the accuracy index of the final solution results.

[0061] The accuracy index of the final solution result is the in-point conformance accuracy, which is quantitatively evaluated based on the variance-covariance matrix obtained from the final 3D coordinate solution. Specifically, according to the type of 3D coordinate solution model used, the variance-covariance matrix of the candidate point at the final 3D geographic coordinates is calculated, and the mean square error of the candidate point in the planar and elevation directions is extracted from the diagonal elements of this matrix. The planar mean square error is represented by the arithmetic square root of the sum of the squares of the mean square errors in the north and east directions, and the elevation mean square error is directly composed of the mean square error of the celestial component. The in-planar conformance accuracy and the elevation conformance accuracy are defined as the numerical values ​​of the planar mean square error and the elevation mean square error, respectively, and are used to quantitatively describe the accuracy level of the solution result itself. This index achieves quality control by comparing with a preset threshold. When both the in-planar conformance accuracy and the elevation conformance accuracy are not greater than the preset threshold, the candidate point is determined to have sufficient positioning reliability and is thus identified as the final SAR control point; otherwise, it is discarded as an invalid point due to insufficient accuracy.

[0062] Optionally, step S160 may include: calculating the variance-covariance matrix at the final calculated three-dimensional coordinates according to the type of three-dimensional coordinate solution model used, and extracting the mean square errors of candidate points in the planar and elevation directions from the diagonal elements of the variance-covariance matrix; defining the intra-point conformance accuracy based on the mean square errors in the planar and elevation directions, wherein the intra-point conformance accuracy includes the intra-planar conformance accuracy and the intra-elevation conformance accuracy; determining whether both the intra-planar conformance accuracy and the intra-elevation conformance accuracy meet the preset threshold conditions; if they do, determining the candidate points as SAR control points and outputting them; if they do not meet the conditions, discarding the candidate points as invalid points.

[0063] In step S160 above, candidate control points are calculated based on the model type used in the solution. The variance-covariance matrix C at point X Based on C X The three-dimensional positional errors of the points in the plane and elevation directions are obtained respectively, and the in-place accuracy σ is defined. H and σ V If σ H and σ V Not greater than the preset threshold σ th If a candidate point is selected, it will be designated as the final SAR control point (SAR GCP); otherwise, it will be discarded. Specifically: From C X Extracting the mean square error in each direction from the diagonal elements: Planar mean square error (horizontal accuracy): ; Elevation error (vertical accuracy): ; If σ is satisfied H ≤σ th，H And σ v ≤σ th,v Typical value σ th,H =σ th,v If the accuracy is 0.5, the candidate point is determined as the final SAR control point; otherwise, it is considered insufficient and is automatically eliminated.

[0064] To facilitate understanding of the working principle of the SAR control point extraction method described above, a specific application example in a particular scenario is provided below. In this application scenario, the SAR control point extraction method mainly includes: Step 1: Based on the image coordinates, precise orbit parameters and image metadata of the matched SAR candidate points, establish 2N range-Doppler observation equations under N (N≥2) images to form a nonlinear observation model for multi-view three-dimensional positioning and form a multi-view geometric constraint system.

[0065] Step 2: Use the initial equal weight matrix P (0)Combined with the initial three-dimensional coordinates r of the candidate points t (0) The first three-dimensional coordinate solution of the candidate control points was completed, and the initial solution result r was obtained. t (1) .

[0066] If a linearized model is used, a weighted least squares solution is performed based on the linearized distance-Doppler equation; if a nonlinear model is used, the weighted residual objective function is minimized using the Levenberg–Marquardt (LMA) algorithm.

[0067] Step 3: Using the initial solution result r t (1) Calculate the current residual vector v lin (1) (Linear case) or v nl (1) (In the nonlinear case), perform geometrically constrained variance component estimation (GC-VCE) to update the variance components of each image observation: .

[0068] This method first constructs a prior variance model based on the intersection angle and orbit combination type. Then, it introduces a regularization term based on this prior model into the VCE objective function and corrects the effective degrees of freedom by incorporating the intersection angle. Finally, it solves for the optimized variance components to generate an adaptive weight matrix P. (1) .

[0069] Step 4: Based on the updated adaptive weight matrix P (1) The distance-Doppler model equations are solved again in three dimensions to obtain the second three-dimensional coordinate solution for the candidate control points, thus yielding the optimization result r. t (2) .

[0070] Step 5: Determine if the convergence condition is met: coordinate update amount ||r t (k+1) -r t (k) ||<ε or the number of iterations reaches the preset maximum value K max If convergence fails, then the current solution result r is used. t (k+1) Based on this, repeat steps 3 and 4, recalculate the residuals and update the weights, and continue iterative optimization. Otherwise, output the final 3D coordinates. (in , (This refers to the final iteration step number).

[0071] Step 6: Calculate candidate control points based on the model type used in the solution process. The variance-covariance matrix C at point X Based on C X The three-dimensional positional errors of the points in the plane and elevation directions are obtained respectively, and the in-place accuracy σ is defined. H and σ V If σ H and σ V Not greater than the preset threshold σ th If the candidate point is selected, it will be determined as the final SAR control point (SARGCP); otherwise, it will be discarded.

[0072] like Figure 3 As shown, based on the same inventive concept, this application also provides a SAR control point extraction device 200, comprising: The image acquisition module 210 is used to acquire multiple SAR images and the observation data and imaging geometric parameters of the matched candidate points in the multiple SAR images; The observation model building module 220 is used to build a multi-view positioning observation model based on data and parameters. The initial solution module 230 is used to perform initial three-dimensional coordinate solution on the observation model and obtain the initial solution results. The adaptive optimization module 240 is used to perform variance component estimation based on the initial solution results and combined with the prior constraints of the imaging geometry configuration, and dynamically optimize the observation weights of each image scene. The secondary solution module 250 is used to recalculate the three-dimensional coordinates based on the optimized observation weights; The control point selection module 260 is used to iteratively execute the weight optimization steps of the adaptive optimization module and the coordinate calculation steps of the secondary calculation module until the convergence condition is met; and to select and determine SAR control points based on the accuracy index of the final calculation result.

[0073] It is understood that the SAR control point extraction device 200 provided in this application embodiment can realize any of the functions of the above-mentioned SAR control point extraction method. For the method embodiment section, please refer to the method embodiment section for the manner and principle of realizing each function. The device embodiment section will not repeat the description.

[0074] Figure 4 A structural block diagram of an electronic device according to an embodiment of the present invention is shown. Figure 4 As shown, the electronic device includes a memory 310 and a processor 320. The memory 310 stores a computer program that can run on the processor 320. When the processor 320 executes the computer program, it implements the SAR control point extraction method in the above embodiments. The number of memories 310 and processors 320 can be one or more.

[0075] The electronic device / terminal / server also includes: The communication interface 330 is used to communicate with external devices and perform data exchange and transmission.

[0076] If the memory 310, processor 320, and communication interface 330 are implemented independently, they can be interconnected via a bus to communicate with each other. This bus can be an Industry Standard Architecture (ISA) bus, a Peripheral Component Interconnect (PCI) bus, or an Extended Industry Standard Architecture (EISA) bus, etc. This bus can be divided into address bus, data bus, control bus, etc. For ease of representation, Figure 4 The bus is represented by a single thick line, but this does not mean that there is only one bus or one type of bus.

[0077] Optionally, in a specific implementation, if the memory 310, processor 320 and communication interface 330 are integrated on a single chip, the memory 310, processor 320 and communication interface 330 can communicate with each other through an internal interface.

[0078] This application provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the method provided in this application.

[0079] This application also provides a chip, which includes a processor for calling and executing instructions stored in a memory, causing a communication device on which the chip is installed to perform the method provided in this application.

[0080] This application also provides a chip, including: an input interface, an output interface, a processor, and a memory. The input interface, output interface, processor, and memory are connected through an internal connection path. The processor is used to execute code in the memory. When the code is executed, the processor is used to execute the method provided in this application.

[0081] It should be understood that the aforementioned processor can be a Central Processing Unit (CPU), or other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. General-purpose processors can be microprocessors or any conventional processor. It is worth noting that the processor can be a processor supporting the Advanced Reduced Instruction Set Computing (RISC) machine (ARM) architecture.

[0082] Further, optionally, the aforementioned memory may include read-only memory and random access memory, and may also include non-volatile random access memory. The memory may be volatile or non-volatile, or may include both. Non-volatile memory may include read-only memory (ROM), programmable read-only memory (PROM), erasable programmable read-only memory (EPROM), electrically erasable programmable read-only memory (EEPROM), or flash memory. Volatile memory may include random access memory (RAM), which serves as an external cache. Many forms of RAM are available by way of example, but not limitation. Examples include static random access memory (SRAM), dynamic random access memory (DRAM), synchronous dynamic random access memory (SDRAM), double data rate synchronous dynamic random access memory (DDR SDRAM), enhanced synchronous dynamic random access memory (ESDRAM), synchronous linked dynamic random access memory (SLDRAM), and direct rambus RAM (DR RAM).

[0083] In the above embodiments, implementation can be achieved, in whole or in part, through software, hardware, firmware, or any combination thereof. When implemented in software, it can be implemented, in whole or in part, as a computer program product. A computer program product includes one or more computer instructions. When the computer program instructions are loaded and executed on a computer, all or part of the flow or function according to this application is generated. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. The computer instructions can be stored in a computer-readable storage medium or transferred from one computer-readable storage medium to another.

[0084] Furthermore, the units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment according to actual needs.

[0085] Furthermore, the functional modules in the various embodiments of the present invention can be integrated together to form an independent part, or each module can exist independently, or two or more modules can be integrated to form an independent part.

[0086] The above description is merely an embodiment of the present invention and is not intended to limit the scope of protection of the present invention. For those skilled in the art, the present invention can have various modifications and variations. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for extracting SAR control points, characterized in that, The method includes: Acquire multiple SAR images and the observation data and imaging geometric parameters of the matched candidate points in the multiple SAR images; A multi-view positioning and observation model is established based on the aforementioned data and parameters; The initial three-dimensional coordinates of the observation model are calculated to obtain the initial calculation results; Based on the initial solution results, variance component estimation is performed in conjunction with the prior constraints of the imaging geometry configuration to dynamically optimize the observation weights of each image scene. Based on the optimized observation weights, the three-dimensional coordinates are recalculated. Iteratively execute the weight optimization step and the coordinate calculation step until the convergence condition is met; SAR control points are selected and determined based on the accuracy index of the final solution results; The prior constraints combined with imaging geometry configuration are used to perform variance component estimation, including: For each image, the intersection angle is calculated based on the spatial relationship between the satellite position vector of the candidate point at the imaging time of the image and the satellite position vectors of the other images at the corresponding times. The geometric robustness level of the stereo image pair formed by the image and the other images is determined according to the intersection angle, and the highest level is selected as the comprehensive geometric robustness level of the image. The orbit combination type between images is determined based on the orbital affiliation of each image. An orbital system error penalty factor is introduced for the combination of images on the same orbit to suppress the influence of the correlation of orbital system errors on the solution results. A prior variance model is constructed based on the comprehensive geometric robustness level and the error penalty factor of the same orbital system, and the prior variance model is introduced as a regularization term into the objective function of variance component estimation.

2. The SAR control point extraction method according to claim 1, characterized in that, The objective function for the variance component estimation includes an observation residual term and a regularization term based on the prior variance model; The dynamic optimization of the observation weights of each image scene includes: obtaining the optimized variance components of each image scene by minimizing the objective function, and generating an adaptive weight matrix based on the optimized variance components to dynamically optimize the observation weights of each image scene.

3. The SAR control point extraction method according to claim 1, characterized in that, The process of selecting and determining SAR control points based on the accuracy index of the final solution results includes: Based on the type of three-dimensional coordinate solution model used, the variance-covariance matrix at the final three-dimensional coordinates is calculated, and the mean square error of the candidate point in the plane and elevation directions is extracted from the diagonal elements of the variance-covariance matrix. The point-positional accuracy is defined based on the mean square error in the plane direction and the mean square error in the elevation direction. The point-positional accuracy includes the plane-positional accuracy and the elevation-positional accuracy. Determine whether the in-plane conformance accuracy and in-elevation conformance accuracy both meet the preset threshold conditions. If they do, the candidate point is determined as a SAR control point and output. If they do not meet the threshold conditions, the candidate point is discarded as an invalid point.

4. The SAR control point extraction method according to any one of claims 1 to 3, characterized in that, The initial three-dimensional coordinate calculation of the observation model includes: Based on the multi-view positioning and observation model, the initial solution is performed by combining the initial equal weight matrix with the initial three-dimensional coordinates of the candidate points. The initial solution is performed using either a linear model or a nonlinear model: If a linearized model is adopted, a first-order Taylor expansion is performed on the observation model at the initial three-dimensional coordinates to construct the design matrix and constant terms, and the weighted least squares problem is solved to obtain the coordinate correction, so as to update the initial solution result. If a nonlinear model is adopted, the initial three-dimensional coordinates are used as initial values, and the weighted residual objective function is minimized by an iterative optimization algorithm with damping factors. The iterative optimization process includes calculating the residual vector and partial derivative matrix, constructing the damping normal equation and solving the coordinate correction. After multiple rounds of iterative optimization, the initial solution result is output.

5. The SAR control point extraction method according to any one of claims 1 to 3, characterized in that, The convergence condition includes at least one of the following: The coordinate update amount between the coordinate solution result of the current iteration step and the coordinate solution result of the previous iteration step is less than the preset convergence threshold; The current iteration step has reached the preset maximum number of iterations.

6. The SAR control point extraction method according to any one of claims 1 to 3, characterized in that, The establishment of a multi-view positioning and observation model based on the data and parameters includes: For each image, based on the image coordinates of the candidate points in the image and the satellite position vector and velocity vector corresponding to the image imaging time, the distance observation equation and the Doppler observation equation are constructed respectively. The distance observation equations of each image are combined with the Doppler observation equations to form a set of nonlinear equations with the three-dimensional geographic coordinates of the candidate points as unknowns, which serves as the multi-view positioning observation model.

7. A SAR control point extraction device, characterized in that, include: The image acquisition module is used to acquire multiple SAR images and the observation data and imaging geometric parameters of the matched candidate points in the multiple SAR images; The observation model building module is used to build a multi-view positioning observation model based on the data and parameters. The initial solution module is used to perform initial three-dimensional coordinate solution on the observation model and obtain the initial solution result; The adaptive optimization module is used to perform variance component estimation based on the initial solution results and in combination with the prior constraints of the imaging geometry configuration, and dynamically optimize the observation weights of each image scene. The secondary solution module is used to recalculate the three-dimensional coordinates based on the optimized observation weights; The control point selection module is used to iteratively execute the weight optimization step of the adaptive optimization module and the coordinate calculation step of the secondary solution module until the convergence condition is met; and to select and determine SAR control points based on the accuracy index of the final solution result. The adaptive optimization module is specifically used for: for each image, calculating the intersection angle based on the spatial relationship between the satellite position vector of the candidate point at the image imaging time and the satellite position vectors of the other images at the corresponding times, and determining the geometric robustness level of the stereo image pair formed by the image and the other images according to the intersection angle, and selecting the highest level as the comprehensive geometric robustness level of the image; determining the orbit combination type between images according to the orbital affiliation relationship of each image, introducing a co-orbital system error penalty factor for co-orbital image combinations to suppress the influence of co-orbital system error correlation on the solution results; constructing a prior variance model based on the comprehensive geometric robustness level and the co-orbital system error penalty factor, and introducing the prior variance model as a regularization term into the objective function of variance component estimation.

8. An electronic device, characterized in that, include: A processor and a memory, wherein instructions are stored in the memory and loaded and executed by the processor to implement the method as claimed in any one of claims 1-6.

9. A computer-readable storage medium storing a computer program that, when executed by a processor, implements the method as described in any one of claims 1-6.

Images