A stereo SAR control point automatic extraction method based on variance component estimation

Abstract

The application provides a stereo SAR control point automatic extraction method based on variance component estimation. The method obtains double-or multi-orbit multi-scene time-series SAR images and DEM data of the same target area, registers the multi-scene time-series SAR images of the same orbit, extracts permanent scatterers on the images to perform geographic coding, obtains rough geographic coordinates of the permanent scatterers, inversely calculates initial homonym pixel coordinates of the permanent scatterers in another orbit, obtains multiple candidate homonym pixels and sub-pixel level coordinates of each candidate homonym pixel, corrects distance and azimuth time observations of the candidate homonym pixels, establishes a function model of adjustment about the time observations and three-dimensional geographic coordinates, iteratively solves the linearized function model through variance component estimation, and obtains a variance-covariance matrix of estimated parameters; and finally determines the homonym pixels and three-dimensional geographic coordinates based on a standard deviation threshold. The application improves the positioning and solving precision of stereo SAR control points, and can extract a large number of high-precision ground control points.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of image processing, and particularly relates to a stereo SAR control point automatic extraction method based on variance component estimation. BACKGROUND

[0002] Ground control points (GCPs) are key elements in photogrammetry, which have known spatial coordinate information. By marking the positions of GCPs in images, the association between images and real-world coordinates can be established, thereby realizing high-precision tracking, positioning and analysis of target objects. This precise spatial coordinate information is of great significance for surveying, remote sensing and three-dimensional modeling applications in aerial photogrammetry. In addition, GCPs can be used to correct images, improve the geometric accuracy of images and the quality of map projections, especially for synthetic aperture radar (SAR) imaging technology. However, because the signal transmission of SAR satellite systems is affected by satellite clock errors, orbit errors, atmospheric delays, and surface changes, the geometric positioning accuracy is reduced, so accurately positioning high-resolution SAR image monitoring points is a difficult and painful point at the present stage. Moreover, time-series InSAR (Interferometric Synthetic Aperture Radar) is a relative technology, and the height estimation error of monitoring points and reference point error are also the main error sources affecting geographic coding. Lower positioning accuracy makes it difficult to meet the fine monitoring needs of InSAR, which brings difficulties to the interpretation of deformation results.

[0003] Stereo SAR (SAR) is a technology that uses double / multi-angle SAR images with a certain intersection angle to retrieve three-dimensional information of the scene, which can be used to extract image ground control points. A prerequisite for stereo SAR positioning is that homonymous points must be visible from multi-view SAR images, such as light poles, corner points and intersection points. However, due to complex overlapping structures and different imaging angles, this type of homonymous points rarely occurs and is difficult to identify with the naked eye. Existing stereo SAR methods mainly include two parts: automatic detection of homonymous points and stereo SAR precise positioning. The extraction of homonymous points relies on external optical data and vectorized road grid data, which is difficult to automatically detect, and the positioning accuracy of stereo SAR solution is limited. SUMMARY

[0004] To solve the above problems, the application provides a stereo SAR control point automatic extraction method based on variance component estimation, which is based on traditional stereo SAR technology, uses additional parameter conditional adjustment and variance component estimation method to solve the three-dimensional geographic coordinates of the image homonym, provides a new idea for homonym detection of double / multi-angle SAR images, can automatically extract a large number of homonym pixels, provides a basis for SAR control point automatic extraction, and thus can solve the problems of difficult automatic detection of homonym and low stereo SAR positioning accuracy when ground control points are extracted by using double / multi-angle SAR images.

[0005] In order to achieve the above technical purposes, the application provides a stereo SAR control point automatic extraction method based on variance component estimation, which is characterized by comprising the following steps:

[0006] S1. Obtain double-orbit or multi-orbit multi-scene time series SAR images and DEM data of the same target area, and register the multi-scene time series SAR images of the same orbit;

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

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

[0009] S4. Use the rough geographic coordinates of the permanent scatterers on the first orbit obtained in step S3 to inversely calculate the SAR pixel coordinates of the permanent scatterers on another orbit by distance-Doppler-ellipsoid equation, that is, the initial homonym pixel coordinates;

[0010] S5. Expand the search range with the initial homonym pixel coordinates of the permanent scatterers on the first orbit on another orbit in step S4 as the center to extract multiple candidate homonym pixels of the permanent scatterers on the first orbit on another orbit;

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

[0012] S7. Correct the azimuth time observation value and the range time observation value of each candidate homonym pixel calculated in step S6 by using external monitoring data;

[0013] S8. For each candidate homonym pixel, use two or more SAR images to simultaneously solve multiple distance-Doppler equations to establish a adjustment function model about time observation value and three-dimensional geographic coordinates;

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

