A geometric calibration method for spaceborne SAR based on control points provided by wind turbines

Abstract

The present invention relates to a spaceborne SAR geometric calibration method based on control points provided by wind turbines, comprising: acquiring SAR images; constructing an initial geometric positioning model to realize mutual mapping between SAR images and DOM based on image coordinates and object coordinates; collecting the same-name point pairs of wind turbines in SAR images and DOM, and dividing them into control point pairs and check point pairs; calculating the geometric positioning error point by point for point data in the control point database; calculating the atmospheric delay value; analyzing the source of the geometric positioning error, and establishing a geometric calibration model considering the atmospheric delay, which is used to solve the error correction parameters of the SAR satellite imaging system; solving the systematic error correction parameters based on the control point database and the geometric calibration model, and updating the initial geometric positioning model; and evaluating the positioning performance of the model. Compared with the prior art, the present invention has the advantages of not requiring the support of a geometric calibration field and having strong versatility.

Description

Technical Field

[0001] The invention relates to a unified geometric calibration method for SAR images, and in particular to a spaceborne SAR geometric calibration method based on control points provided by wind turbines. Background Art

[0002] Synthetic Aperture Radar (SAR) satellites are inevitably affected by factors such as platform stability, equipment aging, observation environment changes, and orbit control during their on-orbit operation. Laboratory calibration parameters cannot meet the requirements of high-precision mapping tasks during the satellite's life cycle. Eliminating the system errors during the on-orbit operation of SAR satellites and improving the positioning accuracy of SAR images are important issues in the field of radar photogrammetry, and play an important role in exploring the application potential of SAR images.

[0003] Regular geometric calibration helps to eliminate the accumulated systematic errors during the life cycle of the satellite, and is the most economical and effective means to improve the accuracy of the uncontrolled geometric positioning of spaceborne SAR. At present, the mainstream SAR geometric calibration models are established using strict geometric positioning models based on the principle of SAR imaging. All publicly published work relies on atmospheric parameter products (European Center for Medium-Term Numerical Weather Forecasts ECMWF and National Center for Environmental Prediction NCEP) or measured atmospheric parameters. Weather numerical forecast products often have unpredictable results. The data generated by the reanalysis system often comes from the assimilation of multi-source meteorological data, often with a time delay of several months. Field measurement of atmospheric parameters requires a lot of manpower and material resources, which is almost impossible to achieve in large-scale geometric calibration tasks. Another dependence of geometric calibration on external data is reflected in the dependence on high-precision ground control points, which often requires researchers to use GNSS equipment on site for high-precision measurements, which is extremely difficult for large-scale geometric calibration tasks. Considering that the positioning accuracy of domestic civilian SAR satellites is only better than 50 meters, for most mapping tasks, it is still necessary to further improve its geometric performance in order to give full play to the advantages of SAR images in all-day and all-weather imaging. Summary of the invention

[0004] The purpose of the present invention is to provide a spaceborne SAR geometric calibration method based on control points provided by wind turbines, which does not require the support of a geometric calibration field and can achieve geometric calibration without relying on reanalysis of meteorological products, thereby improving the ground geometric positioning performance of SAR images.

[0005] The purpose of the present invention can be achieved by the following technical solutions:

[0006] A spaceborne SAR geometric calibration method based on control points provided by wind turbines comprises the following steps:

[0007] Step 1) Acquire SAR images;

[0008] Step 2) constructing an initial geometric positioning model based on the auxiliary file of the SAR image to achieve mutual mapping between the SAR image and the DOM based on the image space coordinates and the object space coordinates;

[0009] Step 3) After preprocessing the SAR image, collect the same-name point pairs of the wind turbine in the SAR image and DOM, and record the image coordinates of the target point in the SAR image and the corresponding object coordinates in the DOM respectively, obtain the elevation information in the DEM according to the object coordinates, generate ground control points of the same-name point pairs, obtain multiple pairs of same-name point pairs through multiple acquisitions, and divide them into control point pairs and check point pairs, and establish a control point pair database and a check point pair database;

[0010] Step 4) Calculate the geometric positioning error of the point data in the check point database point by point;

[0011] Step 5) calculating atmospheric parameters based on a standard atmospheric model, estimating tropospheric atmospheric delay based on the Sastamonen model and the atmospheric parameters, and estimating ionospheric delay based on the Krobucher model;

[0012] Step 6) analyzing the source of geometric positioning error, and establishing a geometric calibration model taking into account atmospheric delay, for solving the error correction parameters of the SAR satellite imaging system, wherein the atmospheric delay value of the geometric calibration model is determined based on the tropospheric atmospheric delay and the ionospheric delay;

