Stereoscopic SAR control point automatic extraction method based on variance component estimation

Abstract

The invention provides a stereo SAR control point automatic extraction method based on variance component estimation. The method comprises the following steps: acquiring double-track or multi-track multi-scene time sequence SAR images and DEM data of the same target area, registering the multi-scene time sequence SAR images of the same track, extracting permanent scatterers on the images for geocoding, acquiring rough geographic coordinates of the permanent scatterers, and inversely calculating initial homonymy pixel coordinates of the permanent scatterers on the other track, so as to obtain a target target area; a plurality of candidate homonymous pixels and sub-pixel-level coordinates of each candidate homonymous pixel are obtained, distance direction and azimuth direction time observation values of the candidate homonymous pixels are corrected, an adjustment function model about the time observation values and three-dimensional geographic coordinates is established, a linearized function model is iteratively solved through variance component estimation, and the time observation values of the candidate homonymous pixels and the three-dimensional geographic coordinates of the candidate homonymous pixels are calculated. Obtaining a variance-covariance matrix of the estimation parameters; and determining the final homonymous pixels and the three-dimensional geographic coordinates thereof based on the standard deviation threshold. According to the invention, the positioning calculation precision of the stereo SAR control points is improved, and a large number of high-precision ground control points can be extracted.

Description

Technical Field

[0001] The present invention belongs to the technical field of image processing, and particularly relates to an automatic extraction method for stereo SAR control points based on variance component estimation. Background Art

[0002] Ground Control Points (GCPs) are key elements in photogrammetry. They have known spatial coordinate information. By marking the positions of GCPs in images, the association between the images and real-world coordinates can be established, thus enabling high-precision tracking, positioning, and analysis of target objects. Such accurate spatial coordinate information is of great significance for applications such as surveying, remote sensing, and 3D modeling in aerial photogrammetry. In addition, GCPs can also be used to correct images, improving the geometric accuracy and map projection quality of images, especially for Synthetic Aperture Radar (SAR) imaging technology. However, due to the influence of various factors such as satellite clock errors, orbital errors, atmospheric delays, and surface changes on the signal transmission of the SAR satellite system, the geometric positioning accuracy is reduced. Therefore, accurately positioning the monitoring points of high-resolution SAR images is a difficulty and pain point at the present stage. Moreover, temporal InSAR (Interferometric Synthetic Aperture Radar) is a relative technology, and the height estimation error of monitoring points and the reference point error are also the main error sources affecting geocoding. The relatively low positioning accuracy is difficult to meet the fine monitoring requirements of InSAR, which brings difficulties to the interpretation of deformation results.

[0003] Stereo SAR (Stereo Synthetic Aperture Radar) is a technology that uses dual / multi-angle SAR images with a certain intersection angle to invert the three-dimensional information of the scene and can be used to extract image ground control points. A prerequisite for stereo SAR positioning is that homologous points must be visible from multi-aspect SAR images, such as lamp posts, corner points, and intersection points. However, due to the complex layover structure and different imaging angles, this type of homologous points rarely occur and are difficult to identify with the naked eye. Existing stereo SAR methods mainly include two parts: automatic detection of homologous points and precise positioning of stereo SAR. The extraction of homologous points depends on external optical data and vectorized road grid data, which are difficult to detect automatically, and the positioning accuracy of stereo SAR calculation is limited. Summary of the Invention

[0004] In view of the above problems, the present invention proposes a method for automatic extraction of stereo SAR control points based on variance component estimation. Based on the traditional stereo SAR technology, the method adopts conditional adjustment of additional parameters and variance component estimation method to solve the three-dimensional geographic coordinates of the same-name points in the image, which provides a new idea for the detection of same-name points in dual / multi-angle SAR images, can automatically extract a large number of same-name pixels, and provides a basis for the automatic extraction of SAR control points, thereby solving the problems of difficulty in automatic detection of same-name points and low stereo SAR positioning accuracy when extracting ground control points using dual / multi-angle SAR images.

[0005] In order to achieve the above technical objectives, the present invention provides a method for automatically extracting control points of a stereo SAR based on variance component estimation, which is characterized by comprising the following steps:

[0006] S1. Obtain dual-track or multi-track multi-view time series SAR images and DEM data of the same target area, and perform registration on the multi-view time series SAR images of the same track;

[0007] S2. Extract permanent scatterers on multi-view time series SAR images of different orbits using amplitude deviation or amplitude threshold method;

[0008] S3. Geocoding the permanent scatterers on track 1 extracted in step S2 based on the DEM data to obtain the rough geographic coordinates of each permanent scatterer;

[0009] S4. Using the rough geographic coordinates of the permanent scatterer on track one obtained in step S3, the SAR pixel coordinates of the permanent scatterer on another track are back-calculated by the range-Doppler-ellipsoid equation, that is, the initial pixel coordinates of the same name;

[0010] S5. Expand the search range with the initial pixel coordinates of the same name of the permanent scatterer on track one on another track in step S4 as the center, and extract multiple candidate pixels of the same name of the permanent scatterer on track one on another track;

[0011] S6. Obtain the sub-pixel coordinates of each candidate pixel with the same name based on the point target analysis, and calculate the azimuth time observation value and the range time observation value of the corresponding sub-pixel coordinates in sequence;

[0012] S7. Correct the azimuth time observation value and the distance time observation value of each candidate pixel with the same name calculated in step S6 using external monitoring data;

[0013] S8. for each candidate pixel with the same name, two or more SAR images are used to combine multiple range-Doppler equations to establish an adjustment function model about time observation values ​​and three-dimensional geographic coordinates;

[0014] S9. Linearize the adjustment function model established in step S8, and iteratively solve the linearized adjustment function model through variance component estimation. Based on the iterative weighted least squares algorithm, obtain the three-dimensional coordinates of the ground target points corresponding to each candidate homologous pixel, and estimate the coordinate variance-covariance matrix;

[0015] S10. Based on the three-dimensional coordinates of the ground target points and the coordinate variance-covariance matrix of all candidate homologous pixels obtained in step S9, determine the final homologous pixels of the permanent scatterers on orbit 1 and their three-dimensional geographical coordinates based on the standard deviation threshold. And the ground points corresponding to the final homologous pixels are the ground control points of the stereo SAR.