[0015] S10. Based on the three-dimensional coordinates of the ground target point of all candidate homonymous pixels obtained in step S9 and the coordinate variance-covariance matrix, determine the final homonymous pixel of the permanent scatterer on the orbit and its three-dimensional geographic coordinates based on a standard deviation threshold, and the ground point corresponding to the final homonymous pixel is the ground control point of the stereo SAR.

[0016] Further technical solutions of the present application: in step S1, the double-orbit or multi-orbit time series SAR images covering the target area are collected, these data are all satellite data, and the number of collected images is not less than two, and the two images come from different orbits; at the same time, external DEM data in the corresponding area range are obtained, including SRTM DEM and TanDEM;

[0017] In step S1, the multi-scene time series SAR images of the same orbit are registered, which includes selecting the SAR image corresponding to the maximum correlation coefficient as the common master image according to the integrated correlation function by comparing the time baseline, spatial baseline and Doppler centroid frequency difference between the SAR images, and registering the remaining SAR images as slave images with the master image, wherein the registration method adopts 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 the two images, so as to realize the registration of the images; specifically, the maximum position of the cross-correlation value is found by calculating the cross-correlation value of two images at each position, so as to determine the best matching of the two images, and the higher the cross-correlation coefficient, the more similar the two images; for two images I and J, the cross-correlation function R(x, y) is defined as:

[0018]

[0019] Wherein, I(i, j) and J(i+x, j+y) 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 amount.

[0021] The preferred technical solution of the present application: in step S2, when the number of images is greater than 15, 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, and specifically, the amplitude deviation threshold method is used to calculate the ratio of amplitude standard deviation to average amplitude to identify stable scatterers; the amplitude deviation DA is expressed by the following formula:

[0022]

[0023] wherein σ φ represents a phase standard deviation;

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

[0025] When the number of images is less than or equal to 15 scenes, the amplitude threshold method is used for high reflection scatterer detection, and the amplitude threshold method is used to find the high coherent point target with strong amplitude sequence representation, i.e. permanent scatterer.

[0026] The further technical scheme of the present application: in step S3, based on 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 specifically as follows:

[0027]

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

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

[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 t a ; V S is the velocity vector of the satellite at t a ;

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

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

[0034] At 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, and the R-D equation after focusing of the SAR raw data refers to zero Doppler frequency, that is, f D = 0, at this time, the position vector and the velocity vector of the satellite are also the values corresponding to zero Doppler frequency.

[0035] A further technical solution of the present application is that in step S5, the search range is expanded with the initial homonymous pixel as the center and a 10*10 rectangular frame as the radius, and a plurality of candidate homonymous pixels with the initial homonymous pixel as the center are extracted;

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

[0037]

[0038] Wherein, t a0 represents the world coordinate time of the azimuth direction starting pixel;

[0039] t r0 is the ratio time of the round-trip distance of the range direction starting pixel and 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] The above parameters can be found in the image parameter file.

[0042] A further technical solution of the present application is that in step S6, the point target analysis includes the following two ways:

[0043] One of them is to apply spectrum zero padding to the image subset containing the target pixel, and to oversample the complex data using a SINC interpolator, and the oversampling factor is equal to the number of inserted zeros;

[0044] Another point target analysis method is to identify sub-pixel coordinates by calculating the centroid coordinates of oversampled images; the specific process is as follows: first, the human eye identifies the corner reflector pixels in the SAR intensity image, then intercepts the n*n window sub-block containing the corner reflector pixels, and adopts SINC function or bilinear interpolation function for m times oversampling on the image sub-block, and finally calculates the centroid coordinates of each row and each column of the oversampled image and the average intensity value

[0045]

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

[0047] (i,j,k) is the pixel coordinate; i',j' = 1,2,…,n;

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

[0049]

[0050] The further technical scheme of the present application is that: the azimuth direction time observation value and the range direction time observation value of the correction candidate homonymic pixels in step S7 are specifically modeled by using external monitoring data to eliminate the influence of solid tide and atmospheric delay error on SAR positioning accuracy;

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

[0052]

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

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

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

[0056] Atmospheric delay mainly includes tropospheric delay and ionospheric delay; using classical geodetic method, tropospheric zenith delay above InSAR monitoring point is estimated by using adjacent GNSS, so as to obtain tropospheric delay of local point, and the influence of tropospheric delay is reduced:

[0057]

[0058] Wherein, I ZPD Tropospheric zenith delay; θ is local incidence angle when imaging;

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

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

[0061] Ionosphere is distributed between 100 to 1500 kilometers above sea level, and is the part of the atmosphere that is partially ionized, and the influence of ionosphere is generally in the order of decimeter; the ionospheric vertical total electron content map I VETC Provided by IGS analysis center is used to compensate the delay of ionosphere in the distance direction:

[0062]

[0063] Wherein, f indicates radar carrier frequency, K=40.28m 3 s -2 ;

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

[0065] Further technical scheme of the application: the step S8 includes:

[0066] According to the geometric relationship between SAR satellite and ground point target, high-precision stereo SAR positioning is realized by the tilt intersection of multi-orbit distance-doppler equation; for one SAR image, time observation value and unknown target point three-dimensional coordinates are implicitly contained in the following two track models:

[0067]

[0068] Here, for target control point P, assume that its three-dimensional geographic coordinates are S P =(XP ,Y P Z P ) T In the direction of time t a At time S, the satellite's position vector is S =(X S ,Y S Z S ) T The velocity vector is c is the speed of light, t r It is the one-way distance propagation time; S S and V S These represent the SAR sensor's position at target imaging time t, respectively. a The position and velocity vectors are calculated using a sixth-order polynomial model of the discrete satellite orbit state vectors provided in the SAR parameter file. The specific calculation formula is as follows:

[0069]

[0070] Where: n represents an n-degree polynomial, and the coefficients of the polynomial model are (a i ,b i ,c i These coefficients are estimated using the least squares method.

[0071] To calculate the three-dimensional coordinates of the target point, image acquisition from two or more different orbits is required.

[0072] A further technical solution of the present invention: In step S9, by linearizing formula (14) in step S8, the observation equation of multiple SAR images is expressed as a conditional adjustment with additional unknown parameters:

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

[0074] Where B is the condition matrix; v is the observation residual; and A represents the coefficient matrix of the unknown parameters.

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

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

[0077] v=B′xl′(17)

[0078] Where, B′=-B -1 A,l′=B -1 w;

[0079] Equation (16) is transformed into the classical least squares formula, and the coordinate correction of the target point is... The specific calculation is based on the following formula:

[0080]

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

[0082] By adding the coordinate correction amount to the initial coordinate iteration Get the three-dimensional geographic coordinates of the ground target point corresponding to each candidate homonymic pixel;The variance-covariance matrix of its coordinate estimation parameter is as follows:

[0083]

[0084] Wherein, Indicates the variance-covariance matrix of the coordinate estimation parameter;

[0085] σ x ,σ y ,σ z Is the coordinate standard deviation;σ xy ,σ xz ,σ yz Indicates the coordinate covariance;

[0086] N is the number of observation values;M represents the number of unknown parameters.

[0087] The further technical scheme of the application: the coordinate internal accuracy after stereoscopic SAR positioning in the step S10 is used as the judgment standard of the homonymic pixel, and the threshold value of the standard deviation and ∑|σ i | is used to determine the final homonymic pixel and its three-dimensional geographic coordinates, wherein the standard deviation of the homonymic pixel is the smallest.

[0088] The stereoscopic SAR can perform absolute three-dimensional positioning on the natural permanent scatterer, which makes it possible to generate ground control points only using SAR data. The prerequisite for obtaining high-precision GCPs by this method is to correctly detect the homologous scatterers in SAR images with different observation geometries. In the application, a processing strategy for automatically detecting the same target in SAR images of the same area taken from different orbits is described, and in addition, a complete GCP automatic generation process based on stereoscopic SAR is provided, which can extract a large number of high-precision ground control points and has high engineering application value.

[0089] The application adopts the conditional adjustment of additional parameters and the variance component estimation method on the basis of the traditional stereo SAR technology, and calculates the three-dimensional geographic coordinates of the image homonymic points, thereby providing a new idea for the homonymic point detection of the dual / multi-angle SAR image, and automatically extracting a large number of homonymic pixels, thereby providing a basis for the automatic extraction of the SAR control points. Compared with the existing stereo SAR control point generation technology, the stereo SAR control point automatic extraction method of the application can automatically detect the SAR image homonymic points without the external optical image and the vectorized road grid data; in addition, the application improves the positioning calculation precision of the stereo SAR control point by the point target analysis, the time observation value correction, the function adjustment and the variance component estimation. In general, the application provides a complete stereo SAR control point automatic generation process, can extract a large number of high-precision ground control points, and has high engineering application value. BRIEF DESCRIPTION OF DRAWINGS

[0090] Figure 1 is a flow chart of the application;

[0091] Figure 2 is a search schematic diagram of the candidate homonymic pixels in the embodiment of the application;

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

[0093] Figure 4-a is a propagation error schematic diagram obtained by projecting the three-dimensional ground surface displacement caused by the solid tide at different times to the azimuth direction and the range direction;

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

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

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

[0097] Figure 5-b is a final homonymic pixel schematic diagram on the corresponding another orbit in the embodiment. DETAILED DESCRIPTION

[0098] It should be noted that the embodiments in the embodiments and the features in the embodiments can be combined with each other without conflict, and the technical solutions in the embodiments will be described clearly and completely in combination with the drawings of the embodiments of the application. Obviously, the described embodiments are only part of the embodiments of the application, not all. Based on the embodiments in the application, all other embodiments obtained by those skilled in the art without creative labor are within the scope of protection of the application.

[0099] The embodiments provide a stereo SAR control point automatic extraction method based on variance component estimation, as shown in Figure 1 The method comprises the following steps:

[0100] S1. According to actual needs, obtain double-orbit or multi-orbit multi-scene time series SAR images of the same target area, DEM data and its parameter file, the number of collected SAR images is not less than two scenes, and the two scenes are from different orbits. In general, the greater the difference in imaging angle, the more stable the spatial intersection, and the higher the positioning accuracy. In addition, collect external DEM data in the corresponding area, such as SRTM DEM and TanDEM, and the data resolution is determined according to the cost and demand. For multi-scene time series SAR images of the same orbit, in order to enable multi-source images to be effectively compared, fused and analyzed, the same orbit images need to be registered 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 realize the registration of the images. The basic principle of intensity cross-correlation registration is to find the position with the maximum cross-correlation value by calculating the cross-correlation value of two images at each position, so as to determine the best matching of two images. The cross-correlation coefficient reflects the similarity of two images in space, and the higher the value, the more similar the two images. 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) represent the pixel values of 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. Extract permanent scatterers on multi-scene time series SAR images of different orbits by using amplitude dispersion or amplitude threshold method;