[0013] Step 7) calculating error correction parameters based on the control point database and the geometric calibration model, and updating the initial geometric positioning model based on the error correction parameters to obtain an updated geometric positioning model;

[0014] Step 8) Evaluate the positioning performance of the initial geometric positioning model and the updated geometric positioning model.

[0015] The step 2) comprises the following steps:

[0016] Step 2-1) Fit the position and velocity during satellite imaging based on the auxiliary data in XML format attached to the first-level SAR image product:

[0017]

[0018] Where (a0, a1, a2, a3, b0, b1, b2, b3, c0, c1, c2, c3) are the coefficients of the fitted position equation; t a is the imaging time corresponding to the j-th row of pixels in the SAR image, i.e., the azimuth imaging time; j is the azimuth image coordinate, i.e., the row coordinate of the observed target in the SAR image; t a0 is the starting imaging time in azimuth, f ais the pulse repetition frequency, (X g ,Y g ,Z g ) is the object coordinate of the ground target point, (X s ,Y s ,Z s ) is the spatial rectangular coordinate of the sensor at the moment of imaging the ground target point, (V x ,V y ,V z ) is the velocity vector of the pixel in the jth row, and R is the slant range measurement between the observed target and the sensor;

[0019] Step 2-2) Based on the SAR satellite imaging principle, construct the initial geometric positioning model for earth observation:

[0020]

[0021] Where R0 is the starting slant distance in the range direction, i is the image coordinate in the range direction, i.e., the column coordinate of the observed target in the SAR image, and P r is the range-slant resolution; f d is the Doppler center frequency between the satellite platform and the observed target at the imaging time, P sc and V sc are the position and velocity vectors of the sensor, P gc and V gc are the position and velocity vector of the observed target, respectively, λ is the radar wavelength; a and b are the major and minor semi-axes of the WGS-84 reference ellipsoid, respectively, and h is the average altitude of the area where the earth observation target is located;

[0022] Step 2-3) Determine the conditional equation of the geometric positioning model based on the error:

[0023]

[0024] Where, t r0 is the distance to the starting time, c is the speed of light, f r is the range sampling frequency of the radar system; t r is the imaging time in the range direction;

[0025] Step 2-4) For each ground target point, the conditional equation of the geometric positioning model is established, and the object space coordinates of the center point provided by the SAR image product are used as the initial value. Based on the least squares principle, the Newton iteration method is used to solve the object space coordinates of the ground target point. During the iteration, the partial derivatives of the object space coordinates of the ground target point are obtained according to the conditional equations of the geometric positioning model to obtain the partial derivative matrix, and the partial derivative matrix is ​​updated. When the correction value of the object space coordinates relative to the previous iteration is less than the pre-configured threshold, the iteration is terminated to obtain the object space coordinates of the ground target point, and the mutual mapping of the image space coordinates and the object space coordinates is completed.

[0026] The step 3) comprises the following steps:

[0027] Step 3-1) Preprocess the SAR image: collect the first-level products of multi-view SAR images in the study area, perform radiation correction, and convert the backscatter information into intensity information;

[0028] Step 3-2) Select a wind turbine target point in a flat area with simple surrounding features in the SAR image, and record its image coordinates;

[0029] Step 3-3) using the initial geometric positioning model constructed in step 2), mapping the image space coordinates recorded in step 3-2) to object space coordinates, narrowing the search space, quickly locating the same-name points in the DOM, and recording the object space coordinates of the same-name points;

[0030] Step 3-4) Obtain elevation information in DEM based on the object coordinates of the same-name points and generate ground control points for the same-name point pairs;

[0031] Step 3-4) Repeat steps 3-2) to 3-4) to complete the collection of same-name point pairs, and divide the same-name point pairs into control point pairs and check point pairs, and establish a control point pair database and a check point pair database.

[0032] The step 4) is specifically as follows:

[0033] The positioning results of the check point pairs are calculated based on the initial geometric positioning model of each image and compared with the object coordinates of the check point pairs. The difference between the two is the geometric positioning error of the initial geometric positioning model.

[0034] The step 5) comprises the following steps:

[0035] Step 5-1) Estimate the atmospheric parameters at any location and time based on the standard atmospheric model:

[0036]

[0037] Where Press is atmospheric pressure, Temp is Kelvin temperature, h is altitude, w is relative humidity, W pressis the water vapor partial pressure;

[0038] Step 5-2) Based on the atmospheric parameters calculated in step 5-1), the tropospheric atmospheric path delay is calculated using the Sastamonen model as follows:

[0039]