[0016] A further technical solution of the present invention: In step S1, collect SAR images of double-track or multi-track time series covering the target area. These data are all data of the same satellite, and the number of collected images is not less than two scenes, and the two scenes of images come from different orbits; at the same time, obtain external DEM data within the corresponding area range, including SRTM DEM and TanDEM;

[0017] In step S1, the multi-scene time series SAR images of the same orbit will be registered. Specifically, for the multi-scene time series SAR images of the same orbit, by comparing the time baseline, spatial baseline, and Doppler centroid frequency difference between the SAR images, select the SAR image corresponding to the maximum correlation coefficient according to the comprehensive correlation function as the common master image, and use the remaining SAR images as slave images to be registered with the master image respectively. The registration method uses intensity cross-correlation registration; the intensity cross-correlation registration is mainly used to evaluate the similarity between two images by calculating the cross-correlation coefficient between them, so as to realize the registration of the images; specifically, by calculating the cross-correlation value at each position of the two images, find the position with the maximum cross-correlation value, so as to determine the best match of the two images. The higher the cross-correlation coefficient, the more similar the two images are; for two images I and J, their cross-correlation function R(x,y) is defined as:

[0018]

[0019] where I(i,j) and J(i+x,j+y) respectively represent the pixel values of the two images at the corresponding positions;

[0020] (i,j) represents the pixel position of the image, and (x,y) represents the displacement.

[0021] A preferred technical solution of the present invention: In step S2, when the number of images is greater than 15 scenes, the amplitude deviation index threshold method is used for permanent scatterer detection. The permanent scatterer corresponds to a high-reflection target that is stable in time. Specifically, the amplitude deviation threshold method is used to identify stable scatterers by calculating the ratio of the amplitude standard deviation to the average amplitude; the amplitude deviation DA It is expressed by the following formula:

[0022]

[0023] where σ φ represents the phase standard deviation;

[0024] m A and σ A respectively represent the amplitude time average of the pixel, i.e., the mathematical expectation, and the amplitude standard deviation;

[0025] When the number of images is less than or equal to 15, the amplitude threshold method is used for high-reflection scatterer detection. The amplitude threshold method is adopted to find the high-coherence point targets with strong reflection in the amplitude sequence, which are the permanent scatterers.

[0026] A further technical solution of the present invention: In step S3, based on the external DEM data, through the range-Doppler-ellipsoid equation, the SAR data is converted from the radar image coordinate system to the geographic coordinate system to obtain the rough three-dimensional geographic coordinates of the permanent scatterers; the range-Doppler-ellipsoid equation in steps S3 and S4 is specifically as follows:

[0027]

[0028] where ct r represents the slant range R from the ground target point P to the satellite; c is the speed of light;

[0029] t r is the one-way range-direction propagation time; λ is the radar signal wavelength;

[0030] m and n are respectively the semi-major axis and semi-minor axis of the reference ellipsoid; h is the geodetic height of the ground target P;

[0031] S S is the position vector of the satellite at time t a ; V S is the velocity vector of the satellite at time t a ;

[0032] S P is the position of the ground target point P at time t a ;

[0033] f D is the Doppler frequency of the target point P;

[0034] At time t a , the position vector of the satellite is S S =(X S , Y S , Z S ) T , and the velocity vector is The position of the ground target point P is S P =(X P , Y P , Z P ) T ; f D is the Doppler frequency of the target point P. After the SAR raw data is focused, the R-D equation refers to the zero Doppler frequency, that is, f D =0. At this time, the position vector and velocity vector of the satellite are also the values corresponding to the zero Doppler frequency.

[0035] A further technical solution of the present invention: In step S5, with the initial homologous pixel as the center, the search range is expanded with a 10×10 rectangular frame as the radius, and multiple candidate homologous pixels centered on the initial homologous pixel are extracted;

[0036] In step S6, the azimuth time observation value t a and the range time observation value t r of each candidate homologous pixel are indirectly calculated through the pixel coordinates (x pt , y pt ), and their calculation formulas are as follows:

[0037]

[0038] where t a0 represents the Coordinated Universal Time of the starting pixel in the azimuth direction;

[0039] t r0 is the ratio time of the round-trip distance of the starting pixel in the range direction to the speed of light;

[0040] Δt r =1 / rsr, Δt a =1 / prf, rsr is the range sampling rate, and prf is the pulse repetition frequency;

[0041] All the above parameters can be found in the image parameter file.

[0042] A further technical solution of the present invention: The point target analysis in step S6 includes the following two methods:

[0043] One is to apply spectral zero-padding to the image subset containing the target pixel, and use a SINC interpolator to perform oversampling on the complex data. The oversampling factor is equal to the number of inserted zeros;

[0044] Another method for point target analysis is to identify sub-pixel coordinates by calculating the centroid coordinates of the oversampled image. The specific process is as follows: First, identify the corner reflector pixels in the SAR intensity image by eye, then intercept the n×n window sub-block containing the corner reflector pixels, and perform m-fold oversampling on the image sub-block using the SINC function or bilinear interpolation function. Finally, calculate the centroid coordinates of each row and each column of the oversampled image and the average intensity value

[0045]

[0046] where I and f represent the pixel intensity values before and after oversampling, respectively;

[0047] (i, j, k) are pixel coordinates; i′, j′ = 1, 2, …, n;

[0048] Finally, calculate the intensity peak center coordinates of the image subset, that is, the sub-pixel positions (x c , y c ) of the main scatterer in the range and azimuth directions:

[0049]

[0050] A further technical solution of the present invention: In step S7, correcting the azimuth time observation value and range time observation value of the candidate homologous pixels specifically means using external monitoring data to model the solid tide and atmospheric delay errors, so as to eliminate the influence of the solid tide and atmospheric delay errors on the SAR positioning accuracy;

[0051] The local area movement of the surface displacement caused by the solid tide in the east, north, and vertical directions is obtained using open-source software, and the three-dimensional deformation is projected onto the azimuth and range directions:

[0052]

[0053] where ξ az , ξ rg represent the azimuth and range delays respectively; β represents the azimuth angle of the satellite flight;