[0106] When the number of images is large, generally greater than 15, the amplitude deviation index threshold method can be used to detect permanent scatterers, which usually correspond to time-stable high-reflective targets such as lamp posts, corner points, etc. 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 is expressed by the following formula:

[0107]

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

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

[0110] When the number of images is small, generally less than or equal to 15, the amplitude threshold method can be used to detect high-reflective scatterers. The amplitude values of each SAR image are calculated to obtain the amplitude mean value. A threshold is set according to the amplitude mean value, and the minimum value of the amplitude mean value is usually selected as the threshold. The pixel points with amplitude values higher than the threshold are identified as permanent scatterers.

[0111] S3. Geocoding the permanent scatterers on the track extracted in step S2 based on DEM data to obtain the rough geographic coordinates of each permanent scatterer; Geocoding means determining the three-dimensional geographic coordinates of the SAR image pixels. According to the time sequence and speed of the radar flight in the azimuth direction during SAR imaging, as well as the echo time delay and speed of light received in the range direction, the two-dimensional geometric coordinates (a, r) of the radar image can be determined. Since the radar beam center intersects the earth's surface to determine the actual geographic 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 (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 propagation time; λ is the radar signal wavelength;

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

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

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

[0118] In t a At time S, the satellite's position vector is S =(X S ,Y S Z S ) T The velocity vector is The location of ground target point P is S P =(X P ,Y P Z P ) T f D Let f be the Doppler frequency of target point P. After focusing, the RD equation of the raw SAR data references the zero Doppler frequency, i.e., f. D =0. At this point, the satellite's position and velocity vectors also correspond to the zero Doppler frequency. It can be seen that azimuth and range timing errors, as well as external DEM errors, are the three factors affecting the accuracy of RD equation positioning. Especially for high-resolution SAR imagery, low-resolution DEMs reduce geocoding accuracy. Since time-series InSAR is a relative measurement technique, reference point errors are also a major source of error affecting positioning accuracy. Stereo SAR positioning only considers the range-Doppler equations for different orbits, effectively avoiding the influence of external DEM errors and providing more accurate ground coordinates.

[0119] S4. Using the rough geographic coordinates of the permanent scatterer on orbit one obtained in step S3, the SAR pixel coordinates of the permanent scatterer on another orbit are calculated by using the distance-Doppler-ellipsoid equation (RD equation). These pixels are regarded as the initial homonymous pixels between different orbits.

[0120] S5. Expand the search range centered on the initial pixel coordinates of the permanent scatterer on track 1 on another track, using a 10×10 rectangle as the radius, and extract multiple candidate pixels with the same name from the permanent scatterer on track 1 on the other track within the rectangle; see details below. Figure 2 ,in Figure 2 (a) represents a permanent scatterer on orbit 1. Figure 2 (b) represents multiple candidate pixels with the same name on another track. Figure 2(c) represents the standard deviation value corresponding to each candidate co-name pixel pair;

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