[0040] Where θ is the SAR satellite incident angle, φ is the satellite elevation angle, and D h is the dry delay, D w is the wet delay, ZTD is the tropospheric atmospheric delay;

[0041] Step 5-3) Estimate the ionospheric path delay based on the Krobucher model:

[0042]

[0043] Where TEC is the ionospheric electron concentration; f is the radar center frequency; ZID is the ionospheric delay;

[0044] Step 5-4) Determine the atmospheric delay value based on the tropospheric atmospheric delay and the ionospheric delay:

[0045] t delay =ZTD+ZID

[0046] The step 6) is specifically as follows:

[0047] Step 6-1) According to the constructed initial geometric positioning model, pre-configured errors are added to the imaging parameters respectively, wherein the imaging parameters include: the orbital coordinates (X s ,Y s ,Z s ), starting slant distance R0 in range, slant distance resolution P in range r , system distance sampling frequency f r and the pulse repetition frequency f a ;

[0048] Step 6-2) According to the working principle of the SAR imaging system, the sources of the SAR image geometric positioning error include the electronic delay inside the system, the time synchronization error between the SAR payload and the GNSS payload, the radar sampling frequency and the pulse repetition frequency error;

[0049] Step 6-3) Determine the SAR geometric calibration model based on the geometric positioning error sources:

[0050]

[0051] Where Δf r and Δf aare the correction values ​​of the range sampling frequency and pulse repetition frequency respectively; t r ,t a is the imaging time in range and azimuth, t r0 ,t a0 is the starting time in range and azimuth, t delay is the atmospheric delay value, Δt r is the electronic delay inside the system, Δt a is the synchronization error between the SAR payload and the GNSS payload;

[0052] Step 6-4) Perform SAR image-to-ground geometric positioning error analysis based on the added pre-configured error, and ignore Δf according to the analysis results. r and Δf a The impact on geometric positioning results in a simplified SAR geometric calibration model:

[0053]

[0054] Step 6-5) Determine the conditional equation of the geometric calibration model based on the SAR imaging geometric characteristics:

[0055]

[0056] The step 7) comprises the following steps:

[0057] Step 7-1) uses Taylor series to expand the calibration parameters in the conditional equation of step 6-5) to the first-order term to obtain the Jacobian matrix:

[0058]

[0059] Where X is the error correction parameter;

[0060] Step 7-2) obtaining a preconfigured number of control point pairs in a control point pair database;

[0061] Step 7-3) establishing a conditional equation of a geometric calibration model for each control point pair;

[0062] Step 7-4) Based on the least squares principle, the Newton iteration method is used to update the Jacobian matrix based on the conditional equation of the geometric calibration model. When X is less than a pre-configured threshold, the iteration process ends and the most likely estimation of the error correction parameter is completed;

[0063] Step 7-5) Update the initial geometric positioning model of the SAR image based on the error correction parameters to obtain an updated geometric positioning model.

[0064] The step 8) is specifically as follows:

[0065] The initial geometric positioning model and the updated geometric positioning model are respectively evaluated for the positioning deviation of each check point pair using the check point pairs. The statistical indicators of all the check point pairs contained in each scene image are counted to compare the performance difference between the two.

[0066] The statistical indicators include minimum deviation, maximum deviation, mean value, and root mean square error.

[0067] In step 3-4), the point pairs with the same name are divided into control point pairs and check point pairs according to the spatial distribution of the control points, wherein the number of control point pairs in each scene image is not less than 4 and the number of check point pairs is not less than 4.

[0068] Compared with the prior art, the present invention has the following beneficial effects:

[0069] (1) The present invention uses wind turbines to provide ground control point information, and can achieve geometric calibration without the support of geometric calibration sites and GNSS equipment. It has strong versatility and a wide range of applications.

[0070] (2) The present invention adopts a standard atmospheric model to obtain atmospheric parameters and combines the Sastamonen model to estimate the tropospheric atmospheric delay and adopts the Krobucher model to estimate the ionospheric delay, thereby realizing near real-time estimation of atmospheric path delay, getting rid of the dependence on ground geometry calibration fields and reanalysis meteorological data, and having high computational efficiency, and can support the geometric calibration method of satellite-borne SAR images based on ground control point information provided by wind turbines.

[0071] (3) The present invention determines the geometric calibration model by analyzing the error sources, estimates the error correction parameters of the SAR satellite imaging system, and updates the geometric positioning model. The positioning accuracy of the updated geometric positioning model is greatly improved compared with the initial geometric positioning model, thereby improving the geometric positioning performance of domestic SAR satellite images. BRIEF DESCRIPTION OF THE DRAWINGS