[0054] θ is the local incident angle at imaging; ξ E , ξ N , ξ U represent the SET displacements in the east, north, and vertical directions respectively;

[0055] The azimuth and range correction amounts calculated based on the above formulas can weaken the timing error caused by SET;

[0056] Atmospheric delay mainly includes tropospheric delay and ionospheric delay. By using classical geodetic methods and the nearby GNSS to estimate the tropospheric zenith delay above the InSAR monitoring point, the tropospheric delay at the local point can be obtained to reduce the influence of tropospheric delay:

[0057]

[0058] Among them, I ZPD is the tropospheric zenith delay; θ is the local incident angle during imaging;

[0059] h SAR represents the ground elevation of the SAR monitoring point;

[0060] h GNSS represents the ground elevation of the GNSS monitoring point near the SAR monitoring point; h 0 = 6000 m;

[0061] The ionosphere is distributed between 100 and 1500 km above sea level and is an ionized atmospheric region in the atmosphere. The influence generated by the ionosphere is generally at the decimeter level. The ionospheric vertical total electron content map I provided by the IGS analysis center is used VETC to compensate for the ionospheric delay in the range direction:

[0062]

[0063] Among them, f represents the radar carrier frequency, K = 40.28 m 3 s -2 ;

[0064] The spatial resolutions of the VTEC grid map provided by the IGS analysis center in longitude and latitude are 2.5° and 5° respectively, and the time resolution is 1 hour. The I at the SAR monitoring point VETC is calculated by the plane interpolation method, and the VTEC at the target point position is interpolated from four adjacent grid points.

[0065] A further technical solution of the present invention: The adjustment function model established for the time observation value and the three-dimensional geographical coordinates in the step S8 includes:

[0066] According to the geometric relationship between the SAR satellite and the ground point target, high-precision stereo SAR positioning is realized by the oblique intersection of the multi-orbit range-Doppler equations; for one scene of SAR image, the time observation value and the three-dimensional coordinates of the unknown target point are implicitly included in the following two trajectory models:

[0067]

[0068] Here, for the target control point P, assuming its three-dimensional geographical coordinates are SP =(X P ,Y P ,Z P ) T , at azimuth time t a moment, the position vector of the satellite is S S =(X S ,Y S ,Z S ) T , and the velocity vector is c is the speed of light, t r is the one-way range propagation time; S S and V S respectively represent the position and velocity vectors of the SAR sensor at the target imaging time t a , which are calculated by using the sixth-order polynomial model of the discrete satellite orbit state vector provided in the SAR parameter file. The specific calculation formula is as follows:

[0069]

[0070] where: n represents the n-th polynomial, and the coefficients of the polynomial model are (a i ,b i ,c i ); these coefficients are estimated by the least squares method;

[0071] In order to calculate the three-dimensional coordinates of the target point, image acquisitions of more than two different orbits are required.

[0072] Further technical solution of the present invention: In step S9, by linearizing formula (14) in step S8, the observation equations of multiple SAR images can be expressed as conditional adjustment of additional unknown parameters:

[0073] Bv + Ax + w = 0 (16)

[0074] where, B is the conditional matrix; v is the residual of the observed value; A represents the coefficient matrix of the unknown parameter;

[0075] x represents the coordinate correction of the target point, and w is the closing error;

[0076] Since the coefficient matrix B is invertible, the above formula can be written as:

[0077] v = B'x - l' (17)

[0078] where, B' = -B -1 A, l' = B -1 w;

[0079] Equation (16) is converted into the classical least squares formula, and the coordinate correction of the target point Specifically, it is calculated according to the following formula:

[0080]

[0081] In the above formula, N represents the normal equation matrix of the coefficient matrix; T is the transpose process; P is the weight;

[0082] By iteratively adding the coordinate correction amount to the initial coordinates The three-dimensional geographic coordinates of the ground target points corresponding to each candidate homologous pixel are obtained; the variance-covariance matrix of its coordinate estimation parameters is as follows:

[0083]

[0084] Among them, represents the variance-covariance matrix of the coordinate estimation parameters;

[0085] σ x ,σ y ,σ z are the coordinate standard deviations; σ xy ,σ xz ,σ yz represents the coordinate covariance;

[0086] n is the number of observations; m represents the number of unknown parameters.

[0087] A further technical solution of the present invention: In the step S10, the internal consistency accuracy of the coordinates after stereo SAR positioning is used as the judgment criterion for homologous pixels, and the final homologous pixels and their three-dimensional geographic coordinates are determined by setting the threshold values of the standard deviation and ∑|σ i |, where the standard deviation of the homologous pixels is the smallest.

[0088] Stereo SAR can perform absolute three-dimensional positioning on natural permanent scatterers, which enables the generation of ground control points using only SAR data. The prerequisite for obtaining high-precision GCPs by this method is the correct detection of homologous scatterers in SAR images with different observation geometries. In the present invention, a processing strategy for automatically detecting the same targets in SAR images of the same area taken from different orbits is described. In addition, a complete automatic generation process of GCPs based on stereo SAR is provided, which can extract a large number of high-precision ground control points and has high engineering application value.

[0089] Based on the traditional stereo SAR technology, the present invention adopts the conditional adjustment with additional parameters and the variance component estimation method to calculate the three-dimensional geographic coordinates of the corresponding points in the images, providing a new idea for the detection of corresponding points in dual / multi-angle SAR images, which can automatically extract a large number of corresponding pixels and lay a foundation for the automatic extraction of SAR control points. Compared with the existing stereo SAR control point generation technology, the automatic extraction method of stereo SAR control points in the present invention can automatically detect the corresponding points in SAR images without external optical images and vectorized road grid data. In addition, through point target analysis, time observation value correction, function adjustment and variance component estimation, the present invention improves the positioning and calculation accuracy of stereo SAR control points. Generally speaking, the present invention provides a complete set of automatic generation processes for stereo SAR control points, which can extract a large number of high-precision ground control points and has high engineering application value. Description of the Drawings

[0090] Figure 1 is the flowchart of the present invention;

[0091] Figure 2 is the search schematic diagram of candidate corresponding pixels in the embodiment of the present invention;

[0092] Figure 3 is the point target analysis schematic diagram in the embodiment of the present invention;