[0122] Referring to Figure 3 For each candidate co-name pixel pair, the point target analysis is used to extract the sub-pixel coordinates, thereby improving the stereoscopic SAR positioning precision; wherein, Figure 3 The left represents the original pixel coordinates of the permanent scatterer and the corresponding post-scattering intensity value, Figure 3 The right represents the sub-pixel coordinates of the permanent scatterer and the corresponding post-scattering intensity value.

[0123] For the existing SAR satellite imaging mode, the pre-processed first-level single-view complex image is converted into a two-dimensional image grid with distance and azimuth radar coordinates, each grid representing a pixel, and the grid coordinates corresponding to the pixel coordinates. The azimuth time observation value t a and the range time observation value t r of each candidate co-name pixel are indirectly calculated through the pixel coordinates (x pt , y pt ), and the calculation formula is as follows:

[0124]

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

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

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

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

[0129] Assuming the given image resolution is 1 m and the SAR measurement accuracy is targeted at 1 cm level, the point target analysis can provide the peak position with about 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 using SINC interpolator to oversample the complex data, and the oversampling factor is equal to the number of inserted zeros, but a large amount of padding will go through the 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 pixel in the SAR intensity image by human eyes, then intercept the n*n window sub-block containing the corner reflector pixel, and use SINC function or bilinear interpolation function to oversample m times for the image sub-block, and 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) is the pixel coordinate; i',j' = 1,2,…,n;

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

[0134]

[0135] S7. Correct the azimuth time observation value and the range time observation value of each candidate homonym pixel calculated in step S6 by using external monitoring data;

[0136] In order to correct the errors of the azimuth time observation value and the range time observation value, the present application models the solid tide and the atmospheric delay error by using external monitoring data, so as to eliminate the influence of these factors and improve the stereoscopic SAR positioning accuracy. Referring to Figures 4-a to 4-c , Figure 4-a represents the three-dimensional surface displacement caused by the solid tide at different times, and the propagation error can be obtained by projecting it to the azimuth direction and the range direction, 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 achievable SAR observation accuracy currently mainly needs to correct the solid earth tides (SET). SET is the deformation of the solid earth surface caused by the gravity of the sun and the moon. The periodic deformation caused by SET within a day can reach several decimeters. The local area motion of the ground surface displacement caused by SET in the east, north and vertical directions can be obtained by using open source software, and the three-dimensional deformation is projected to the azimuth and range directions:

[0138]

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

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

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

[0142] The azimuth and range corrections calculated based on the above formula can weaken the timing errors caused by SET.

[0143] Atmospheric delay mainly includes tropospheric delay and ionospheric delay. Tropospheric delay is the layer of atmosphere closest to the earth's surface, which is related to the terrain height and is located in the lowermost layer of the atmosphere. The troposphere has the maximum density in the atmosphere, and its influence on positioning accuracy is also the most serious, and the error produced is even as high as several meters. In order to reduce the influence of tropospheric delay, the classical geodetic method can be used to estimate the zenith path delay (ZPD) of the troposphere above the InSAR monitoring point by using the GNSS close to the monitoring point, so as to obtain the tropospheric delay of the local point:

[0144]

[0145] wherein I ZPD is the zenith path delay of the troposphere; θ is the local incidence angle at 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; h0=6000 meters.

[0148] This method assumes that the troposphere within a certain range has spatial correlation, so the premise of calculating the tropospheric delay by 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 deviation of the estimated tropospheric delay. After years of research, more advanced tropospheric delay estimation methods have emerged. When there is no GNSS station near the study area, the ZPD map provided by the GACOS can be used to calculate the tropospheric delay. The advantage of the GACOS tropospheric map is global availability, near real-time and high spatial resolution.

[0149] The ionosphere is located between 100 and 1500 kilometers above sea level, and is a region of the atmosphere that is partially ionized. The impact of the ionosphere is generally on the order of decimeters. Here, the ionosphere vertical total electron content (VTEC) map provided by the IGS analysis center (University of Bern, Switzerland) is used to compensate for the delay of the ionosphere in the range direction: VETC

[0150]

[0151] where f represents the radar carrier frequency, K = 40.28 m 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 time resolution of 1 hour. The VTEC at the SAR monitoring point is calculated by a planar interpolation method, and the VTEC at the target point position is interpolated by four adjacent grid points. VETC

[0153] S8. For each candidate homonymous pixel, two or more SAR images are used to solve multiple range-Doppler equations to establish a function model of the time observation value and three-dimensional geographic coordinates; according to the geometric relationship between the SAR satellite and the ground point target, high-precision stereo SAR positioning can be realized by the inclined intersection of multiple track range-Doppler equations. For one SAR image, the time observation value and the unknown three-dimensional coordinates of the target point are implicitly contained in the two track models:

[0154]

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

[0156]

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

[0158] In order to calculate the three-dimensional coordinates of the target point, more than twice image acquisition of different orbits is required. Two orbits provide four equations, and four equations can solve three unknown coordinate parameters. The more the number of images, the more the number of equations, and the more robust the parameter solution.

[0159] S9. Linearize the adjustment function model established in step S8, and iteratively solve the linearized adjustment function model through variance component estimation, to obtain the three-dimensional coordinates of the ground target point corresponding to each candidate homonymous pixel based on the iterative weighted least squares algorithm, and estimate the coordinate variance-covariance matrix; the specific process is as follows: through linearization equation (14), the observation equation of the multi-scene SAR image can be represented as a conditional adjustment with additional unknown parameters:

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

[0161] where B is the condition matrix; v is the observation value 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 formula 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 a classic 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 is the transpose processing, and P is the weight.

[0169] By iterating the initial coordinates and adding the coordinate correction The three-dimensional geographic coordinates of the target point are obtained.

[0170] One special advantage of least squares estimation is to provide complete error statistics for the estimation result. One of the quality indicators used in the present application 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] wherein, represents the variance-covariance matrix of the coordinate estimation parameters, which is used to describe the internal accuracy of the coordinate estimation of the target point;

[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 errors of the two types of observations cannot be accurately known, the weight P is determined iteratively using the variance component estimation method. First, a reasonable prior mean error of the range-time and azimuth-time observations is given, and the first weight is determined, such as assuming that the initial weight of the range-time observation is The initial weight of the azimuth-time observation is After the first calculation, the unit weight variance of the two types of observations is estimated using the residual

[0176]

[0177] wherein, v ρ and v aThese represent the corrections for the range and azimuth observations, respectively; p ρ and p a Assign their weights;

[0178] n ρ and n a The number of observations in both directions; trace represents the trace of the matrix.

[0179] Reweighting based on the estimated unit weight variance:

[0180]

[0181] Where k represents the number of iterations.

[0182] Iterate through the adjustment and weighting steps until... If the difference between the two is less than a predetermined threshold, the iteration operation stops, and the estimated 3D geographic coordinates of the candidate pixels with the same name are output.

[0183] S10. Based on the 3D coordinates and coordinate variance-covariance matrix of all candidate homonymous pixels obtained in step S9, determine the final homonymous pixels and their 3D geographic coordinates of the permanent scatterer on orbit 1 based on the standard deviation threshold. The ground point corresponding to the final homonymous pixel is the ground control point of the stereo SAR. (See details...) Figure 5-a and Figure 5-b The final corresponding pixel on orbital 5-a is shown. Figure 5-b For the final pixel with the same name on another track.

[0184] The standard deviation threshold represents the standard deviation of the coordinates of a target pixel and its candidate counterparts after stereo SAR localization, and the threshold value is used. Experiments have shown that the standard deviation of the true counterparts is the smallest. Therefore, this invention uses the in-coordinate matching accuracy after stereo SAR localization as the criterion for judging counterparts:

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

[0186] Where i = x, y, z. By setting the standard deviation and ∑|σ i The threshold (e.g., set to 3) can determine the final pixel with the same name and its three-dimensional geographic coordinates. The target point corresponding to the pixel with the same name is regarded as the ideal ground control point.

[0187] The above description is merely one embodiment of the present invention, and while it is detailed and specific, it should not be construed as limiting the scope of the invention. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these modifications and improvements all fall within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the appended claims.

Claims

1. An automatic method for extracting control points in stereo SAR based on variance component estimation, characterized in that, Includes the following steps: S1. Acquire dual-track or multi-track time-series SAR images and DEM data of the same target area, and register the multiple time-series SAR images of the same track; S2. Permanent scatterers in multiple time-series SAR images with different orbits are extracted using amplitude deviation or amplitude thresholding methods; S3. Based on the DEM data, geocode the permanent scatterers on orbit 1 extracted in step S2 to obtain the approximate geographic coordinates of each permanent scatterer; S4. Using the rough geographic coordinates of the permanent scatterer on orbit 1 obtained in step S3, the SAR pixel coordinates of the permanent scatterer on another orbit are calculated in reverse using the distance-Doppler-ellipsoid equation, i.e., the initial pixel coordinates of the same name. S5. Expand the search range with the initial corresponding pixel coordinates of the permanent scatterer on track 1 on another track as the center in step S4, and extract multiple candidate corresponding pixels of the permanent scatterer on track 1 on another track. S6. Based on point target analysis, obtain the sub-pixel level coordinates of each candidate pixel with the same name, and calculate the azimuth time observation and range time observation of the corresponding sub-pixel level coordinates in turn. S7. Use external monitoring data to correct the azimuth time observation and range time observation of each candidate pixel with the same name calculated in step S6; S8. For each candidate pixel with the same name, use two or more SAR images to establish multiple range-Doppler equations and build an adjustment function model for time observations and three-dimensional geographic coordinates. S9. Linearize the adjustment function model established in step S8, and iteratively solve the linearized adjustment function model by estimating the variance components. Based on the iterative weighted least squares algorithm, obtain the three-dimensional coordinates of the ground target point corresponding to each candidate pixel with the same name, and estimate the coordinate variance-covariance matrix. S10. Based on the three-dimensional coordinates and coordinate variance-covariance matrix of all candidate synodic pixels obtained in step S9, determine the final synodic pixels and their three-dimensional geographic coordinates of the permanent scatterer on orbit 1 based on the standard deviation threshold, and the ground point corresponding to the final synodic pixel is the ground control point of the stereo SAR.

2. The method for automatic extraction of stereo 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 acquired. These images are all from the same satellite, and the number of acquired images is no less than two, with the two images coming from different orbits. At the same time, external DEM data within the corresponding area are acquired, including SRTM DEM and TanDEM. Step S1 involves matching multiple time-series SAR images from the same orbit. This includes comparing the temporal baseline, spatial baseline, and Doppler centroid frequency difference between the SAR images, and selecting the SAR image with the highest correlation coefficient as the common master image based on the comprehensive correlation function. The remaining SAR images are then used as slave images, and intensity cross-correlation registration is performed between them and the master image. Intensity cross-correlation registration is primarily used to assess the similarity between two images by calculating their cross-correlation coefficient, thereby achieving image registration. Specifically, it involves calculating the cross-correlation value of the two images at each location, finding the location with the highest cross-correlation value, and thus determining the best match between the two images. A higher cross-correlation coefficient indicates greater similarity between the two images. For two images I and J, their cross-correlation function R(x,y) is defined as: Where I(i,j) and J(i+x,j+y) represent the pixel values ​​of the two images at corresponding positions, respectively; (i,j) represents the pixel position of the image, and (x,y) represents the displacement.

3. The method for automatic extraction of stereo 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 reflectivity targets that are stable over 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 D A Expressed as follows: Where, σ φ Indicates the phase standard deviation; m A and σ A These represent the pixel's amplitude over time, i.e., the mathematical expectation and the amplitude standard deviation, respectively. When the number of images is less than or equal to 15, the amplitude thresholding method is used to detect high reflectivity scatterers. The amplitude thresholding method is used to find high coherence point targets that exhibit strong reflectivity in the amplitude sequence, which are permanent scatterers.

4. The method for automatic extraction of stereo SAR control points based on variance component estimation according to claim 1, characterized in that: In step S3, based on external DEM data, the SAR data is transformed from the radar image coordinate system to the geographic coordinate system using 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 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 radar signal wavelength. m and n are the semi-major and semi-minor axes of the reference ellipsoid, respectively; h is the geodetic height of the ground target P; S S For t a The position vector of the satellite at any given moment; V S For t a The velocity vector of the satellite at any given moment; S P For t a The position of ground target point P at any given time; f D The Doppler frequency of the target point P; In t a At time S, the satellite's position vector is S =(X S ,Y S Z S ) T The velocity vector is The location of ground target point P is S P =(X P ,Y P Z P ) T ;f D The Doppler frequency of the target point P, i.e., f D =0, at which point the satellite's position vector and velocity vector are also the values ​​corresponding to the zero Doppler frequency.

5. The method for automatic extraction of stereo 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 homonym pixel as the center and a 10×10 rectangle as the radius to extract multiple candidate homonyms centered on the initial homonym pixel. In step S6, the azimuth time observation value t for each candidate pixel with the same name a and distance to time observation t r By pixel coordinates (x) pt ,y pt It is obtained indirectly, and its calculation formula is as follows: Among them, t a0 Indicates the world-coordinated time (UTC) from the azimuth to the starting pixel; t r0 Time 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, where rsr is the distance sampling rate and prr is the pulse repetition frequency.

6. The method for automatic extraction of stereo 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 approach is to apply spectral zero-padding to a subset of the image containing the target pixels, using a SINC interpolator to oversample the complex data, with the oversampling factor equal to the number of zeros inserted. Another point target analysis method identifies sub-pixel coordinates by calculating the centroid coordinates of oversampled images. The specific process is as follows: First, corner reflector pixels are visually identified in the SAR intensity image. Then, an n×n window containing the corner reflector pixels is extracted as a sub-block. This sub-block is then oversampled m times using the SINC function or bilinear interpolation function. Finally, the centroid coordinates of each row and column of the oversampled image are calculated. and average strength value Where 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 coordinates of the intensity peak center of the image subset were calculated, i.e., the sub-pixel positions of the main scatterer in the range and azimuth directions (x, y, y). c ,y c ):