[0072] Figure 1 is a flow chart of the method of the present invention;

[0073] Figure 2 Schematic diagram of the imaging process of wind turbines in SAR images;

[0074] Figure 3 It is a schematic diagram of the overall design and execution process of the present invention;

[0075] Figure 4 A schematic diagram of a research area of ​​an embodiment;

[0076] Figure 5Characteristic images of wind turbines in optical images and SAR images, where (a), (b), (c), (g), (h), and (i) are characteristic images in optical images, and (d), (e), (f), (j), (k), and (l) are characteristic images in SAR images corresponding to the optical images. DETAILED DESCRIPTION

[0077] The present invention is described in detail below in conjunction with the accompanying drawings and specific embodiments. This embodiment is implemented based on the technical solution of the present invention, and provides a detailed implementation method and specific operation process, but the protection scope of the present invention is not limited to the following embodiments.

[0078] Before establishing a SAR image geometric calibration method based on ground control point information provided by wind turbines, it is necessary to verify the feasibility of using wind turbines as ground control points in SAR geometric calibration through theoretical analysis. The specific steps include:

[0079] a) Analyze the SAR imaging process and SAR image product processing flow, where the imaging process of wind turbines in SAR images is as follows: Figure 2 As shown in the figure, compared with the traditional trihedral corner reflector, the wind turbine has the advantages of large radar backscattering cross-section and wide spatial distribution, which can fully tap its application potential in the field of SAR geometric calibration and is expected to build a global geometric calibration field;

[0080] b) Use Solid Works to build a 3D geometric model of a common wind turbine and analyze the transmission path of the microwave signal emitted by the SAR sensor;

[0081] c) Analyze and establish the correspondence between the reflection center of the wind turbine and the features in the image, use the position of the center of the wind turbine base in the DOM as the ground control point, obtain the elevation of this position in the DEM, and determine the location and method of wind turbine coordinate collection.

[0082] The analysis described in c) to establish the correspondence between the reflection center of the wind turbine and the features in the image mainly includes: analyzing the geometric structure of the wind turbine and analyzing its backscattering intensity. Wind turbines are generally installed perpendicular to the ground in locations with flat terrain and relatively simple surrounding landforms and environmental elements. The superior installation environment provides a lower background backscattering coefficient. The wind turbine and the ground form an equivalent huge dihedral reflector with strong backscattering capability. Therefore, the microwave signal actively transmitted by the SAR satellite can receive strong echo signals in two forms: satellite SAR antenna-ground-wind turbine tower-satellite SAR antenna or satellite SAR antenna-wind turbine tower-ground-satellite SAR antenna, thereby forming a "cross" bright spot in the SAR image, which is sharply different from the surrounding environment.

[0083] Based on the above analysis, the feasibility of using wind turbines as ground control points is determined, and then a spaceborne SAR geometric calibration method based on wind turbines providing control points is proposed. Figure 1 As shown, the following steps are included:

[0084] Step 1) Acquire SAR images;

[0085] In this embodiment, the regional location is as follows Figure 4 , the Yangtze River Delta of China was used as the implementation area, and the Gaofen-3 SAR images from October to November 2018 were used as the research objects. The basic information of the acquired SAR images is shown in Table 1.

[0086] Table 1 Basic information of SAR images

[0087]

[0088] Step 2) constructing an initial geometric positioning model based on the auxiliary file of the SAR image to achieve mutual mapping between the SAR image and the DOM based on the image space coordinates and the object space coordinates;

[0089] Step 2-1) Use a cubic polynomial to fit the position and velocity of the satellite during imaging based on the auxiliary data in XML format attached to the first-level SAR image product:

[0090]

[0091] Where (a0, a1, a2, a3, b0, b1, b2, b3, c0, c1, c2, c3) are the coefficients of the fitted position equation; t a is the imaging time corresponding to the j-th row of pixels in the SAR image, i.e., the azimuth imaging time; j is the azimuth image coordinate, i.e., the row coordinate of the observed target in the SAR image; t a0 is the starting imaging time in azimuth, f a is the pulse repetition frequency, (X g ,Y g ,Z g ) is the object coordinate of the ground target point, (X s ,Y s ,Z s ) is the spatial rectangular coordinate of the sensor at the moment of imaging the ground target point, (V x ,V y ,V z ) is the velocity vector of the pixel in the jth row, and R is the slant range measurement between the observed target and the sensor;

[0092] Step 2-2) Based on the SAR satellite imaging principle, construct the initial geometric positioning model for earth observation:

[0093]