[0093] Figure 4-a is the three-dimensional surface displacement caused by solid tides at different times. Projecting it onto the azimuth and range directions can obtain the propagation error schematic diagram;

[0094] Figure 4-b is the range direction error schematic diagram caused by tropospheric delay;

[0095] Figure 4-c is the range direction error schematic diagram caused by ionospheric delay;

[0096] Figure 5-a is the final corresponding pixel schematic diagram on orbit 1 in the embodiment;

[0097] Figure 5-b is the final corresponding pixel schematic diagram on another corresponding orbit in the embodiment. Detailed Embodiment

[0098] It should be noted that, without conflict, the embodiments in this example and the features in the embodiments can be combined with each other. Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0099] The embodiment provides a method for automatically extracting stereo SAR control points based on variance component estimation, as Figure 1 shown, which includes the following steps:

[0100] S1. According to actual needs, obtain multi-scene time series SAR images, DEM data and their parameter files of the same target area on double orbits or multiple orbits. The number of collected SAR images is not less than two scenes, and the two scenes of images come from different orbits. Generally, the greater the imaging angle difference, the more robust the spatial intersection and the higher the positioning accuracy. In addition, collect external DEM data within the corresponding area range, such as SRTM DEM and TanDEM, and the data resolution depends on cost and requirements; for multi-scene time series SAR images on the same orbit, in order to enable effective comparison, fusion and analysis of multi-source images, it is necessary to register the same-orbit images to the sub-pixel level.

[0101] Intensity cross-correlation registration is a commonly used image registration method, which is mainly used to evaluate the similarity between two images by calculating the cross-correlation coefficient between them, so as to achieve image registration. The basic principle of intensity cross-correlation registration is to calculate the cross-correlation value at each position of the two images, find the position with the largest cross-correlation value, and thus determine the best match of the two images. The cross-correlation coefficient reflects the similarity degree of the two images in space, and the higher its value, the more similar the two images are. For two images I and J, their cross-correlation function R(x, y) can be defined as:

[0102]

[0103] where I(i, j) and J(i + x, j + y) respectively represent the pixel values of the two images at the corresponding positions;

[0104] (i, j) represents the pixel position of the image, and (x, y) represents the displacement amount.

[0105] S2. Use the amplitude deviation or amplitude threshold method to extract the permanent scatterers on the multi-scene time series SAR images of different orbits;

[0106] When the number of images is large, generally when the number is greater than 15 scenes, the amplitude deviation index threshold method can be used for permanent scatterer detection. Permanent scatterers usually correspond to high-reflection targets that are stable over time, such as lamp posts, corners, etc. Specifically, the amplitude deviation threshold method is used to identify stable scatterers by calculating the ratio of the amplitude standard deviation to the average amplitude; the amplitude deviation D A is expressed by the following formula:

[0107]

[0108] where σ φ represents the phase standard deviation;

[0109] m A and σ A respectively represent the average amplitude over time (mathematical expectation) and the amplitude standard deviation of the pixel.

[0110] When the number of images is small, generally less than or equal to 15 scenes, the amplitude threshold method can be used for high-reflection scatterer detection. Calculate the amplitude value of each SAR image to obtain the amplitude mean. Set a threshold based on the amplitude mean. Usually, the minimum value of the amplitude mean is selected as the threshold, and the pixel points with amplitude values higher than this threshold are identified as permanent scatterers.

[0111] S3. Geocode the permanent scatterers on orbit 1 extracted in step S2 based on DEM data to obtain the approximate geographical coordinates of each permanent scatterer; Geocoding is to determine the three-dimensional geographical coordinates of SAR image pixels. According to the time sequence and speed of the radar flight in the azimuth direction during SAR imaging, and the echo delay and the speed of light received in the range direction, the two-dimensional geometric coordinates (a, r) of the radar image can be determined. Since the intersection of the radar beam center and the Earth's surface determines the actual geographical location of any pixel in the image, the ground coordinates (x, y, z) corresponding to the radar image coordinates can be calculated using a strict mathematical model. The most commonly used model is the range-Doppler-ellipsoid equation (abbreviated as the R-D equation):

[0112]

[0113] where ct r represents the slant range R from the ground target point P to the satellite; c is the speed of light;

[0114] t r is the one-way range-direction propagation time; λ is the radar signal wavelength;

[0115] m and n are respectively the semi-major axis and semi-minor axis of the reference ellipsoid; S S is the position vector of the satellite at time t a ;

[0116] VS The velocity vector of the satellite at time t; S a The position of the ground target point P at time t; P For t a The position of the ground target point P at time t;

[0117] f D The Doppler frequency of the target point P; h is the geodetic height of the ground target P.

[0118] At time t a The position vector of the satellite is S S =(X S , Y S , Z S ) T , and the velocity vector is The position of the ground target point P is S P =(X P , Y P , Z P ) T . f D The Doppler frequency of the target point P. After focusing the SAR raw data, the R-D equation refers to the zero Doppler frequency, i.e., f D =0. At this time, the position vector and velocity vector of the satellite are also the values corresponding to the zero Doppler frequency. It can be seen that the azimuth and range timing errors and the external DEM error are three factors affecting the positioning accuracy of the R-D equation. Especially for high-resolution SAR images, the low-resolution DEM reduces the accuracy of geocoding. Since the temporal InSAR is a relative measurement technique, the reference point error is also the main error source affecting the positioning accuracy. The stereo SAR positioning only considers the range-Doppler equations of different orbits, effectively avoiding the influence of the external DEM error and providing more accurate ground coordinates.

[0119] S4. Use the rough geographical coordinates of the permanent scatterers on orbit 1 obtained in step S3 to inversely calculate the SAR pixel coordinates of the permanent scatterers on another orbit through the range-Doppler-ellipsoid equation (R-D equation). These pixels are regarded as the initial homologous pixels between different orbits.

[0120] S5. Expand the search range centered on the initial homologous pixel coordinates of the permanent scatterers on orbit 1 in step S4 on another orbit. Expand the search range with a 10×10 rectangular box as the radius, and extract multiple candidate homologous pixels of the permanent scatterers on orbit 1 within the rectangular box on another orbit; specifically refer to Figure 2 , where Figure 2 (a) represents the permanent scatterers on orbit 1, Figure 2 (b) represents multiple candidate homologous pixels on another orbit, Figure 2(c) represents the standard deviation value corresponding to each candidate homologous pixel;