7. The method for automatic extraction of stereo SAR control points based on variance component estimation according to claim 1, characterized in that: In step S7, correcting the azimuth and range time observations of candidate pixels with the same name specifically involves using external monitoring data to model solid tide and atmospheric delay errors, thereby eliminating the impact of solid tide and atmospheric delay errors on SAR positioning accuracy. Localized surface displacements caused by solid tides in the east, north, and vertical directions were obtained using open-source software, and the three-dimensional deformations were projected onto the azimuth and range directions: Where, ξ az ,ξ rg These represent the azimuth and range delays, respectively; β represents the azimuth angle of the satellite's flight. θ is the local incident angle during imaging; ξ E ,ξ N ,ξ U These represent the solid tide SET displacements in the east, north, and vertical directions, respectively; The azimuth and range delays calculated based on the above formulas can reduce the timing errors caused by SET. Atmospheric delay mainly includes tropospheric delay and ionospheric delay. Using classical geodetic methods, the tropospheric zenith delay above the InSAR monitoring point is estimated using nearby GNSS data, thus obtaining the local tropospheric delay and reducing its impact. The local tropospheric delay ξ... tro The calculation is as follows: Among them, I ZPD The tropospheric zenith delay is θ; the local incident angle during imaging is θ. h SAR Indicates the ground elevation of the SAR monitoring point; h GNSs This represents the ground elevation of the GNSS monitoring point near the SAR monitoring point; h0 = 6000 meters; The ionosphere, located 100 to 1500 kilometers above sea level, is a partially ionized region of the atmosphere. The effects of the ionosphere are generally on the decimeter level. The vertical total electron content of the ionosphere is shown in Figure I provided by the IGS analysis center. VETC To compensate for the ionospheric delay in the distance direction: Where f represents the radar carrier frequency, and 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 time resolution of 1 hour. The I... VETC It is calculated using a planar interpolation method, which interpolates the VTEC of the target point position using four neighboring grid points.