[0094] Where R0 is the starting slant distance in the range direction, i is the image coordinate in the range direction, i.e., the column coordinate of the observed target in the SAR image, and P r is the range-slant resolution; f d is the Doppler center frequency between the satellite platform and the observed target at the imaging time, P sc and V sc are the position and velocity vectors of the sensor, P gc and V gc are the position and velocity vector of the observed target, respectively, λ is the radar wavelength; a and b are the major and minor semi-axes of the WGS-84 reference ellipsoid, respectively, and h is the average altitude of the area where the earth observation target is located;

[0095] Step 2-3) Determine the conditional equation of the geometric positioning model based on the error:

[0096]

[0097] Where, t r0 is the distance to the starting time, c is the speed of light, f r is the range sampling frequency of the radar system; t r is the imaging time in the range direction;

[0098] Step 2-4) For each ground target point, the conditional equation of the geometric positioning model is established, and the object space coordinates of the center point provided by the SAR image product are used as the initial value. Based on the least squares principle, the Newton iteration method is used to solve the object space coordinates of the ground target point. During the iteration, the partial derivatives of the object space coordinates of the ground target point are obtained according to the conditional equations of the geometric positioning model to obtain the partial derivative matrix, and the partial derivative matrix is ​​updated. When the correction value of the object space coordinates relative to the previous iteration is less than the pre-configured threshold, the iteration is terminated to obtain the object space coordinates of the ground target point, and the mutual mapping of the image space coordinates and the object space coordinates is completed.

[0099] Step 3) After preprocessing the SAR image, collect the same-name point pairs of the wind turbine in the SAR image and DOM, and record the image coordinates of the target point in the SAR image and the corresponding object coordinates in the DOM respectively, obtain the elevation information in the DEM according to the object coordinates, generate ground control points of the same-name point pairs, obtain multiple pairs of same-name point pairs through multiple acquisitions, and divide them into control point pairs and check point pairs, and establish a control point pair database and a check point pair database;

[0100] Step 3-1) Preprocess the SAR image: collect the first-level products of multi-view SAR images in the study area, perform radiation correction, and convert the backscatter information into intensity information;

[0101] Step 3-2) Select a wind turbine target point in a flat area with simple surrounding features in the SAR image, and record its image coordinates;

[0102] Step 3-3) using the initial geometric positioning model constructed in step 2), mapping the image space coordinates recorded in step 3-2) to object space coordinates, locating the same-name points in the DOM, and recording the object space coordinates of the same-name points;

[0103] Step 3-4) Obtain elevation information in DEM based on the object coordinates of the same-name points and generate ground control points for the same-name point pairs;

[0104] Step 3-4) Repeat steps 3-2) - 3-4) to complete the collection of same-name point pairs, and divide the same-name point pairs into control point pairs and check point pairs. According to the spatial distribution of control point pairs, ensure that each scene image has no less than 4 control point pairs and 4 check point pairs, and establish a control point pair database and a check point pair database.

[0105] In this embodiment, the characteristic images of the wind turbine in the optical image and the SAR image are as follows: Figure 5 shown.

[0106] Step 4) Calculate the geometric positioning error of the point data in the check point database point by point;

[0107] The positioning results of the check point pairs are calculated based on the initial geometric positioning model of each image and compared with the object coordinates of the check point pairs. The difference between the two is the geometric positioning error of the initial geometric positioning model.

[0108] Step 5) calculating atmospheric parameters based on a standard atmospheric model, estimating tropospheric atmospheric delay based on the Sastamonen model and the atmospheric parameters, and estimating ionospheric delay based on the Krobucher model;

[0109] Step 5-1) Estimate the atmospheric parameters at any location and time based on the standard atmospheric model:

[0110]

[0111] Where Press is atmospheric pressure, Temp is Kelvin temperature, h is altitude, w is relative humidity, W press is the water vapor partial pressure;

[0112] Step 5-2) Based on the atmospheric parameters calculated in step 5-1), the tropospheric atmospheric path delay is calculated using the Sastamonen model as follows:

[0113]

[0114] Where θ is the SAR satellite incident angle, φ is the satellite elevation angle, and D h is the dry delay, Dw is the wet delay, ZTD is the tropospheric atmospheric delay;

[0115] Step 5-3) Estimate the ionospheric path delay based on the Krobucher model:

[0116]

[0117] Where TEC is the ionospheric electron concentration, which can be provided by the daily 2-hour ionospheric electron concentration product provided by the Center for Orbit Determination in Europe (CODE); f is the radar center frequency; ZID is the ionospheric delay;

[0118] Step 5-4) Determine the atmospheric delay value based on the tropospheric atmospheric delay and the ionospheric delay:

[0119] t delay =ZTD+ZID

[0120] In this embodiment, the difference between the tropospheric atmospheric delay estimated by the Sastamonen model and the reference value provided by the International GNSS Service (IGS) is shown in Table 2.

[0121] Table 2 Comparison of tropospheric delay estimation bias

[0122]

[0123]

[0124] Step 6) analyzing the source of geometric positioning error, and establishing a geometric calibration model taking into account atmospheric delay, for solving the error correction parameters of the SAR satellite imaging system, wherein the atmospheric delay value of the geometric calibration model is determined based on the tropospheric atmospheric delay and the ionospheric delay;

[0125] Step 6-1) According to the constructed initial geometric positioning model, pre-configured errors are added to the imaging parameters respectively, wherein the imaging parameters include: the orbital coordinates (X s ,Y s ,Z s ), starting slant distance R0 in range, slant distance resolution P in range r , system distance sampling frequency f r and the pulse repetition frequency f a ;

[0126] Step 6-2) According to the working principle of the SAR imaging system, the sources of the SAR image geometric positioning error include the electronic delay inside the system, the time synchronization error between the SAR payload and the GNSS payload, the radar sampling frequency and the pulse repetition frequency error;

[0127] Step 6-3) Determine the SAR geometric calibration model based on the error sources:

[0128]

[0129] Where Δf r and Δf a are the correction values ​​of the range sampling frequency and pulse repetition frequency respectively; t r ,t a is the imaging time in range and azimuth, t r0 ,t a0 is the starting time in range and azimuth, t delay is the atmospheric delay value, Δt r is the electronic delay inside the system, Δt a is the synchronization error between the SAR payload and the GNSS payload;

[0130] Step 6-4) At present, the radar sampling frequency, sensor crystal oscillator and other components can achieve high accuracy. Based on the added pre-configured error, the SAR image to ground geometric positioning error is analyzed. According to the analysis results, Δf r and Δf a The effect on geometric positioning is small, so Δf is ignored. r and Δf a The impact on geometric positioning results in a simplified SAR geometric calibration model:

[0131]

[0132] Step 6-5) Determine the conditional equation of the geometric calibration model based on the SAR imaging geometric characteristics:

[0133]

[0134] Step 7) calculating error correction parameters based on the control point database and the geometric calibration model, and updating the initial geometric positioning model based on the error correction parameters to obtain an updated geometric positioning model;

[0135] Step 7-1) uses Taylor series to expand the calibration parameters in the conditional equation of step 6-5) to the first-order term to obtain the Jacobian matrix:

[0136]

[0137] Where X is the error correction parameter;

[0138] Step 7-2) obtaining a preconfigured number of control point pairs in a control point pair database;

[0139] Step 7-3) establishing a conditional equation of a geometric calibration model for each control point pair;

[0140] Step 7-4) Based on the least squares principle, the Newton iteration method is used to update the Jacobian matrix based on the conditional equation of the geometric calibration model. When X is less than a pre-configured threshold, the iteration process ends and the most likely estimation of the error correction parameter is completed;

[0141] Step 7-5) Update the initial geometric positioning model of the SAR image based on the error correction parameters to obtain an updated geometric positioning model.

[0142] Step 8) Evaluate the positioning performance of the initial geometric positioning model and the updated geometric positioning model.

[0143] The initial geometric positioning model and the updated geometric positioning model are respectively evaluated for the positioning deviation at each check point pair using the check point pair, and the statistical indicators of all the check point pairs contained in each scene image are counted to compare the performance difference between the two, wherein the statistical indicators include the minimum deviation, the maximum deviation, the mean deviation, and the root mean square error. In this embodiment, the positioning results obtained by the geometric calibration method (joint calibration strategy) using the wind turbine to provide ground control point information and the positioning results obtained by the initial RD model (single scene strategy) are shown in Table 3.

[0144] Based on the above process, the overall design and implementation process of the present invention is as follows: Figure 3 shown.

[0145] Table 3 Geometric calibration results of two calibration strategies

[0146]

[0147] The preferred specific embodiments of the present invention are described in detail above. It should be understood that a person skilled in the art can make many modifications and changes based on the concept of the present invention without creative work. Therefore, any technical solution that can be obtained by a person skilled in the art through logical analysis, reasoning, or limited experiments based on the concept of the present invention on the basis of the prior art should be within the scope of protection determined by the claims.

Claims