[0121] S6. Obtain the sub-pixel coordinates of each candidate homologous pixel based on point target analysis, and sequentially calculate the azimuth time observation value and the range time observation value corresponding to the sub-pixel coordinates;

[0122] Refer to Figure 3 , for each pair of candidate homologous pixels, extract the sub-pixel coordinates by point target analysis to improve the stereo SAR positioning accuracy; where Figure 3 Left represents the original pixel coordinates of the permanent scatterer and its corresponding backscattering intensity value, Figure 3 Right represents the sub-pixel coordinates of the permanent scatterer and its corresponding backscattering intensity value.

[0123] For the existing SAR satellite imaging mode, the preprocessed level 1 single-look complex image is converted into a two-dimensional image grid with range and azimuth radar coordinates. Each grid represents a pixel, and the grid coordinates correspond to the pixel coordinates. The azimuth time observation value t a and the range time observation value t r of each candidate homologous pixel are indirectly calculated through the pixel coordinates (x pt , y pt ), and their calculation formulas are as follows:

[0124]

[0125] Among them, t a0 represents the Coordinated Universal Time (UTC) time of the starting pixel in the azimuth direction;

[0126] t r0 is the ratio time of the round-trip distance of the starting pixel in the range direction and the speed of light;

[0127] Δt r = 1 / rsr, Δt a = 1 / prf, rsr is the range sampling rate, and prf is the pulse repetition frequency;

[0128] These parameters can all be found in the image parameter file.

[0129] Assume that the given image resolution is 1 m and the target is an SAR measurement accuracy at the centimeter level. Then, point target analysis can provide a peak position with approximately 1 / 100 pixel. An effective way to achieve this accuracy is to apply spectral zero-padding to the image subset containing the target pixel, which is equivalent to oversampling the complex data using a SINC interpolator. The oversampling factor is equal to the number of inserted zeros. However, a large amount of padding will involve the computationally time-consuming inverse Fourier transform in spectral calculation. Another point target analysis method is to identify the sub-pixel coordinates by calculating the centroid coordinates of the oversampled image. First, identify the corner reflector pixels by eye in the SAR intensity image. Then, intercept the sub-block of the n×n window containing the corner reflector pixels, and perform m-fold oversampling on the image sub-block using the SINC function or the bilinear interpolation function. Finally, calculate the centroid coordinates of each row and each column of the oversampled image and the average intensity value

[0130]

[0131] where I and f represent the pixel intensity values before and after oversampling, respectively;

[0132] (i,j,k) are the pixel coordinates; i′,j′ = 1,2,…,n;

[0133] Finally, calculate the intensity peak center coordinates of the image subset, that is, the sub-pixel positions (x c ,y c ) of the main scatterer in the range and azimuth directions:

[0134]

[0135] S7. Use external monitoring data to correct the azimuth time observation value and the range time observation value of each candidate corresponding pixel calculated in step S6;

[0136] In order to correct the errors of the azimuth time observation value and the range time observation value, the present invention models the solid tide and atmospheric delay errors using external monitoring data, thereby eliminating the influence of these factors and improving the stereo SAR positioning accuracy. Refer to Figures 4-a to 4-c , Figure 4-a represents the three-dimensional surface displacement caused by the solid tide at different times. The propagation error can be obtained by projecting it onto the azimuth and range directions, Figure 4-b represents the range error caused by the tropospheric delay, Figure 4-c represents the range error caused by the ionospheric delay.

[0137] The currently achievable SAR observation accuracy mainly requires correction of the Solid Earth Tides (SET). SET is the deformation of the solid Earth's surface caused by the gravitational forces of the sun and the moon, and the periodic deformation caused within a day can reach several decimeters. The local area movement of the ground displacement caused by SET in the east, north, and vertical directions can be obtained using open-source software, and the three-dimensional deformation is projected onto the azimuth and range directions:

[0138]

[0139] Among them, ξ az , ξ rg represent the azimuth and range delays respectively;

[0140] β represents the azimuth angle of the satellite flight; θ is the local incident angle during imaging;

[0141] ξ E , ξ N , ξ U represent the SET displacements in the east, north, and vertical directions respectively.

[0142] The azimuth and range correction amounts calculated based on the above formulas can weaken the timing error caused by SET.

[0143] The atmospheric delay mainly includes tropospheric delay and ionospheric delay. The tropospheric delay is the layer of the atmosphere closest to the Earth's surface, which is related to the terrain height and is located in the lowest layer of the atmosphere. The troposphere has the highest density in the atmosphere and has the most serious impact on the positioning accuracy, and the resulting error can even reach several meters. To reduce the influence of tropospheric delay, classical geodetic methods can be used to estimate the tropospheric zenith delay (ZPD) above the InSAR monitoring point using adjacent GNSS, so as to obtain the tropospheric delay at the local point:

[0144]

[0145] Among them, I ZPD is the tropospheric zenith delay; θ is the local incident angle during imaging;

[0146] h SAR represents the ground elevation of the SAR monitoring point;

[0147] h GNSS represents the ground elevation of the GNSS monitoring point near the SAR monitoring point; h 0 = 6000 meters.

[0148] This method assumes that the troposphere within a certain range has spatial correlation. Therefore, the premise for calculating the tropospheric delay using this method is that the target point needs to be within 50 kilometers of the GNSS observation station. The farther away from the observation station, the greater the estimated tropospheric delay deviation. After years of research, more advanced tropospheric delay estimation methods have now emerged. When there is no GNSS station near the study area, the ZPD map provided by the Generic Atmospheric Correction Online Service for InSAR (GACOS) can be considered to calculate the tropospheric delay. The advantages of GACOS tropospheric maps are global availability, near real-time and high spatial resolution.

[0149] The ionosphere is located between 100 and 1500 km above sea level. It is a part of the atmosphere that is ionized. The impact of the ionosphere is generally at the decimeter level. The ionosphere vertical total electron content (VTEC) provided by the IGS Analysis Center (University of Bern, Switzerland) is used here. VETC To compensate for the ionospheric delay in the range upwards:

[0150]

[0151] Where, f represents the radar carrier frequency, K = 40.28m 3 s -2 ;

[0152] The VTEC grid map provided by the IGS Analysis Center has a spatial resolution of 2.5° and 5° in longitude and latitude, respectively, and a temporal resolution of 1 hour. VETC It is calculated by plane interpolation method, and the VTEC at the target point is interpolated through four adjacent grid points.

[0153] S8. For each candidate pixel with the same name, two or more SAR images are used to combine multiple range-Doppler equations to establish an adjustment function model for time observations and three-dimensional geographic coordinates; based on the geometric relationship between the SAR satellite and the ground point target, high-precision stereo SAR positioning can be achieved through the oblique intersection of multi-track range-Doppler equations. For one of the SAR images, the time observations and the three-dimensional coordinates of the unknown target point are implicitly included in the two trajectory models:

[0154]

[0155] Here, for the target control point P, assume that its three-dimensional geographic coordinates are S P =(X P,Y P ,Z P ) T , at time t a (azimuth time), the position vector of the satellite is S S = (X S ,Y S ,Z S ) T , and the velocity vector is where c is the speed of light and t r is the one-way range propagation time. S S and V S respectively represent the position and velocity vectors of the SAR sensor at the target imaging time t a , and can be calculated by using a sixth-order polynomial model of the discrete satellite orbit state vector provided in the SAR parameter file:

[0156]

[0157] where n represents the n-th polynomial, and the coefficients of the polynomial model are (a i ,b i ,c i ), and these coefficients can be estimated by the least squares method.

[0158] To calculate the three-dimensional coordinates of the target point, more than two image acquisitions of different orbits are required. Two orbits provide 4 equations, and 4 equations can solve for three unknown coordinate parameters. The more images there are, the more equations there are, and the more robust the parameter solution is.

[0159] S9. Linearize the adjustment function model established in step S8, and iteratively solve the linearized adjustment function model through variance component estimation. Based on the iterative weighted least squares algorithm, obtain the three-dimensional coordinates of the ground target point corresponding to each candidate homologous pixel, and estimate the coordinate variance-covariance matrix; specifically as follows: Through the linearized equation (14), the observation equations of multiple SAR images can be expressed as conditional adjustment with additional unknown parameters:

[0160] Bv + Ax + w = 0 (16)

[0161] where B is the conditional matrix; v is the observation residual; A represents the coefficient matrix of the unknown parameters;

[0162] x represents the coordinate correction of the target point, and w is the closure error;

[0163] Since the coefficient matrix B is invertible, the above equation can be written as:

[0164] v = B'x - l' (17)

[0165] where B ′=-B- 1 A, l ′ =B- 1 w;

[0166] Equation (16) is converted into the classical least - squares formula, and the coordinate correction of the target point can be calculated according to the following formula:

[0167]

[0168] Here, N represents the normal equation matrix of the coefficient matrix, T represents the transpose operation, and P represents the weight.

[0169] By iteratively adding the coordinate correction to the initial coordinates the three - dimensional geographical coordinates of the target point are obtained.

[0170] A special advantage of least - squares estimation is that it provides complete error statistics for the estimation results. One of the quality indicators used in the present invention is the variance - covariance of the coordinate estimation parameters. By calculating the square root of the diagonal elements, the standard deviation of the coordinates in each direction can be derived:

[0171]

[0172] Among them, represents the variance - covariance matrix of the coordinate estimation parameters, which is used to describe the internal precision of the target point coordinate estimation;

[0173] σ x , σ y , σ z is the coordinate standard deviation; σ xy , σ xz , σ yz represents the coordinate covariance;

[0174] n is the number of observations; m represents the number of unknown parameters.

[0175] Since the prior mean square errors of the two types of observations are not precisely known, the weight P is iteratively determined using the variance component estimation method. First, reasonable prior mean square errors of the range - time and azimuth - time observations are given for the first weight determination. For example, assume the initial weight of the range - time observation The initial weight of the azimuth - time observation After the first solution, the unit weight variance of the two types of observations is estimated using the residuals

[0176]

[0177] Among them, v ρ and v aDenote the corrections of the range - direction and azimuth - direction observations respectively; p ρ and p a are their weights;

[0178] n ρ and n a are the numbers of observations in two directions; trace represents the trace of a matrix.

[0179] Re - weight based on the estimated variance of unit weight:

[0180]

[0181] where k represents the number of iterations.

[0182] Iteratively perform the adjustment and re - weighting steps until the difference between the two is less than a predetermined threshold, stop the iterative operation, and output the three - dimensional geographical coordinates corresponding to the estimated candidate homologous pixels.

[0183] S10. Based on the three - dimensional coordinates of the ground target points and the coordinate variance - covariance matrix of all candidate homologous pixels obtained in step S9, determine the final homologous pixels and their three - dimensional geographical coordinates of the permanent scatterers on orbit 1 based on the standard deviation threshold, and the ground points corresponding to the final homologous pixels are the ground control points of the stereo SAR. For specific reference, see Figure 5-a and Figure 5-b where is the final homologous pixel on orbit 1 in 5 - a, Figure 5-b is the final homologous pixel on another orbit.

[0184] The standard deviation threshold represents the standard deviation sum and threshold of the coordinates after stereo SAR positioning for a certain target pixel and candidate homologous pixel. Experiments prove that the standard deviation sum of the true homologous pixel is the smallest. Therefore, the present invention uses the internal consistency accuracy of the coordinates after stereo SAR positioning as the judgment criterion for homologous pixels:

[0185] ∑|σ i | = min(24)

[0186] where i = x, y, z. By setting the standard deviation sum ∑|σ i | threshold (such as set to 3), the final homologous pixels and their three - dimensional geographical coordinates can be determined, and the target points corresponding to the homologous pixels are regarded as ideal ground control points.

[0187] The above is only an embodiment of the present invention, and its description is relatively specific and detailed, but it should not be construed as a limitation on the scope of the present invention. It should be noted that for those of ordinary skill in the art, without departing from the concept of the present invention, several deformations and improvements can be made, and these all belong to the protection scope of the present invention. Therefore, the protection scope of the present invention should be subject to the appended claims.

Claims