8. The method for automatic extraction of stereo SAR control points based on variance component estimation according to claim 1, characterized in that: The step S8, which establishes the adjustment function model for time observations and three-dimensional geographic coordinates, includes: Based on the geometric relationship between SAR satellites and ground point targets, high-precision stereo SAR positioning is achieved through the tilted intersection of multiple orbital range-Doppler equations. For a single SAR image, the temporal observations and the three-dimensional coordinates of the unknown target point are implicitly contained in the following two trajectory models: Here, for the target control point P, let its three-dimensional geographic coordinates be S. P =(X P ,Y P Z P ) T In the direction of time t a At time S, the satellite's position vector is s =(X s ,Y s Z s ) T The velocity vector is c is the speed of light, t r It is the one-way distance propagation time; S s and V s These represent the SAR sensor's position at target imaging time t, respectively. a The position and velocity vectors are calculated using a sixth-order polynomial model of the discrete satellite orbit state vectors provided in the SAR parameter file. The specific calculation formula is as follows: Where: n represents an nth-degree polynomial, and the coefficients of the polynomial model are (a i ,b i ,c i These coefficients are estimated using the least squares method. To calculate the three-dimensional coordinates of the target point, image acquisition from two or more different orbits is required.

9. The method for automatic extraction of stereo 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 equations for multiple SAR images are expressed as conditional adjustment with additional unknown parameters: Bv + Ax + w = ​​0 (16) Where B is the condition matrix; v is the observation residual; and A represents the coefficient matrix of the 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 equation can be written as: v=B′xl′(17) Where, B′=-B -1 A,l′=B -1 w; Equation (16) is transformed into the classical least squares formula, and the coordinate correction of the target point is... The specific calculation is based on the following formula: In the above formula, N represents the normal equation matrix of the coefficient matrix; T is the transpose; and P is the weight. By iteratively adding coordinate corrections to the initial coordinates The 3D 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, This represents the variance-covariance matrix of the coordinate estimation parameters; σ x ,σ y ,σ z σ represents the standard deviation of the coordinates. xy ,σ xz ,σ yz Indicates coordinate covariance; n obs m is the number of observations. unk This indicates the number of unknown coordinate parameters, which is 3 in this case.

10. The method for automatic extraction of stereo 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 criterion for judging pixels with the same name. This is achieved by setting the standard deviation and ∑|σ i The threshold is used to determine the final corresponding pixel and its 3D geographic coordinates, where the standard deviation of the corresponding pixel is minimized.

Images