1. A spaceborne SAR geometric calibration method based on control points provided by wind turbines, characterized in that: The following steps are involved: Step 1) Acquire SAR images; Step 2) constructing an initial geometric positioning model based on the auxiliary file of the SAR image to achieve mutual mapping between the SAR image and the DOM based on the image space coordinates and the object space coordinates; Step 3) After preprocessing the SAR image, collect the same-name point pairs of the wind turbine in the SAR image and DOM, and record the image coordinates of the target point in the SAR image and the corresponding object coordinates in the DOM respectively, obtain the elevation information in the DEM according to the object coordinates, generate ground control points of the same-name point pairs, obtain multiple pairs of same-name point pairs through multiple acquisitions, and divide them into control point pairs and check point pairs, and establish a control point pair database and a check point pair database; Step 4) Calculate the geometric positioning error of the point data in the check point database point by point; Step 5) calculating atmospheric parameters based on a standard atmospheric model, estimating tropospheric atmospheric delay based on the Sastamonen model and the atmospheric parameters, and estimating ionospheric delay based on the Krobucher model; Step 6) analyzing the source of geometric positioning error, and establishing a geometric calibration model taking into account atmospheric delay, for solving the error correction parameters of the SAR satellite imaging system, wherein the atmospheric delay value of the geometric calibration model is determined based on the tropospheric atmospheric delay and the ionospheric delay; Step 7) calculating error correction parameters based on the control point database and the geometric calibration model, and updating the initial geometric positioning model based on the error correction parameters to obtain an updated geometric positioning model; Step 8) evaluating the positioning performance of the initial geometric positioning model and the updated geometric positioning model based on the check point pair database; The step 2) comprises the following steps: Step 2-1) Fit the position and velocity during satellite imaging based on the auxiliary data in XML format attached to the first-level SAR image product: Where (a0, a1, a2, a3, b0, b1, b2, b3, c0, c1, c2, c3) are the coefficients of the fitted position equation; t a is the imaging time corresponding to the j-th row of pixels in the SAR image, i.e., the azimuth imaging time; j is the azimuth image coordinate, i.e., the row coordinate of the observed target in the SAR image; t a0 is the starting imaging time in azimuth, f a is the pulse repetition frequency, (X g ,Y g ,Z g ) is the object coordinate of the ground target point, (X s ,Y s ,Z s ) is the spatial rectangular coordinate of the sensor at the moment of imaging the ground target point, (V x ,V y ,V z ) is the velocity vector of the pixel in the jth row, and R is the slant range measurement between the observed target and the sensor; Step 2-2) Based on the SAR satellite imaging principle, construct the initial geometric positioning model for earth observation: Where R0 is the starting slant distance in the range direction, i is the image coordinate in the range direction, i.e., the column coordinate of the observed target in the SAR image, and P r is the range-slant resolution; f d is the Doppler center frequency between the satellite platform and the observed target at the imaging time, P sc and V sc are the position and velocity vectors of the sensor, P gc and V gc are the position and velocity vector of the observed target, respectively; λ is the radar wavelength; a and b are the major and minor semi-axes of the WGS-84 reference ellipsoid, respectively; h is the average altitude of the area where the earth observation target is located; Step 2-3) Determine the conditional equation of the geometric positioning model based on the error: Where, t r0 is the distance to the starting time, c is the speed of light, f r is the range sampling frequency of the radar system; t r is the imaging time in the range direction; Step 2-4) Establish the conditional equation of the geometric positioning model for each ground target point, use the object space coordinates of the center point provided by the SAR image product as the initial value, and use the Newton iteration method based on the least squares principle to solve the object space coordinates of the ground target point. During the iteration, the partial derivatives of the object space coordinates of the ground target point are obtained according to the conditional equations of the geometric positioning model to obtain the partial derivative matrix, and update it. When the correction value of the object space coordinates relative to the previous iteration is less than the pre-configured threshold, the iteration is terminated to obtain the object space coordinates of the ground target point, and the mutual mapping of the image space coordinates and the object space coordinates is completed; The step 3) comprises the following steps: Step 3-1) Preprocess the SAR image: collect the first-level products of multi-view SAR images in the study area, perform radiation correction, and convert the backscatter information into intensity information; Step 3-2) Select a wind turbine target point in a flat area with simple surrounding features in the SAR image, and record its image coordinates; Step 3-3) using the initial geometric positioning model constructed in step 2), mapping the image space coordinates recorded in step 3-2) to object space coordinates, narrowing the search space, quickly locating the same-name points in the DOM, and recording the object space coordinates of the same-name points; Step 3-4) Obtain elevation information in DEM based on the object coordinates of the same-name points and generate ground control points for the same-name point pairs; Step 3-4) Repeat steps 3-2) to 3-4) to complete the collection of same-name point pairs, and divide the same-name point pairs into control point pairs and check point pairs, and establish a control point pair database and a check point pair database; The step 6) is specifically as follows: Step 6-1) According to the constructed initial geometric positioning model, pre-configured errors are added to the imaging parameters respectively, wherein the imaging parameters include: the orbital coordinates (X s ,Y s ,Z s ), starting slant distance R0 in range, slant distance resolution P in range r , system distance sampling frequency f r and the pulse repetition frequency f a ; Step 6-2) According to the working principle of the SAR imaging system, the sources of the SAR image geometric positioning error include the electronic delay inside the system, the time synchronization error between the SAR payload and the GNSS payload, the radar sampling frequency and the pulse repetition frequency error; Step 6-3) Determine the SAR geometric calibration model based on the geometric positioning error sources: Where Δf r and Δf a are the correction values ​​of the range sampling frequency and pulse repetition frequency respectively; t r ,t a is the imaging time in range and azimuth, t r0 ,t a0 is the starting time in range and azimuth, t delay is the atmospheric delay value, Δt r is the electronic delay inside the system, Δt a is the synchronization error between the SAR payload and the GNSS payload; Step 6-4) Perform SAR image-to-ground geometric positioning error analysis based on the added pre-configured error, and ignore Δf according to the analysis results. r and Δf a The impact on geometric positioning results in a simplified SAR geometric calibration model: Step 6-5) Determine the conditional equation of the geometric calibration model based on the SAR imaging geometric characteristics: The step 7) comprises the following steps: Step 7-1) Since the nonlinear equations are complex in form and usually difficult to directly obtain analytical solutions, the Taylor series is used to expand the calibration parameters in the conditional equations of step 6-5) to first-order terms to obtain a linear equation system. The Jacobian matrix of the linear equation system is shown below: Where X is the error correction parameter; Step 7-2) obtaining a preconfigured number of control point pairs in a control point pair database; Step 7-3) establishing a conditional equation of a geometric calibration model for each control point pair; Step 7-4) Based on the least squares principle, the Newton iteration method is used to update the Jacobian matrix based on the conditional equation of the geometric calibration model. When X is less than a pre-configured threshold, the iteration process ends and the most likely estimation of the error correction parameter is completed; Step 7-5) Update the initial geometric positioning model of the SAR image based on the error correction parameters to obtain an updated geometric positioning model.