1. A method for automatic extraction of stereo SAR control points based on variance component estimation, characterized in that: The following steps are involved: S1. Obtain dual-track or multi-track multi-view time series SAR images and DEM data of the same target area, and perform registration on the multi-view time series SAR images of the same track; S2. Extract permanent scatterers on multi-view time series SAR images of different orbits using amplitude deviation or amplitude threshold method; S3. Geocoding the permanent scatterers on track 1 extracted in step S2 based on the DEM data to obtain the rough geographic coordinates of each permanent scatterer; S4. Using the rough geographic coordinates of the permanent scatterer on track one obtained in step S3, the SAR pixel coordinates of the permanent scatterer on another track are back-calculated by the range-Doppler-ellipsoid equation, that is, the initial pixel coordinates of the same name; S5. Expand the search range with the initial pixel coordinates of the same name of the permanent scatterer on track one on another track in step S4 as the center, and extract multiple candidate pixels of the same name of the permanent scatterer on track one on another track; S6. Obtain the sub-pixel coordinates of each candidate pixel with the same name based on the point target analysis, and calculate the azimuth time observation value and the range time observation value of the corresponding sub-pixel coordinates in sequence; S7. Correct the azimuth time observation value and the distance time observation value of each candidate pixel with the same name calculated in step S6 using external monitoring data; S8. for each candidate pixel with the same name, two or more SAR images are used to combine multiple range-Doppler equations to establish an adjustment function model about time observation values ​​and three-dimensional geographic coordinates; The adjustment function model of setting up in S9. linearization step S8, and by variance component estimation iterative solution linearized adjustment function model, obtain the corresponding ground target point three-dimensional coordinates of each candidate pixel of the same name based on iterative weighted least squares algorithm, and estimate coordinate variance-covariance matrix; S10. According to the three-dimensional coordinates of the ground target points and the coordinate variance-covariance matrix of all candidate pixels with the same name obtained in step S9, the final pixels with the same name and their three-dimensional geographic coordinates of the permanent scatterer on track one are determined based on the standard deviation threshold, and the ground point corresponding to the final pixel with the same name is the ground control point of the stereo SAR.

2. The method for automatically extracting 3D SAR control points based on variance component estimation according to claim 1, characterized in that: In step S1, dual-orbit or multi-orbit time series SAR images covering the target area are collected, and these data are all from the same satellite, and the number of collected images is not less than two, and the two images are from different orbits; at the same time, external DEM data within the corresponding area is obtained, including SRTM DEM and TanDEM; In step S1, the multi-view time series SAR images of the same track are matched, including comparing the time baseline, spatial baseline and Doppler centroid frequency difference between the SAR images, selecting the SAR image corresponding to the maximum correlation coefficient as the common main image according to the comprehensive correlation function, and matching the remaining SAR images with the main image as slave images respectively, and the matching method adopts intensity cross-correlation matching; the intensity cross-correlation matching is mainly used to evaluate the similarity of two images by calculating the cross-correlation coefficient between them, so as to achieve image matching; specifically, by calculating the cross-correlation value of the two images at each position, finding the position with the maximum cross-correlation value, so as to determine the best match of the two images, the higher the cross-correlation coefficient, the more similar the two images are; for two images I and J, their cross-correlation function R(x, y) is defined as: Among them, I(i,j) and J(i+x,j+y) represent the pixel values ​​of the two images at corresponding positions; (i,j) represents the pixel position of the image, and (x,y) represents the displacement.

3. The method for automatically extracting 3D SAR control points based on variance component estimation according to claim 1, characterized in that: In step S2, when the number of images is greater than 15, the amplitude deviation index threshold method is used to detect permanent scatterers. Permanent scatterers correspond to high-reflection targets that are stable in time. Specifically, the amplitude deviation threshold method is used to calculate the ratio of the amplitude standard deviation to the average amplitude to identify stable scatterers; the amplitude deviation D A It is expressed as follows: Among them, σ φ represents the phase standard deviation; m A and σ A They represent the time average value of pixel amplitude, i.e. mathematical expectation, and amplitude standard deviation respectively; When the number of images is less than or equal to 15, the amplitude threshold method is used to detect high-reflection scatterers. The amplitude threshold method is used to find high-coherence point targets that show strong reflection in the amplitude sequence, which are permanent scatterers.

4. The method for automatically extracting 3D SAR control points based on variance component estimation according to claim 1, characterized in that: In step S3, based on the external DEM data, the SAR data is converted from the radar image coordinate system to the geographic coordinate system through the range-Doppler-ellipsoid equation to obtain the rough three-dimensional geographic coordinates of the permanent scatterer; the range-Doppler-ellipsoid equation in steps S3 and S4 is as follows: Among them, ct r represents the slant distance R from the ground target point P to the satellite; c is the speed of light; t r is the one-way distance propagation time; λ is the wavelength of the radar signal; m and n are the semi-major axis and semi-minor axis of the reference ellipsoid, respectively; h is the geodetic height of the ground target P; S S t a The satellite position vector at this moment; V S t a The satellite's velocity vector at this moment; S P t a The position of the ground target point P at the moment; f D is the Doppler frequency of the target point P; In t a At this moment, the satellite position vector is S S =(X S ,Y S ,Z S ) T , the velocity vector is The position of the ground target point P is S P =(X P ,Y P ,Z P ) T ;f D is the Doppler frequency of the target point P. After the SAR raw data is focused, the RD equation refers to the zero Doppler frequency, i.e., f D =0, at this time, the satellite's position vector and velocity vector are also the values ​​corresponding to zero Doppler frequency.

5. The method for automatically extracting 3D SAR control points based on variance component estimation according to claim 1, characterized in that: In step S5, the search range is expanded with the initial pixel of the same name as the center and a 10×10 rectangular box as the radius to extract multiple candidate pixels of the same name centered on the initial pixel of the same name; The azimuth time observation value t of each candidate pixel with the same name in step S6 a and the distance-to-time observation t r is the pixel coordinate (x pt ,y pt ) is calculated indirectly, and the calculation formula is as follows: Among them, t a0 Indicates the UTC time of the starting pixel of the azimuth; t r0 is the ratio of the round-trip distance from the starting pixel to the speed of light; Δt r =1 / rsr,Δt a =1 / prf, rsr is the range sampling rate, prf is the pulse repetition frequency; The above parameters can be found in the image parameter file.