2. The method for geometric calibration of satellite-borne SAR based on control points provided by wind turbines according to claim 1, characterized in that: The step 4) is specifically as follows: The positioning results of the check point pairs are calculated based on the initial geometric positioning model of each image and compared with the object coordinates of the check point pairs. The difference between the two is the geometric positioning error of the initial geometric positioning model.

3. The method for geometric calibration of satellite-borne SAR based on control points provided by wind turbines according to claim 1, characterized in that: The step 5) comprises the following steps: Step 5-1) Estimate the atmospheric parameters at any location and time based on the standard atmospheric model: Where Press is atmospheric pressure, Temp is Kelvin temperature, h is altitude, w is relative humidity, W press is the water vapor partial pressure; Step 5-2) Based on the atmospheric parameters calculated in step 5-1), the tropospheric atmospheric path delay is calculated using the Sastamonen model as follows: Where θ is the SAR satellite incident angle, φ is the satellite elevation angle, and D h is the dry delay, D w is the wet delay, ZTD is the tropospheric atmospheric delay; Step 5-3) Estimate the ionospheric path delay based on the Krobucher model: Where TEC is the ionospheric electron concentration; f is the radar center frequency; ZID is the ionospheric delay; Step 5-4) Determine the atmospheric delay value based on the tropospheric atmospheric delay and the ionospheric delay: t delay = ZTD+ZID.

4. The method for geometric calibration of satellite-borne SAR based on control points provided by wind turbines according to claim 1, characterized in that: The step 8) is specifically as follows: The initial geometric positioning model and the updated geometric positioning model are respectively evaluated for the positioning deviation of each check point pair using the check point pairs. The statistical indicators of all the check point pairs contained in each scene image are counted to compare the performance difference between the two.

5. The method for geometric calibration of satellite-borne SAR based on control points provided by wind turbines according to claim 4, characterized in that: The statistical indicators include minimum deviation, maximum deviation, mean value, and root mean square error.

6. The method for geometric calibration of satellite-borne SAR based on control points provided by wind turbines according to claim 1, characterized in that: In step 3-4), the point pairs with the same name are divided into control point pairs and check point pairs according to the spatial distribution of the control points, wherein the number of control point pairs in each scene image is not less than 4 and the number of check point pairs is not less than 4.

Images