6. The method for automatic extraction of 3D SAR control points based on variance component estimation according to claim 1, characterized in that: The point target analysis in step S6 includes the following two methods: One of them is to apply spectral zero padding to the image subset containing the target pixels, using a SINC interpolator to oversample the complex data with an oversampling factor equal to the number of inserted zeros; Another point target analysis method is to identify sub-pixel coordinates by calculating the centroid coordinates of the oversampled image; the specific process is as follows: first, the human eye identifies the corner reflector pixels in the SAR intensity image, then intercepts the sub-blocks of the n×n window containing the corner reflector pixels, and uses the SINC function or bilinear interpolation function to perform m-fold oversampling on the image sub-blocks, and finally calculates the centroid coordinates of each row and column of the oversampled image. and the average intensity value Among them, I and f represent the pixel intensity values ​​before and after oversampling respectively; (i, j, k) are pixel coordinates; i′, j′ = 1, 2, …, n; Finally, the center coordinates of the intensity peak of the image subset are calculated, that is, the sub-pixel position (x c ,y c ):

7. The method for automatically extracting 3D SAR control points based on variance component estimation according to claim 1, characterized in that: In step S7, the azimuth time observation value and the range time observation value of the candidate pixels with the same name are corrected by using external monitoring data to model the solid tide and atmospheric delay errors, thereby eliminating the influence of the solid tide and atmospheric delay errors on the SAR positioning accuracy; The local regional motion of the surface displacement caused by the solid tide in the east, north and vertical directions is obtained using open source software, and the three-dimensional deformation is projected into the azimuth and range directions: Among them, ξ az ,ξ rg They represent the delay in azimuth and range respectively; β represents the azimuth angle of satellite flight; θ is the local incident angle during imaging; ξ E ,ξ N ,ξ U denote the SET displacements in the east, north, and vertical directions, respectively; The azimuth and range corrections calculated based on the above formula can reduce the timing error caused by SET; Atmospheric delay mainly includes tropospheric delay and ionospheric delay. The classical geodetic method is used to estimate the tropospheric zenith delay above the InSAR monitoring point using the nearby GNSS to obtain the tropospheric delay of the local point and reduce the impact of tropospheric delay: Among them, I ZPD is the tropospheric zenith delay; θ is the local incident angle during imaging; h SAR Indicates the ground elevation of the SAR monitoring point; h GNSS Indicates the ground elevation of the GNSS monitoring point near the SAR monitoring point; h0 = 6000 meters; The ionosphere is distributed between 100 and 1500 kilometers above sea level. It is an atmospheric region where the atmosphere is partially ionized. The impact of the ionosphere is generally at the decimeter level. The vertical total electron content of the ionosphere provided by the IGS Analysis Center is used. VETC To compensate for the ionospheric delay in the range upwards: Where, f represents the radar carrier frequency, K = 40.28m 3 s -2 ; The VTEC grid map provided by the IGS Analysis Center has a spatial resolution of 2.5° and 5° in longitude and latitude, respectively, and a temporal resolution of 1 hour. VETC It is calculated by plane interpolation method, and the VTEC at the target point is interpolated through four adjacent grid points.

8. The method for automatically extracting 3D SAR control points based on variance component estimation according to claim 1, characterized in that: The adjustment function model for time observation values ​​and three-dimensional geographic coordinates established in step S8 includes: Based on the geometric relationship between the SAR satellite and the ground point target, high-precision stereo SAR positioning is achieved by the tilt intersection of the multi-track range-Doppler equation; for one of the SAR images, the time observation value and the three-dimensional coordinates of the unknown target point are implicitly included in the following two trajectory models: Here, for the target control point P, assume that its three-dimensional geographic coordinates are S P =(X P ,Y P ,Z P ) T , at time t in azimuth a At this moment, the satellite position vector is S S =(X S ,Y S ,Z S ) T , the velocity vector is c is the speed of light, t r is the one-way distance propagation time; S S and V S They represent the SAR sensor at the target imaging time t a The position and velocity vectors are calculated using the sixth-order polynomial model of the discrete satellite orbit state vector provided in the SAR parameter file. The specific calculation formula is as follows: Where: n represents an n-order polynomial, and the coefficients of the polynomial model are (a i ,b i ,c i ); these coefficients are estimated by the least squares method; In order to calculate the three-dimensional coordinates of the target point, more than two image acquisitions on different tracks are required.

9. The method for automatically extracting 3D SAR control points based on variance component estimation according to claim 8, characterized in that: In step S9, by linearizing formula (14) in step S8, the observation equation of the multi-view SAR image is expressed as a conditional adjustment with additional unknown parameters: Bv+Ax+w=0(16) Among them, B is the conditional matrix; v is the observed value residual; A represents the coefficient matrix of unknown parameters; x represents the coordinate correction of the target point, and w is the closure error; Since the coefficient matrix B is invertible, the above formula can be written as: v=B′xl′(17) Where B′=-B -1 A,l′=B -1 w; Equation (16) is converted into the classical least squares formula, and the coordinate correction of the target point is It is calculated according to the following formula: In the above formula, N represents the normal equation matrix of the coefficient matrix; T is the transposition process; P is the weight; By iteratively adding coordinate corrections to the initial coordinates The three-dimensional geographic coordinates of the ground target point corresponding to each candidate pixel with the same name are obtained; the variance-covariance matrix of its coordinate estimation parameters is as follows: in, represents the variance-covariance matrix of the coordinate estimation parameters; σ x ,σ y ,σ z is the coordinate standard deviation; σ xy ,σ xz ,σ yz represents the coordinate covariance; n is the number of observations; m is the number of unknown parameters.

10. The method for automatically extracting 3D SAR control points based on variance component estimation according to claim 1, characterized in that: In step S10, the coordinate accuracy after stereo SAR positioning is used as the judgment standard for pixels with the same name, and the standard deviation and Σ|σ are set. i The final homonymous pixels and their 3D geographic coordinates are determined by the threshold value of |, where the sum of the standard deviations of homonymous pixels is the smallest.

Images