CN106526593B - Sub-pixel-level corner reflector automatic positioning method based on the tight imaging model of SAR - Google Patents

Abstract

The invention discloses a kind of sub-pixel-level corner reflector automatic positioning methods based on the tight imaging model of SAR, it specifically includes acquisition initial data, the initial pixel coordinate of calculating corner reflector point, the Pixel-level coordinate of resolving corner reflector point and calibration excluding gross error is measured by interferometric parameter, obtain the sub-pixel-level coordinate of accurate corner reflector point.The present invention has taken into account the image feature and geometrical characteristic of corner reflector, automation sub-pixel-level corner reflector can be achieved to extract, the method increase point extraction efficiencies, while ensure that the precision of position coordinate, realize the extraction corner reflector point coordinate of sub-pixel-level automation.

Description

Sub-pixel-level corner reflector automatic positioning method based on the tight imaging model of SAR

Technical field

The present invention relates to radar interference processing technology field, in particular to a kind of sub- picture based on the tight imaging model of SAR Plain grade corner reflector automatic positioning method.

Background technique

Synthetic aperture radar (SAR) technology is born in nineteen fifty-one, it is by collecting and recording target, background in microwave band Interior electromagenetic wave radiation, scattering energy extract the useful of identification target object or phenomenon using signal processing technology processing and sorting Information is widely used in surveying and drawing terrain and its features, detection ocean, detection earth resource and reconnaissance and surveillance military target etc..Corner reflection Device (CR) has the characteristics that high reflection echo, stiff stability in Synthetic Aperture Radar images, has been applied to the geometry of SAR image In correction and synthetic aperture radar interferometry (InSAR) calibration.

The selection such as corner reflector face, design, machining tolerance, angle are concentrated mainly on to the design work of corner reflector at present Error, incidence angle alignment error etc. are spent, but is related to the less of automatic positioning to corner reflector, especially in the prior art not There is the automatic positioning of the reflector of sub-pixel-level.Existing reflector location technology is mostly using corner reflector point with stronger Echo-signal and surrounding atural object have relatively weak echo-signal characteristic, using a certain size template window go detection, With image, corresponding related coefficient is calculated, the high point of related coefficient is the position where corner reflector, and this method calculates speed Degree is very fast, however its detection accuracy can only arrive Pixel-level.In addition, there are also identify angle based on the method for extracting PS point in the prior art Reflector points calculate coherence factor around corner reflector point in certain window, amplitude deviation value, and there is high signal-to-noise ratio point to be For corner reflector point position, this method is better than correlation coefficient process, and precision is can be controlled within 1-2 pixel, however image number Easily tens of scapes, cost are excessively high；In addition, there are also the corner reflector extracting method for being based on signal-to-noise ratio (SCR), the party in the prior art Method estimates the peak excursion of initial alignment first, determines the approximate range of triplane target, then draws determining corner reflection The SCR value figure of device target zone, the maximum point of SCR value is triplane target.However, these three methods of the prior art by Quality evaluation and mode is checked in lacking, therefore the quality of positioning result can not be assessed, once there are rough error, in high-precision In mapping, biggish resolution error can be introduced.

Influence of the SAR image vulnerable to speckle noise, if used template window and practical corner reflector image feature It is inconsistent, then identify, the precision that detects it is not high；In addition, position, incidence angle, processing size etc. that corner reflector is put all can be right SAR image feature generates apparent influence, for example the incidence angle of radar sensor deviates ideal position, then three faces of corner reflector Corner effect is unobvious, and the corner reflector feature actually obtained is unpredictable, and existing detection corner reflector algorithm can not be known Not, it either still goes to extract position using the method for PS point using template matching method, detection accuracy is not also high, therefore Without versatility, and positioning accuracy is unable to reach the precision of sub-pixel-level.In addition, extracting corner reflector using template window Point can not effectively provide the extraction accuracy of corner reflector.Corner reflector point identification to or not, precision Gao Bugao, Wu Faping Valence.Interference of certain corner reflectors by strong backscattering echo signal around can easily lead to positioning mistake, to influence Using effect of corner reflector during high-acruracy survey calibration.

It is, therefore, desirable to provide a kind of novel reflector automatic positioning method that can overcome disadvantages mentioned above.

Summary of the invention

The present invention proposes a kind of sub-pixel-level corner reflector automatic positioning method based on the tight imaging model of SAR, this hair It is bright use strict SAR imaging model carry out corner reflector automatically extract, to avoid in artificial participation process poor efficiency and Low precision problem.Meanwhile the defect of intuitive quality evaluating method, this hair can not be provided for corner reflector automatic Extraction Algorithm It is bright that point-by-point error is provided, and can excluding gross error point, improve the robustness that point extracts result, can fully eliminate due to existing skill One or more problem caused by the limitation and defect of art.

Additional advantages, objects, and features of the invention, a part will be set forth in the following description, and another portion Divide will be apparent or from implementation of the invention by the investigation to following explanation for those of ordinary skill in the art It acquires.The present invention can be realized and obtained by the structure particularly pointed out in the specification and claims of text and attached drawing Objects and advantages.

The present invention provides a kind of sub-pixel-level corner reflector automatic positioning method based on the tight imaging model of SAR, Be characterized in that, the method specifically includes the following steps:

S1 acquires initial data, and the initial data includes coordinate, main image and the main image of corner reflector point Imaging parameters；

S2 constructs the tight imaging model of SAR, and the tight imaging model of SAR based on building using the initial data of acquisition Calculate the initial pixel coordinate (x of all corner reflector points₀,y₀), obtain the set of the initial pixel coordinate of all corner reflector points {P_init}；

S3, the initial pixel coordinate set { P based on the step S2 corner reflector point obtained_init, it is calculated using peak swing Method resolves Pixel-level coordinate of all corner reflector points in main image, and obtains picture of all corner reflector points in main image Plain grade coordinate set { (P_pix)}；

S4, Pixel-level the coordinate set { (P of diagonal reflector points_pix) in all corner reflector points Pixel-level coordinate into The measurement verification of row interferometric parameter, obtains the sub-pixel-level coordinate of corner reflector point.

The present invention also provides a kind of computer-readable medium, there are computer program, the computer program in medium For executing method of the present invention.

The image feature and geometrical characteristic that the present invention has taken into account corner reflector are, it can be achieved that automation sub-pixel-level corner reflection Device extracts, and the method increase point extraction efficiencies, while ensure that the precision of position coordinate, realizes sub-pixel-level automation Extraction corner reflector point coordinate.

Detailed description of the invention

Fig. 1 is to be automatically positioned according to the sub-pixel-level corner reflector based on the tight imaging model of SAR of the embodiment of the present invention The flow chart of method.

Fig. 2 is according to landform three-dimensional reconstruction geometrical relationship under the earth's core rectangular coordinate system in space of the embodiment of the present invention.

Specific embodiment

With reference to the accompanying drawings to invention is more fully described, wherein illustrating exemplary embodiment of the present invention.

The tight imaging model of SAR utilizes the three-dimensional of range equation, Doppler equation and interferometric phase equation solution target point Coordinate describes the relationship between interference SAR fundamental measurement value and ground target point position, has merged related ground target position All information, be generally acknowledged tight imaging model.Sub-pixel-level angle proposed by the present invention based on the tight imaging model of SAR is anti- Emitter automatic positioning method can comprehensively consider the back scattering feature and interference geometrical characteristic of image, obtain corner reflector Be accurately positioned as a result, and provide the position error of each corner reflector, auxiliary carries out the screening of corner reflector, is further ensured that dry Relate to the validity of measurement calibration scheme.

As shown in Figure 1, the sub-pixel-level corner reflector automatic positioning method based on the tight imaging model of SAR, can realize The high-precision automatic of corner reflector extracts.This method specifically includes the following steps:

Step S1, acquires initial data, the initial data include the coordinate of corner reflector point, main image, from image with And the main image and the imaging parameters from image.

According to an embodiment of the invention, can be adapted for any spaceborne/airborne imagery；Preferably, implementation case of the invention Example is by taking airborne Ka wave band interference SAR image as an example, and but it is not limited to this.Data source is in the space flight Ke Gong Second Academy, flight time On March 22nd, 2015, flight range are Shaanxi Province Yan Liangshi.Preferably, the present invention is by taking two scape images in air strips as an example.It is led It is image, as shown in table 1 from video imaging parameter.It is well known that in spaceborne/airborne SAR image Data processing, in order to obtain height The data processed result of precision is needed using control point data purification SAR treatment process, therefore often through laying trihedral angle angle Reflector points come stablized, be reliable, easy to measure control point coordinates, in subsequent statement use " corner reflector point " table Up to " control point ", the corner reflector point coordinate format provided is indicated in the form of the most commonly used longitude, latitude and elevation, is mapping The most commonly used coordinate expression-form of remote sensing fields.In aircraft flight, at the same laid corner reflector and using GPS into Measurement of coordinates is gone, GPS measures that coordinate is as shown in table 2, and corner reflector point coordinate includes: the name that first row indicates corner reflector point Claim, secondary series indicates that longitude, third column indicate that latitude and the 4th column indicate elevation information.

The master and slave video imaging parameter of table 1

The geodetic coordinates of 2 corner reflector point of table

Step S2 constructs the tight imaging model of SAR using the initial data of acquisition, and the SAR based on building is tightly imaged Model calculates the initial pixel coordinate (x of all corner reflector points₀,y₀), obtain the initial pixel coordinate of all corner reflector points Gather { P_init}；Below with the main image data instance for acquisition.

Step S2 relates generally to two parts, and first part is that the geographical coordinate of corner reflector point is transformed into the earth's core space is straight Angular coordinate, second part are the initial pixel coordinates for obtaining corner reflector point, and the expression of this coordinate is retouched under radar fix system It states, describes image space coordinate, the initial pixel coordinate of corner reflector point is the basis for carrying out subsequent coordinates and refining.Initially The calculation accuracy of pixel coordinate require be not it is very high, in general, corner reflector positioning accuracy can achieve dozens of pixel ?.Oblique distance-Doppler equation is generally used in solution process.It describes oblique distance in imaging process, Doppler, the earth The geometrical relationship of ellipsoid is the tight imaging model of SAR.

Specifically, step S2 specifically includes following sub-step:

The geodetic coordinates of the corner reflector point is transformed into the earth's core rectangular space coordinate by step S2.1；

The geodetic coordinates of corner reflector point is the ground point geographical coordinate of field measurement personnel actual measurement, it is transformed into the earth's core Rectangular space coordinate is because in subsequent processing, and many processing steps are completed under the rectangular coordinate system in space of the earth's core, because This needs to carry out coordinate conversion.

Earth coordinates are the coordinate informations that a point is described using longitude and latitude height, and wherein geodetic longitude L was ground point The angle of ellipsoid meridian plane and Greenwich Observatory meridian plane, geodetic latitude B be ground point ellipsoid normal and ellipsoid it is red The angle in road face, geodetic height H are distance of the ground point along ellipsoid normal to ground ellipsoid.And in the tight imaging model of SAR, It is generally adopted by the earth's core rectangular coordinate system in space, coordinate origin is overlapped with earth centroid, and Z axis is directed toward the earth arctic, and X-axis refers to To the average angle point from roofing and terrestrial equator of Greenwich, Y-axis constitutes right-handed coordinate system perpendicular to XOZ plane.It can be seen that two sit The definition for marking system is entirely different, generally requires geodetic coordinates being converted to the earth's core rectangular space coordinate, so as in same coordinate frame It is calculated under frame.

Specifically, the geodetic coordinates of corner reflector point is converted to the earth's core rectangular space coordinate according to following formula (1):

Wherein, (X_P,Y_P,Z_P) be corner reflector point the earth's core rectangular space coordinate, N is radius of curvature in prime vertical, (B, L, H) be respectively corner reflector point geodetic coordinates latitude, longitude and geodetic height, e be reference ellipsoid the first eccentricity.

Illustrate that the geodetic coordinates by corner reflector point is converted to the earth's core space right-angle seat by taking corner reflector point D002 as an example Target process.Be by the geodetic coordinates (B, L, H) of the step S1 D002 point obtained (110.0559,35.10238, 667.631), corresponding ellipsoid be WGS-84 ellipsoid, semimajor axis of ellipsoid a be 6378137m, the first eccentric ratio e be 6.6944 × 10-3, corresponding radius of curvature in prime vertical N can be expressed as N=a (1-e²sin²B)^-1/2=6385208.2m.Obtain D002 point Geodetic coordinates (B, L, H) and every ellipsoidal parameter after, substituting into above-mentioned Formula of Coordinate System Transformation can obtain:

Step S2.2 calculates the corresponding satellite imagery moment t of main image；Specifically, calculating satellite according to following formula (2) Imaging moment t:

T=t₀+Δt·y (2)

Wherein, t₀For the starting imaging time (that is, imaging time of the first row image) of main image, Δ t is imaging time Interval, y are orientation coordinate of the corner reflector point in radar fix system.

It should be noted that radar fix system is using radar flight path direction as orientation coordinate y, with radar scanning side Two-dimensional imaging is carried out to coordinate x to for distance.

Likewise, illustrating the calculating process at satellite imagery moment by taking D002 as an example.By the starting obtained by step S1 at As time (12079.8304) and imaging time interval (0.00779883) substitution above-mentioned formula (2), can obtain:

T=12079.8304+0.00779883y.

Step S2.3, the second order polynomial being fitted between satellite position vectors and the satellite imagery moment t, the second order Shown in multinomial such as formula (3):

Wherein, (X_S,Y_S,Z_S) it is satellite position vectors to be fitted,It is multinomial that second order is fitted for position Formula coefficient matrix.It, in a short time being capable of perfect fitting satellite with second order polynomial since satellite generally flies along regular tracks Position.

Specifically, coefficient matrix of polynomial in order to obtain, needs to obtain at least three groups of observed quantities.By taking main image as an example, To in the image obtained by step S1 altogether include 13 groups of observed quantities.Each group of observation time can be originated by state vector Time (12079.8304s) and corresponding time interval (2.193425s), are calculated: t_i=12079.8304+ (2.193425*i i=0,1,2..., 12)

By taking X-direction coordinate as an example, 13 equations can be constructed using 13 state position vectors that step S1 is obtained, i.e.,

X=(a₀ a₁ a₂)T

Wherein, (- 1790516.1-1790657.7-1790799.2-1790940.7-1791082.3-X= 1791223.8 -1791365.3 -1791506.9 -1791648.4 -1791789.9 -1791931.4 -1792073 - 1792214.5) ', T=(t₀ t₁ … t₁₂)′.It resolves and can get multinomial coefficient

(a_x0 a_x1 a_x2)=(- 992845.37-67.5383 0.00012459).

It repeats the above steps respectively for Y and Z, the relational matrix between satellite position and imaging moment can be constructed:

Step S2.4, the second order polynomial being fitted between satellite velocity vector and the satellite imagery moment t, the second order Shown in multinomial such as formula (4):

Wherein,For satellite velocity vector to be fitted, t is the satellite imagery moment, For velocity fitting second order polynomial coefficient matrix.

Specifically, coefficient matrix of polynomial in order to obtain, needs to obtain at least three groups of observed quantities.By taking main image as an example, To in the image obtained by step S1 altogether include 13 groups of observed quantities.Each group of observation time can be originated by state vector Time (12079.8304s) and corresponding time interval (2.193425s), are calculated: t_i=12079.8304+ (2.193425*i i=0,1,2..., 12).

With V_XFor the speed of direction, 13 equations can be constructed using 13 state position vectors that step S1 is obtained, i.e.,

V_X=(b₀ b₁ b₂)T

Wherein, V_X=(- 64.552448-64.551767-64.579169-64.578488-64.577807- 64.549033 -64.548349 -64.547665 -64.546981 -64.684219 -64.655433 -64.654748 - 64.544246) ', T=(t₀ t₁ … t₁₂)′.It resolves and can get multinomial coefficient:

(b_x0 b_x1 b_x2)=(- 5693.1802 0.93345-0.0000387).

For V_YAnd V_ZIt repeats the above steps respectively, the relational matrix between satellite velocities and imaging moment can be constructed:

Step S2.5 constructs the tight imaging model of SAR, shown in the tight imaging model of SAR such as formula (5):

Wherein, Vx_S·(X_S-X_P)+Vy_S·(Y_S-Y_P)+Vz_S·(Z_S-Z_P)=- f_d·λ·(R₀+ Δ Rx)/2 to be how general Equation is strangled,For range equation, (X_S,Y_S,Z_S) it is satellite position Vector is set, (X_P,Y_P,Z_P) be corner reflector point the earth's core rectangular space coordinate,For satellite velocity vector, x is The distance of corner reflector point is to coordinate, R₀For perigee oblique distance, Δ R is range resolution, f_dFor doppler centroid, λ is Radar wavelength.

For main image shown in the table 1 obtained by step S1, perigee oblique distance R₀(4885.9532m), distance To resolution ax R (0.5353m), doppler centroid f_d(560.9245) and radar wavelength (λ=c/f=3 × 10⁸/(3 ×10⁵)=0.008565m) it is known quantity.

Step S2.6, the tight imaging model of SAR constructed using iterative method process of solution S2.5 solve corner reflector point Initial pixel coordinate.

Step S2.6 specifically includes following sub-step:

Step S2.6.1, the initial orientation of the Doppler equation iterative calculation corner reflector point based on formula (5) is to coordinate y；The step S2.6.1 includes following sub-step:

S2.6.1.1, the pixel coordinate initial value of set angle reflector points, it is preferred that the pixel coordinate initial value set is SAR Central point (the x of main image_center,y_center)；

According to an embodiment of the invention, in the main image parameter list that table 1 provides, line number 3376, columns 5400, Therefore central point (the x of main image_center,y_center)=(2700,1688).

S2.6.1.2, according to the orientation y of corner reflector point_centerCalculate satellite imagery time t_y, and by polynomial fitting Calculate corresponding satellite sensor position vector (X_S,Y_S,Z_S) and velocity vector

Wherein, satellite imagery time t is calculated according to above-mentioned steps S2.2_y, according to the satellite of step S2.3 and S2.4 acquisition Position and speed vector multinomial, details are not described herein.

According to an embodiment of the invention, orientation coordinate y_center=1688, it is t that the satellite imagery time, which is calculated,_y= 12092.9948s the position vector of satellite sensor is (X_S,Y_S,Z_S)=(- 1791365.573,4907370.081, 3652941.435), the velocity vector of satellite sensor is

S2.6.1.3 utilizes satellite sensor position vector (X_S,Y_S,Z_S), velocity vectorAnd corner reflection Device point position vector (X_P,Y_P,Z_P), doppler values f is calculated according to Doppler equation_d, while Doppler's variation is calculated Rate f_d', in conjunction with Doppler measurement f_dop, the time variable Δ t for needing to adjust is calculated according to the following formula:

Δ t=- (f_d-f_dop)/f_d'

According to an embodiment of the invention, the time variable Δ t=-0.1403s for needing to adjust is calculated.

Wherein, doppler changing rate f_d' be to be found out by the method approximation of numerical differentiation, it is assumed that calculate present orientation to The corresponding time t of coordinate y_yF_d', if dt=0.01s, t is calculated separately using Doppler equation_yAnd t_yThe doppler values at the place+dt f_d(t_y) and f_d(t_y+ dt), then

f_d'=(f_d(t_i+dt)-f_d(t_i))/dt

S2.6.1.4, using time variable Δ t to satellite imagery time t_yIt modifies, obtains new imaging moment t_y+1 =t_y+Δt；

According to an embodiment of the invention, obtaining new imaging moment

t_y+1=t_y+ Δ t=12092.9948-0.1403=12092.8544s.

S2.6.1.5 utilizes t_y+1The iterative calculation for carrying out step S2.6.1.2-S2.6.1.4, until Δ t absolute value is all small In threshold epsilon, iteration terminates, and exports main image imaging moment t_y；Preferred embodiment in accordance with the present invention, threshold epsilon 10^-7。

According to an embodiment of the invention, main image imaging moment t is calculated_y=12092.8535s.

S2.6.1.6, by the main image imaging moment t resolved_y, the orientation coordinate of corner reflector point is calculated according to the following formula Y:y=(t_y-t₀)/Δt_y, wherein t₀For initial time, Δ t_yFor orientation time interval (as shown in table 1).

According to an embodiment of the invention, calculating the orientation coordinate y=1669.8915 of corner reflector point.

Step S2.6.2 calculates the initial distance of corner reflector point to coordinate x according to the range equation of formula (5).

By the orientation coordinate y of corner reflector point, satellite sensor position is calculated, is sat in conjunction with known corner reflector point geography Mark calculates satellite sensor to the oblique distance value of corner reflector point, and then calculates the distances of ground angle reflector points to coordinate x.

The step S2.6.2 includes following sub-step:

S2.6.2.1, the initial orientation of the corner reflector point obtained from step S2.6.1 is to coordinate y, according to following formula meters Calculate satellite sensor position vector (X_S,Y_S,Z_S):

Wherein, t is current calculating main image corresponding satellite imagery moment, t₀For main image starting imaging time (that is, The imaging time of the first row image), Δ t is imaging time interval, and y is that orientation of the corner reflector point in radar fix system is sat Mark, (X_S,Y_S,Z_S) it is satellite position vectors to be fitted,Second order polynomial coefficient matrix is fitted for position.

According to an embodiment of the invention, the orientation coordinate y=1669.89 of known corner reflector point, is calculated satellite Sensor position vector (X_S,Y_S,Z_S)=(- 1791356.456,4907373.417,3652941.425).

S2.6.2.2, by range equationIn conjunction with angle Position vector (the X of reflector points_P,Y_P,Z_P), perigee oblique distance R₀, range resolution Δ R, calculate corner reflection point it is initial away from Descriscent coordinate x.

Preferably, by the position vector (X of corner reflector point_P,Y_P,Z_P)=(- 1791657.3,4907631.1, 3647548.8), perigee oblique distance R₀=4885.9532m, range resolution Δ R=0.5353m, are calculated corner reflection point Initial distance to coordinate x=977.72736.

As described above, finally being acquired according to the orientation coordinate y and distance for calculating the corner reflector point obtained to coordinate x The initial pixel coordinate (x, y) of corner reflector point.By taking corner reflector point D002 as an example, the initial pixel coordinate result that resolves For (977.72736,1669.8915).

By step 2.2-step 2.6 it is found that in the tight imaging model of SAR of building, satellite position vectors (X_S,Y_S,Z_S) (13 discrete point satellite position informations), satellite velocity vector(13 discrete point satellite velocities information), satellite Imaging time section (initial time t₀, interlude t₁, terminate time t₂)t₀-t₂, corner reflector point geodetic coordinates (B, L, H), Ellipsoidal parameter (such as the first eccentric ratio e, radius of curvature in prime vertical N), radar wavelength λ, perigee oblique distance R₀, range resolution Δ R, doppler centroid f_dEtc. parameters or intermediate auxiliary variable pass through step S1 and obtain, design parameter is as shown in table 1,2.? Formula is related to satellite position vectors and satellite velocity vector in (5), they are can be calculated according to formula (2) and (3), Y is orientation coordinate of the corner reflector in radar fix system in formula (2).So to be solved is unknown in formula (5) Only there are two (x, y) for parameter, due to being the initial picpointed coordinate obtained for the first time, so (x, y) that is calculated is corner reflection The initial coordinate of device, as shown in table 3, with (x₀,y₀) indicate the initial pixel of corner reflector point being calculated according to formula (5) Coordinate.

In addition, the geodetic coordinates (B, L, H) of the corner reflector point of step S1 acquisition, is fieldwork geographical coordinate, not The initial pixel coordinate of corner reflector point, so being converted to rectangular space coordinate (X by coordinate_P,Y_P,Z_P), participate in step S2.5 Solution process.Using two equations of step S2.5, the two unknown parameters (x, y) can be resolved.

Step S2.7 repeats step S2.1-S2.6, calculates the initial pixel coordinate of all corner reflector points, owned Set { the P of the initial pixel coordinate of corner reflector point_init}.The initial pixel coordinate of calculated each corner reflector point is shown in Table 3。

The initial pixel coordinate of 3 corner reflector point of table

Step S3, the initial pixel coordinate set { P based on the step S2 corner reflector point obtained_init, it is shaken using maximum Width algorithm resolves Pixel-level coordinate of all corner reflector points in main image, and obtains all corner reflector points in main image Pixel-level coordinate set { (P_pix)}；

For the precision of step S2 generally in dozens of pixel, this is far from enough for the positioning of accurate corner reflector 's.The characteristic using SAR image itself is generally required at this time, by coordinate precision control within several pixels.In view of angle is anti- Emitter belongs to highlighted point target on image, therefore image amplitude can be used and carry out Pixel-level positioning.In position fixing process, need Interpolation is carried out using Fast Fourier Transform, therefore all window sizes used are 2 index.In solution process, still By taking D002 as an example, other corner reflector Pixel-level calculating coordinate methods are consistent with D002.

Step S3 specifically includes following sub-step:

Step S3.1, with the initial pixel coordinate (x of the corner reflector point obtained by step S2₀,y₀) centered on, it obtains Window size is W₁×W₁Complex matrix；

By taking D002 as an example, with centered on the initial pixel coordinate (977.72736,1669.8915) that is obtained by step S2, Acquisition window size is W₁×W₁Complex matrix.Preferably, W₁=128.

Step S3.2 calculates the W₁×W₁Complex matrix amplitude, obtain the first amplitude value matrix.

Specifically, calculating the W according to following formula₁×W₁Complex matrix amplitude:

Wherein, Re (i, j) indicates W₁×W₁Complex matrix pixel value (i, j) real part, Im (i, j) indicate W₁×W₁ Complex matrix pixel value (i, j) imaginary part.

Preferably, with W₁For=128,128 × 128 plural image is obtained, to each pixel value in complex matrix (i, j) carries out amplitude calculating, and calculated amplitude value formula isFinally obtain 128 × 128 amplitude value matrix.

Step S3.3 finds out the peak swing value in the first amplitude value matrix, and calculates the peak swing value institute Putting the pixel coordinate in main image

It is and surplus with first pixel g (1,1) for peak swing value point for example, in 128 × 128 amplitude value matrix Remaining all amplitude points are compared g (i, j), wherein 128, j=2 i=2 ... ... 128, if g (i, j) is greater than g (1,1), then G (i, j) amplitude is assigned to g (1,1), continues to participate in subsequent calculating, until all the points traversal terminates, exports peak swing value Point.

To the first amplitude matrix obtained by step S3.2, pixel coordinate of the maximum of points in raw video is calculatedAmplitude maximum point provides the coordinate of further refining of corner reflector.If do not considered thick The influence of difference, the corner reflector positioning accuracy that amplitude maximum point provides can reach a pixel.In order to reach sub-pixel-level essence Degree, is generally required and is further refined to coordinate, i.e., obtained sub-pixel-level coordinate by the way of interpolation.

Step S3.4, with the pixel coordinateCentered on, acquisition window size is W2 × W₂Complex matrix.

For example, centered on the D002 coordinate (1001,1638) of step S3.3 acquisition, acquisition window size is W₂× W₂Complex matrix.Because there is deviation in the result that too small window size will lead to interpolation algorithm, the window chosen at this time Mouth size should be less than the window size of step S3.3 selection, that is, W₂<W₁.Preferably, W₂=32.

Step S3.5 calculates the W₂×W₂Complex matrix amplitude, obtain the second amplitude value matrix；

Specifically, calculating the W according to following formula₂×W₂Complex matrix amplitude:

Wherein, Re (i, j) indicates W₂×W₂Complex matrix pixel value (i, j) real part, Im (i, j) indicate W₂×W₂ Complex matrix pixel value (i, j) imaginary part.

Preferably, with W₂For=32, obtain 32 × 32 plural image, to each pixel value in complex matrix (i, J) amplitude calculating is carried out, calculated amplitude value formula isFinally obtain 32 × 32 vibration Amplitude matrix.

Step S3.6 carries out N to the amplitude in second amplitude matrix₁Times over-sampling, obtains inverse-transform matrix, and The amplitude for calculating each pixel in the inverse-transform matrix obtains third amplitude matrix；

In order to guarantee computing speed and calculation accuracy, over-sampling is generally carried out using Fast Fourier Transform.It is basic Principle is that the mutability information (such as texture, noise, spot etc.) in image is mostly high-frequency information, after Fourier transform, These high-frequency informations are in image edge, and low-frequency information is in image center.It is converted anyway in fourier space, as long as guaranteeing high frequency And low-frequency information is not lost, signal characteristic can remain.Step S3.6 specifically includes following sub-step:

Step S3.6.1 carries out Fast Fourier Transform (FFT) to the amplitude in second amplitude matrix, goes to plural sky Between；

Step S3.6.2 carries out N to the amplitude in the second amplitude matrix for going to complex number space₁Times over-sampling；

The complex number space is divided into four parts, four parts of the upper left corner and bottom right angular coordinate are respectively as follows: A [(1,1), (W₂/ 2,W₂/2)],B[(W₂/2+1,1),(W₂,W₂/2)],C[(1,W₂/2+1),(W₂/2,W₂)],D[(W₂/2+1,W₂/2+1),(W₂, W₂)], while create an empty matrix M, for store the frequency-region signal through Fast Fourier Transform (FFT), finally to pass through assignment it Matrix M afterwards carries out inverse fast Fourier transform, obtains the complex matrix after inverse transformation, this process is to the second amplitude Amplitude in matrix carries out N₁Times over-sampling；

Step S3.6.2 specifically includes following sub-step:

The complex number space is divided into four parts by S3.6.2.1, and four parts of the upper left corner and bottom right angular coordinate are respectively as follows: A [(1,1),(W₂/2,W₂/2)],B[(W₂/2+1,1),(W₂,W₂/2)],C[(1,W₂/2+1),(W₂/2,W₂)],D[(W₂/2+1, W₂/2+1),(W₂,W₂)]；

Preferably, by taking D002 corner reflector point as an example, four groups of coordinates are respectively [(1,1), (16,16)] A, B [(17, 1),(32,16)],C[(1,17),(16,32)],D[(17,17),(32,32)]。

S3.6.2., 2, empty matrix M is created, is used to store the frequency-region signal through Fast Fourier Transform (FFT), wherein described The size of empty matrix M is (W₂×N₁)×(W₂×N₁), wherein N₁For over-sampling multiple, by the upper left of the empty matrix, upper right, a left side Under, (the W of bottom right four corners₂/2)×(W₂/ 2) pixel is assigned a value of A, B, C, D respectively, and the matrix M after assignment participates in next step Inverse Fourier transform.

Empty matrix is the frequency-region signal being used to store by Fast Fourier Transform (FFT).Its cardinal principle is, frequency domain with During spatial domain is mutually converted, the dimension of matrix can't change.Therefore, in order to guarantee the matrix after over-sampling Size is met the requirements, and the size of empty matrix M must be consistent with the matrix size after over-sampling.

Preferably, the block of pixels that the upper left of the empty matrix M, upper right, lower-left, each angle in the quadrangle of bottom right occupy is 16 × 16, in addition, over-sampling multiple N₁=4, then (W₂×N₁)×(W₂×N₁)=128 × 128.A, the size of B, C, D block is 16 × 16, contain 256 data values, by taking A block as an example, since A matrix-block is complex data, so there are [(3.14,4.48)_Ii=1, (5.01,1.70)_Ii=2,(0.43,-1.94)_Ii=3,…,(-4.51,1.70)_Ii=255,(0.02,0.09)_Ii=256], wherein ii= 1 to 256.

S3.6.2.3 carries out inverse fast Fourier transform to by the empty matrix M after assignment, obtains inverse-transform matrix, Wherein, the inverse-transform matrix is still complex matrix form；

The treatment process of above three sub-step as carries out N to the amplitude in the second amplitude matrix₁Times over-sampling.

Step S3.6.3 calculates the amplitude of each pixel in the inverse-transform matrix, obtains third amplitude matrix.

Specifically, calculating the amplitude of the inverse-transform matrix according to following formula:

Wherein, Re (i, j) indicates that the real part of the pixel value (i, j) of the inverse-transform matrix, Im (i, j) indicate the inversion Change the imaginary part of the pixel value (i, j) of matrix.

Step S3.7 finds out the peak swing value in the third amplitude value matrix, and calculates the peak swing value institute Putting the Pixel-level coordinate (x in main image₁,y₁)。

Still by taking corner reflector point D002 as an example, the pixel coordinate (x that is calculated through the above steps₁,y₁)= (1000.5,1638.0625)。

Step S3.8 repeats step S3.1-S3.7 to all corner reflector points, obtains all corner reflector points in main image In Pixel-level coordinate set { (P_pix), it is specific as shown in table 4.

4 corner reflector point Pixel-level image coordinate of table

Step S4, Pixel-level the coordinate set { (P of diagonal reflector points_pix) in all corner reflector points Pixel-level sit Mark carries out interferometric parameter measurement verification, obtains the sub-pixel-level coordinate of corner reflector point；

The SAR that this step is constructed according to step S2 tight imaging model is simultaneously handled, three in conjunction with interferometric phase equation A equation group (two equation groups and interferometric phase equation of the tight imaging model of SAR) can be under the rectangular coordinate system in space of the earth's core The three-dimensional coordinate of ground target point, referred to as landform Three-dimension Reconstruction Model are described.And the point of the corner reflector obtained in step S2.1 For coordinate equally under the rectangular coordinate system in space of the earth's core, the angle that can be obtained accordingly with landform Three-dimension Reconstruction Model evaluation procedure S3 is anti- The Pixel-level image coordinate positioning accuracy of emitter point obtains the sub-pixel-level coordinate (x of final corner reflector point₂,y₂)。

Under normal conditions, influence of the landform Three-dimension Reconstruction Model vulnerable to interference system parameter error needs to utilize high-precision Corner reflector point (i.e. control point) carry out interferometry calibration, obtain the higher interference system parameter of precision.However, corner reflection Influence of the device point vulnerable to SAR speckle noise, there is a certain error for the image pixel coordinate of acquisition, can use landform Three-dimensional Gravity Established model carries out the evaluation of corner reflector point pixel precision and rejects the biggish point of error.Based on above-mentioned analysis, step S4 is first The basic relational expression of landform Three-dimension Reconstruction Model is constructed, and Three-dimension Reconstruction Model is derived to interference based on landform Three-dimension Reconstruction Model The partial derivative of system parameter constructs interferometry calibration model accordingly；Then the Pixel-level image obtained by step S3 is sat Mark set { P_pixLeast-squares iteration calculating is carried out to interferometry calibration model, optimize and export interference system parameter, then The precision before and after each corner reflector point adjustment is counted, the biggish point of error is rejected, is iterated calculating again, until all the points Meet sets requirement, the sub-pixel-level positioning result of output angle reflector points.

Step S4 specifically includes following sub-step:

Step S4.1 is based on the tight imaging pattern of SAR and interferometric phase equation, determines three based on visual direction amount Orthogonal Decomposition Tie up reconstruction model；

In this step, it is first determined the basic relational expression of Three-dimension Reconstruction Model, the basic pass of the Three-dimension Reconstruction Model Be formula be such as formula (1) shown in including range equation, Doppler equation and interferometric phase equation equation group, then, according to Range equation, Doppler equation and interferometric phase equation, derive the Three-dimension Reconstruction Model based on visual direction amount Orthogonal Decomposition.

Wherein, Vx_S·(X_S-X_P)+Vy_S·(Y_S-Y_P)+Vz_S·(Z_S-Z_P)=- f_d·λ·(R₀+ Δ Rx)/2 to be how general Equation is strangled,For range equation,For Interferometric phase equation, (X_S,Y_S,Z_S) it is satellite position vectors, (X_P,Y_P,Z_P) be corner reflector point the earth's core rectangular space coordinate,For satellite velocity vector, x is the distance of corner reflector point to coordinate, R₀For perigee oblique distance, Δ R is distance To resolution ratio, f_dFor doppler centroid, λ is radar wavelength,For absolute interferometric phase, φ is winding interferometric phase, and k is Complete cycle issue, Q are interference pattern,For the distance of satellite primary antenna to target point,For satellite slave antenna to target point away from From.

Range equation, Doppler equation and interferometric phase equation, the equation group that they are constituted is nonlinearity, directly It is very difficult for being solved by numerical solution to them, therefore is just needed through analysis InSAR geometrical model (InSAR geometry Include three equations: range equation, Doppler equation and interferometric phase equation in model) their (three equations) are passed through three Dimension reconstruction model expression-form is solved.And visual direction amount Orthogonal Decomposition method avoids the solution of the Nonlinear Numerical between them Relationship has definite conception, solves the features such as easy.Based on this, the present invention uses visual direction amount Orthogonal Decomposition rule, passes through visual direction The form of Orthogonal Decomposition and the coordinate conversion of amount obtains the three-dimensional coordinate analytic solutions of target point.

As shown in Fig. 2, the present invention is based on the tight imaging pattern of SAR and interferometric phase equation, it is identified to be based on visual direction amount The Three-dimension Reconstruction Model of Orthogonal Decomposition are as follows:

Wherein, P is coordinate of ground point, S₁For satellite position vectors, r₁For visual direction amount,Indicate the corresponding satellite of main image Position to ground target distance (wherein,That is,For visual direction amount r₁Mould it is long),For visual direction amount r₁After normalizing Obtained unit visual direction amount.

Wherein, coordinate of ground point P is by the corresponding satellite position vectors S of main image₁With visual direction amount r₁The triangle constituted it And expression, the corresponding satellite position vectors S of main image need to be calculated₁, the distance r of the corresponding satellite of main image to ground target_s1, with And unit visual direction amount

It is well known that Three-dimension Reconstruction Model can describe corner reflector point image picpointed coordinate and ground point geographical coordinate it Between mathematical relationship, however as required for Three-dimension Reconstruction Model solution parameter there are error, need to model by Three-dimensional Gravity Type determines error term, and constructs error model (that is, interferometry calibration model) by Three-dimension Reconstruction Model, is missed with this to eliminate Poor item provides accurate interference system parameter information and participates in subsequent calculating.In the following step, step S4.1 is determined first Each variable of Three-dimension Reconstruction Model derived, obtain its analytic solutions expression-form (specifically including step S4.2-S4.4)； Secondly, the parameter for constituting the analytic solutions expression-form of Three-dimension Reconstruction Model certainly exists error, it is necessary to be modeled by Three-dimensional Gravity Type determines main error source, and partial derivative of the Three-dimension Reconstruction Model function to main error source, structure are derived on tri- directions X, Y, Z Interferometry calibration model is built, is optimized interference system parameter (specifically including step S4.5-S4.8)；Finally modeled using Three-dimensional Gravity Type and interferometry calibration model, optimize interference system parameter, excluding gross error, repeat step S4.5-S4.8, until institute There is corner reflector point to be all satisfied threshold requirement, the sub-pixel-level coordinate of output angle reflector points (specifically includes step S4.9- S4.12).Each step is described in detail below.

Step S4.2 calculates the corresponding satellite position vectors S of main image₁, calculation formula isIts In, (X_S,Y_S,Z_S) it is satellite position vectors to be fitted,Second order polynomial coefficient square is fitted for position Battle array, t are the main image corresponding satellite imagery moment.

Preferably, the satellite position vectors after fitting are calculated according to step S2.3

Step S4.3, the distance of the corresponding satellite position of calculating main image to ground targetWith unit visual direction amount

According to visual direction amount proper orthogonal decomposition, in VNW moving coordinate system, by unit visual direction amountBe decomposed into V axis, N axis and Component on W axis, is embodied as:Wherein, r_vFor unit visual direction amountComponent on V axis, r_nFor unit visual direction AmountComponent on N axis, r_wFor unit visual direction amountComponent on W axis.

That is, step S4.3 calculating isWith unit visual direction amount

The basic thought of visual direction amount proper orthogonal decomposition be using antenna phase center as under the VNW moving coordinate system of origin, Orthogonal Decomposition is carried out to unit visual direction amount.

The definition of VNW moving coordinate system: using phase center where primary antenna as coordinate origin, V axis is defined as carrier aircraft speed Direction vector, W axis are defined as carrier aircraft velocity vector and interfere the cross product of baseline vector, and N axis is jointly true by the cross product of V axis and W axis It is fixed, constitute right-handed coordinate system.In interference SAR Data processing, moving coordinate system is generally used for solving unit visual direction amount.In general, Three orthogonal basis are named as v, n and w, successively indicate V axis, N axis and W axis, referred to as VNW system.

Three orthogonal basis under VNW moving coordinate systemIts expression formula are as follows:

In above formula, b is interference baseline vector, and v is carrier aircraft flight velocity vector.

Step S4.3 specifically includes following sub-step:

Step S4.3.1 is based on unit visual direction amountUnit of account visual direction amountComponent r on V axis_v, table Up to formula are as follows:

Wherein, v is satellite flight velocity vector, | v | long for speed mould, r is visual direction amount,Mould for visual direction amount r is long.

Step S4.3.2 is based on unit visual direction amountUnit of account visual direction amountComponent r on N axis_n, table Up to formula are as follows:

Wherein,Indicate primary antenna to target point distance, | b_pv| for the baseline vector mould perpendicular to satellite velocities direction It is long, b_vProjection value for baseline vector in satellite velocities direction, r_vFor unit visual direction amountComponent on V axis, B are baseline mould Long, λ is radar wavelength,For absolute interferometric phase, Q is antenna operating mode, wherein when the operating mode of antenna is attached most importance to rerail When road mode, Q=2, when the operating mode of antenna is a hair pair receipts mode, Q=1.Preferably, in the present embodiment, Q=1.

Step S4.3.3, according to the pixel coordinate { P of corner reflector point_pixObtain corner reflector point absolute interferometric phaseWherein φ₁To wind phase, φ₂For level land phase.Wind phase pacifically phase by corresponding data file It provides, is read according to the pixel coordinate information of the step S3 corner reflector obtained.Obtain the absolute interference of corner reflector point The range difference that phase is used to calculate major and minor antenna phase center to ground target pointIn step It participates in calculating in the formula (4) of S4.3.2.

Step S4.3.4 describes unit visual direction amount in step S4.3Due toFor unit visual direction amount, Expression formula meets right-hand rule, unit visual direction amountComponent r on W axis_wIt indicates are as follows:

It is (including main image initial time, more to the imaging parameters obtained by step S1 by taking D002 corner reflector point as an example General Le multinomial, nearly oblique distance, radar wavelength, 13 satellite position vectors and velocity vector, azimuth resolution, distance are to resolution Rate, baseline vector, semimajor axis of ellipsoid, semiminor axis of ellipsoid), the step S2.3 fitting satellite position vectors and satellite imagery moment Second order polynomial, step S2.4 fitting satellite velocity vector between t and the second order between the satellite imagery moment t are multinomial Formula is computed

It is defined according to VNW coordinate system it is found that when W axis is downwardly directed, takes positive sign；When W axis points up, negative sign is taken, here It points up or lower section, and non-critical direction zenith or the earth's core direction.Preferably, positive sign is taken in the present embodiment.

Step S4.3.5 has been described that the tight imaging model of SAR is in the earth's core rectangular coordinate system in space in step S2.2 Lower completion, the geographical coordinate of the target point P provided is also to indicate under the rectangular coordinate system in space of the earth's core, and the earth's core is empty in order to obtain Between target point P under rectangular coordinate system unit visual direction amount, need the corner reflector point coordinate conversion under VNW moving coordinate system To the earth's core rectangular coordinate system in space, transition matrixAre as follows:

Wherein, b is interference baseline vector, and v is satellite flight velocity vector,For the unit vector of satellite flight speed.

Preferably, the imaging parameters (including baseline vector, satellite velocity vector), the step S2 that are obtained by step S1 are obtained The satellite velocity vector obtained, has been computed:

Step S4.4, the corresponding satellite position vectors S of main image obtained based on step S4.2₁, step S4.3 obtains WithDerive the Three-dimension Reconstruction Model are as follows:

Wherein, P is coordinate of ground point, S₁For satellite position vectors, r₁For visual direction amount, r_s1Indicate the corresponding satellite of main image To the distance of ground target point, For transition matrix,For unit visual direction amount, r_vFor unit visual direction amount Component on V axis, r_nFor unit visual direction amountComponent on N axis, r_wFor unit visual direction amountComponent on W axis, b are dry Baseline vector is related to, v is carrier aircraft flight velocity vector,For the unit vector of flying speed, λ is radar wavelength,It is absolutely dry Phase is related to, Q is antenna operating mode, | b_pv|, b long for the baseline vector mould perpendicular to satellite velocities direction_vExist for baseline vector The projection value in satellite velocities direction, f_dopFor doppler centroid, B is that interference baseline mould is long.

Illustrate the process calculated, the corresponding satellite position vectors of main image by taking corner reflector point D002 as an exampleOblique distance r_s1=5421.521, spin matrixAnd Unit visual direction amountHave:

Step S4.2-S4.4 has carried out the derivation of equation and description to the basic representation that step S4.1 is provided, it is determined that three The Analytical Expression form for tieing up reconstruction model, provides some theoretical reference foundations for subsequent step S4.5-S4.8, will be according to step The Three-dimension Reconstruction Model that S4.2-S4.4 is determined carries out the calculating of partial derivative according to the Three-dimension Reconstruction Model to interference system parameter, The building of interferometry calibration model (error model) and least square solve.Step S4.5- is detailed below The treatment process of S4.8.

Step S4.5 rewrites the Three-dimension Reconstruction Model of step S4.4, obtains the observational equation of Three-dimension Reconstruction Model；

During the Three-dimension Reconstruction Model that S4.1-S4.4 is established, every system parameter is inevitably by error It influences, in order to obtain accurate every system parameter, needs to carry out calibration processing, optimization interference system ginseng using corner reflector point Number, exports high-precision interference system parameter, provides reliable interference system ginseng finally to obtain high-precision corner reflector coordinate Number.Interferometry calibration model derives the partial derivative of each system parameter based on Three-dimension Reconstruction Model, and then utilizes minimum two Multiply principle solving system parameter, optimization and the pixel coordinate for improving corner reflector point, and excluding gross error.Specifically, to formula (7) It is rewritten, obtains its function representation form, the i.e. observational equation of Three-dimension Reconstruction Model:

F in formula_x(X)、F_y(X)、F_zIt (X) is respectively the function of X, Y, Z-direction,For to excellent The interference system parameter of change, b_lFor baseline length, b_αFor baseline angle, r₀For nearly oblique distance,For interferometric phase, f_dFor Doppler's frequency Rate, these parameters both participate in the calculating of Three-dimension Reconstruction Model, and there are certain error, need to carry out these parameters excellent Change processing；S_1x,S_1y,S_1zFor the corresponding coordinate of the satellite position of main image,For main image to ground corner reflector (i.e. appearance Punctuate) oblique distance, u₁₁,u₁₂,u₁₃,u₂₁,u₂₂,u₂₃,u₃₁,u₃₂,u₃₃For spin matrixCalculating variable；r_vFor unit visual direction AmountComponent on V axis, r_nFor unit visual direction amountComponent on N axis, r_wFor unit visual direction amountComponent on W axis； P_x,P_y,P_zFor the earth's core rectangular space coordinate of corner reflector, the major and minor image parameters information acquisition provided by step S1.

Step S4.6 calculates the partial derivative of the interference system parameter in the observational equation of the Three-dimension Reconstruction Model；

According to the formula (8) of step S4.5, formula (8) is to interference system parameter (including b_lFor baseline length, b_αFor baseline Inclination angle, r₀For nearly oblique distance,For interferometric phase, f_dFor Doppler frequency) derivation linearized to obtain error equation (that is, interference Measure calibration model), it is therefore desirable to the partial derivative for calculating interference system parameter specifically includes following sub-step:

Step S4.6.1 calculates separately corner reflector point F on tri- directions X, Y, Z according to formula (8)_x、F_yAnd F_zFunction To baseline length b_lPartial derivative.We are with F_xTo baseline length b_lIt asks explanation for partial derivative to seek the process of partial derivative, is computed HaveWhereinFor the oblique of main image to ground corner reflector (i.e. ground target point) Away from (distance), spin matrixUnit visual direction amount r_vFor unit visual direction amountComponent on V axis, r_nFor unit visual direction amountComponent on N axis, r_wFor unit visual direction amountIn W Component on axis, finds through analysis:And u₁₁,u₁₂,u₁₃Calculating process be not related to baseline length b_l, however's Calculating needs to use baseline length b_l, soAnd u₁₁,u₁₂,u₁₃To baseline length b_lPartial derivative be 0, unit visual direction amountTo baseline length b_lPartial derivative be not 0, be calculative, last F_xTo baseline length b_lPartial derivative is asked just to convert For unit visual direction amountTo baseline length b_lSeek partial derivative.F_yAnd F_zTo baseline length b_lAsk partial derivative and F_xTo baseline length b_lIt asks Partial derivative is similar.Therefore, the F on tri- directions X, Y, Z_x、F_yAnd F_zFunction is to baseline length b_lAsking partial derivative conversion is that unit regards VectorTo baseline length b_lIt differentiates, obtains:

Wherein, b_lFor baseline length；Unit visual direction amountr_vFor unit visual direction amountComponent on V axis, r_nFor Unit visual direction amountComponent on N axis, r_wFor unit visual direction amountComponent on W axis；For main image to ground corner reflection The oblique distance of device point, Δ r are the distance satellite corresponding with main image of the corresponding satellite of main image to corner reflector point to corner reflection The difference of the distance of device point.

Step S4.6.2 calculates separately F on tri- directions X, Y, Z according to formula (8)_x、F_yAnd F_zFunction is to baseline angle b_α Partial derivative, the unit visual direction amount in moving coordinate systemTo baseline angle b_αDerivative be all 0, baseline angle b_αOnly establishing It is had an impact when moving coordinate system, to b_αDerivative are as follows:

Wherein, b_αFor baseline angle, unit visual direction amountr_vFor unit visual direction amountComponent on V axis, r_nFor Unit visual direction amountComponent on N axis, r_wFor unit visual direction amountComponent on W axis；,For main image to ground corner reflection The oblique distance of device,For spin matrix.

In formula (11)To b_αPartial derivative is sought, is had:

Wherein,WithFor three orthogonal basis under VNW moving coordinate system, b_αFor baseline angle,For spin moment Battle array.

Step S4.6.3 calculates separately F on tri- directions X, Y, Z according to formula (8)_x、F_yAnd F_zFunction is to initial oblique distance Partial derivative.Oblique distance of the primary antenna to ground angle reflector pointsWith initial oblique distance r₀Between relationship:Along away from The image coordinate of descriscent Shangdi millet cake is n_rg, the pixel resolution in oblique distance is r_ps.Formula (9) is to r₀It differentiates, obtains:

Wherein, F_x、F_yAnd F_zThe function of respectively X, Y, Z-direction, unit visual direction amountr_vFor unit visual direction amount Component on V axis, r_nFor unit visual direction amountComponent on N axis, r_wFor unit visual direction amountComponent on W axis, u₁₁, u₁₂,u₁₃,u₂₁,u₂₂,u₂₃,u₃₁,u₃₂,u₃₃For spin matrixCalculating variable, r₀For nearly oblique distance,For primary antenna to ground The oblique distance of corner reflector point.

Unit visual direction amount is sought againTo r₀Derivative:

In above formula, the baseline vector component vertical with satellite velocities directionFormula (14) are substituted into formula (13) In seek r₀Partial derivative.Wherein,For primary antenna to the oblique distance of ground angle reflector points, unit visual direction amountr_vFor list Position visual direction amountComponent on V axis, r_nFor unit visual direction amountComponent on N axis, r_wFor unit visual direction amountOn W axis Component, b_pvFor the baseline vector component vertical with satellite velocities direction, Δ r is main subpictures satellite position to corner reflector point Apart from its difference, b_lFor baseline length.

Step S4.6.4 calculates separately F on tri- directions X, Y, Z according to formula (8)_x、F_yAnd F_zFunction is to interferometric phase Partial derivative, obtain:

Wherein, unit visual direction amountIt is rightDerivative are as follows:

Wherein,For interferometric phase,For primary antenna to the oblique distance of ground angle reflector points, unit visual direction amount r_vFor unit visual direction amountComponent on V axis, r_nFor unit visual direction amountComponent on N axis, r_wFor unit visual direction amountIn W Component on axis, b_pvFor the baseline vector component vertical with satellite velocities direction,For absolute interferometric phase, b_lFor base length Degree；Q is antenna operating mode, wherein when the operating mode of antenna is repeat track mode, Q=2, when the Working mould of antenna When formula is a hair pair receipts mode, Q=1；b_pvFor the baseline vector component vertical with satellite velocities direction, λ is radar wavelength.

Step S4.6.5 calculates separately F on tri- directions X, Y, Z according to formula (8)_x、F_yAnd F_zFunction is to Doppler center Frequency f_dPartial derivative, obtain:

Wherein, unit visual direction amountTo doppler centroid f_dDerivative are as follows:

In formula, f_dFor doppler centroid, b_vFor the interference baseline projection components parallel with satellite velocities direction；Unit Visual direction amountr_vFor unit visual direction amountComponent on V axis, r_nFor unit visual direction amountComponent on N axis, r_wFor Unit visual direction amountComponent on W axis, b_pvFor the baseline vector component vertical with satellite velocities direction, λ is radar wavelength, | v | it is long for the mould of satellite velocity vector.

Step S4.7, the partial derivative of the interference system parameter based on the obtained observational equation of step S4.5 and S4.6 derivation, Nonlinear least squares fitting is carried out to interference system parameter using several corner reflector points；

Interferometry calibration model is constructed using corner reflector point.Assuming that in SAR image, there are m_GA corner reflector point, There is error equation (that is, interferometry calibration model) at each corner reflector point:

In formula, v_xG(i),v_yG(i),v_zG(i)Respectively corner reflector point G_x(i),G_y(i),G_z(i) X, Y, Z-direction it is residual Difference,Respectively F_xG(i),F_yG(i),F_zG(i)Function is to baseline length b_lIn X, Y, the local derviation of Z-direction Number；Respectively F_xG(i),F_yG(i),F_zG(i)Function is to baseline angle b_αIn X, Y, the local derviation of Z-direction Number；Respectively F_xG(i),F_yG(i),F_zG(i)Function is to nearly oblique distance r₀In X, Y, the partial derivative of Z-direction；Respectively F_xG(i),F_yG(i),F_zG(i)Function is to interferometric phaseIn X, Y, the partial derivative of Z-direction；Respectively F_xG(i),F_yG(i),F_zG(i)Function is to Doppler frequency f_dopIn X, Y, the local derviation of Z-direction Number；l_xG(i),l_yG(i),l_zG(i)For the difference structure for the geocentric coordinates that geocentric coordinates and Three-dimension Reconstruction Model known to corner reflection point calculate At constant term, corner reflector point G (i) indicate i-th of corner reflector point, i=1,2 ..., n_G, add up to n_GA corner reflector Point.

Step S4.8 resolves interference system parameter based on the principle of least square, specifically includes following sub-step:

Step S4.8.1 lists error equation according to error model listed by formula (18)

V=A Δ x-l (19)

Wherein,

Matrix V is residual error item matrix, wherein v_xG(1), v_yG(1), v_zG(1)Indicate a corner reflector point of G (1) in X, Y, the side Z To residual error,Indicate G (m_G) a corner reflector point is in X, Y, the residual error of Z-direction；Matrix A is Partial derivative coefficient matrix, whereinRespectively F_xG(i),F_yG(i),F_zG(i)Function is to baseline length b_l? X, Y, Z-direction partial derivative；Respectively F_xG(i),F_yG(i),F_zG(i)Function is to baseline angle b_αX, Y, the partial derivative of Z-direction；Respectively F_xG(i),F_yG(i),F_zG(i)Function is to nearly oblique distance r₀In X, Y, Z The partial derivative in direction；Respectively F_xG(i),F_yG(i),F_zG(i)Function is to interferometric phaseIn X, Y, the side Z To partial derivative；Respectively F_xG(i),F_yG(i),F_zG(i)Function is to Doppler frequency f_dopIn X, Y, Z The partial derivative in direction, i=1,2 ..., n_G；Matrix l is constant term matrix l_xG(i),l_yG(i),l_zG(i)Respectively L_xG(i)-L_0,xG(i), L_yG(i)-L_0,yG(i), L_zG(i)-L_0,zG(i)It is calculated, i=1,2 ..., n_GFor i-th of corner reflector point of calculating, n_GFor corner reflection Device point total number；Δ x is that interference system parameter to be solved corrects moment matrix, baseline length reduction Δ b_l, baseline angle changes Positive quantity Δ b_α, nearly oblique distance reduction Δ r₀, interferometric phase reductionWith Doppler frequency reduction Δ f_d。

Step S4.8.2, the video imaging parameter (shown in table 1) obtained by step S1, design factor matrix A, in step Its calculation formula has been described in S4.7 and S4.8, brings calculating into；The corner reflection point the earth's core space obtained by step S2 is straight Angular coordinate, computational constant item l=L₀- L, L are calculated by Three-dimension Reconstruction Model formula (7), L₀Angle is surveyed by field measurement personnel Reflector points position obtains.

Design factor matrix A, calculation formula are shown below:

We calculate separately D002 point F on tri- directions X, Y, Z by taking corner reflector point D002 as an example_x、F_yAnd F_zFunction To baseline length b_lPartial derivative, illustrate specific calculating process.Three-dimensional Gravity modeling is described in the formula (8) of step S4.5 The function representation form of type, is shown below:

Calculate D002 point F on tri- directions X, Y, Z_x、F_yAnd F_zFunction is to baseline length b_lPartial derivative, haveBe related to main image to ground corner reflector oblique distanceSpin matrixUnit visual direction amountTo baseline length b_lPartially fall Number

It calculatesAccording to step S2.5 it is found thatAccording to provided in table 4 D002 point pixel coordinate (x, Y)=(1000.5,1638.0625)；Nearly oblique distance R in the table 1 that step S is provided₀=4885.9532m, distance is to sampling interval Δ R =0.536m；It finally obtains

Calculating spin matrix can according to step S4.3.5 Know, in conjunction with velocity vector and baseline vector that step S1 is provided, obtains

Unit of account visual direction amount is to baseline length b_lCalculating partially reciprocal, the number that can be obtained according to step S4.3.1, S4.3.2 and S4.3.3 Value and b_l=0.312816m, It can obtain

For baseline angle b_α, nearly oblique distance r₀, interferometric phaseDoppler frequency f_dopIt repeats the above steps respectively Obtain coefficient matrices A.D002 point F on tri- directions X, Y, Z is listed below_x、F_yAnd F_zFunction is to baseline length b_l, baseline inclines Angle b_α, nearly oblique distance r₀, interferometric phaseDoppler frequency f_dopPartial derivative, it is as follows:

Computational constant item matrix l, calculation formula are shown below:

We have by taking corner reflector point D002 as an exampleWherein L_x,L_y,L_zBy the three-dimensional of step S4.4 Reconstruction model formula (7) is calculated；L_0,x,L_0,y,L_0,zThe longitude that is there is provided by table 2, latitude, elevation i.e. (110.0559, 35.10238,667.631) geodetic coordinates at the control point of corner reflector, is converted to the earth's core space right-angle through step S2.1 to sit Mark；It finally obtains

l_x,D002=L_0,x-L_x=-1791657.022- (- 1791633.082)=- 23.940

l_y,D002=L_0,y-L_y=4907631.293-4907598.887=32.406.

l_z,D002=L_0,z-L_z=3647548.734-3647534.327=14.406

Step S4.8.3, adjustment criterion select common criterion of least squares: V^TV=min；Solve interference system parameter Correction has:

Δ x=(A^TA)^-1A^Tl (20)

Wherein, Δ x is that interference system parameter to be solved corrects moment matrix, and A is the coefficient that step S4.8.2 is calculated Matrix, A^TFor the transposed matrix of matrix A, l is constant term matrix.

Step S4.8.4 calculates the root mean square of the coordinate difference of each iteration front-and-back angle reflector points, judge whether accordingly into The iterative calculation of row next step, formula are as follows:

In formula, eps is the root mean square of the coordinate difference of each iteration front-and-back angle reflector points,Indicate i-th of corner reflection The coordinate that kth time after the error model optimization that device point is constructed according to step S4.8 iterates to calculate out,Indicate i-th jiao it is anti- Kth -1 time coordinate iterated to calculate out after the error model optimization that emitter point is constructed according to step S4.8, n_GFor corner reflector Point total number, when eps is less than the threshold value of setting, adjustment terminates, and counts the residual error before and after each corner reflector point adjustment.If adjustment In the process, abnormal end is then needed to check relevant parameter, be recalculated, until meeting the requirements.

Preferably, eps takes 10^-4。

In step S4.5-S4.8, the calculating of partial derivative is carried out to interference system parameter according to the Three-dimension Reconstruction Model, is done Relate to measurement calibration model (error model) building and least square solve and etc., be fully described pushing away for error model It leads, construct, solving, exporting the treatment processes such as result, providing corner reflector point coordinate for subsequent step S4.9-S4.12.And In step S4.9-S4.12, the pixel precision of diagonal reflector points evaluate simultaneously excluding gross error, exports final corner reflection Device point sub-pixel-level coordinate.

Step S4.9, the interference system parameter obtained using step S4.8, to the corner reflector picpointed coordinate in step S3 Interferometric parameter measurement verification is carried out, the biggish corner reflector point of error is rejected；

By step S4.8 it is found that the Pixel-level image coordinate { P that step S3 is obtained_pix, since SAR image is made an uproar vulnerable to spot Certainly there is the biggish point of error in the influence of sound, these points will affect the precision of interferometry calibration, accordingly can be according to interference It measuring calibration model and carries out Detection of Gross Errors and rejecting, the automatic rough error point for detecting corner reflector iterates to calculate, excluding gross error point, It realizes sub-pixel-level automatic positioning, specifically includes following sub-step:

The Pixel-level image coordinate set { P that step S4.9.1, step S3 are obtained_pixIn all corner reflector picpointed coordinates Interferometry calibration compensating computation is both participated in, until meeting least square constraint condition, iteration terminates, at this time will be anti-by model The interference system parameter drilled calculates the corresponding the earth's core 3 d space coordinate of corner reflector point, knot according to the formula (7) of step S4.4 Fruit is denoted as { P_{XYZ_m}}。

The corresponding the earth's core 3 d space coordinate fundamental formular of corner reflector point is calculated to be shown below:

So according to interference system parameter (including the baseline length b after optimization_l, baseline angle b_α, nearly oblique distance r₀, interference phase PositionWith Doppler frequency f_d) and major and minor image parameters list (including baseline vector, 13 satellite positions arrow that provides of step S1 Amount and velocity vector, radar wavelength, nearly oblique distance) and step S3.8 provide 4 corner reflector point pixel coordinate of table.One by one by angle Reflector points bring formula calculating into, are as a result denoted as { P_{XYZ_m}}。

The earth's core rectangular space coordinate collection known to the corner reflection point obtained by step S2.1 is combined by step S4.9.2 {P_XYZ, calculate two set { P_XYZAnd { P_{XYZ_m}In the plane mean square error of a point { P of X-direction and Y-direction_{XYZ_err}, and obtain flat The mean value and standard deviation μ, σ of millet cake position error of coordinate, there is following calculation formula:

P_{XY_err}=SQRT ((P_{XYZ_m}(X)-P_XYZ(X))²+(P_{XYZ_m}(Y)-P_XYZ(Y))²)

Wherein, P_{XYZ_err}For the plane mean square error of a point of X-direction and Y-direction, SQRT () is to open radical sign operation, P_{XYZ_m}(X) The earth's core 3 d space coordinate { the P obtained for step S4.9.1_{XYZ_m}X-direction, P_XYZIt (X) is the corner reflector point of field operation actual measurement The earth's core 3 d space coordinate { P_XYZX-direction, P_{XYZ_m}(Y) the earth's core 3 d space coordinate obtained for step S4.9.1 {P_{XYZ_m}Y-direction, P_XYZIt (Y) is the earth's core 3 d space coordinate { P of the corner reflector point of field operation actual measurement_XYZY-direction.μ is The mean value of plane point error of coordinate, P_{XY_err}It (i) is { P_{XYZ_err}Element in set, n_GFor corner reflector point total number.σ is The standard deviation of plane point error of coordinate.

Step S4.9.3, according to the distance of the step S1 image obtained to, azimuth resolution be (x_Res,y_Res), it can be with Obtain corresponding picture point ground resolutionIf the standard deviation sigma of the plane point of step S4.9.2 is greater than Picture point ground resolution (xy)_Res, then there are the biggish corner reflector points of error for explanation.

Preferably, (xy)_Res=0.75m.

Step S4.9.4, the standard deviation sigma and picture point ground resolution k* (xy) of Calculation Plane point_ResDifference, it is contemplated that angle The positioning accuracy of reflector points will meet the requirement of sub-pixel-level, it is necessary to be better than a pixel or the corresponding ground of a pixel Unit-sized is differentiated, k indicates the proportionality coefficient k ∈ (0.1,1) for being better than a ground resolution cell.According to 1 times of standard deviation sigma threshold value Excluding gross error, shown in following formula, if σ > k* (xy)_Res, then retain the point that numberical range in { e } meets following condition, be discontented with The point of sufficient condition is automatically set as Null Spot, is no longer participate in compensating computation, until all corner reflector points for participating in adjustment meet The threshold condition of setting；

e∈(μ-σ,μ+σ)

Preferably, k=0.5, k* (xy)_Res=0.375m.

Step S4.9.5 carries out interferometry calibration adjustment and elimination of rough difference operation by step S4.7-S4.9.4, passes through Operation is terminated three times, meets sub-pixel-level positioning accuracy request, is obtained model optimization parameter (that is, interference system parameter), system It is as shown in table 5 below to count result:

5 statistical model operation result of table

Step S4.10, the model optimization parameter obtained using step S4.9.5 (including baseline length b_l, baseline angle b_α、 Nearly oblique distance r₀, interferometric phaseAnd Doppler frequency), angle is recalculated according to the tight imaging model of SAR that abovementioned steps S2 is constructed The image picpointed coordinate of reflector points, the sub-pixel-level corner reflector point coordinate finally obtained are (x₂,y₂), and by all coordinates Set is denoted as { P_{pix_f}It is shown in Table 6.

Corner reflector point image coordinate after the optimization of 6 geometrical relationship of table

The above is only a preferred embodiment of the present invention, for those of ordinary skill in the art, according to the present invention Thought, there will be changes in the specific implementation manner and application range, and the content of the present specification should not be construed as to the present invention Limitation.

Claims (10)

1. a kind of sub-pixel-level corner reflector automatic positioning method based on the tight imaging model of SAR, which is characterized in that the side Method specifically includes the following steps:

S1, acquires initial data, the initial data include the coordinate of corner reflector point, main image and the main image at As parameter；

S2 constructs the tight imaging model of SAR using the initial data of acquisition, and the tight imaging model of the SAR based on building calculates Initial pixel coordinate (the x of all corner reflector points₀,y₀), obtain the set of the initial pixel coordinate of all corner reflector points {P_init}；

S3, the initial pixel coordinate set { P based on the step S2 corner reflector point obtained_init, it is resolved using peak swing algorithm Pixel-level coordinate x1, y1 of all corner reflector points in main image, and obtain picture of all corner reflector points in main image Plain grade coordinate set { (P_pix)}；

S4, Pixel-level the coordinate set { (P of diagonal reflector points_pix) in the Pixel-level coordinates of all corner reflector points done Parameter measurement verification is related to, sub-pixel-level the coordinate x2, y2 of corner reflector point are obtained.

2. the method according to claim 1, wherein step S2 specifically includes following sub-step:

The geodetic coordinates of the corner reflector point is transformed into the earth's core rectangular space coordinate by S2.1；Specifically, according to following formula (1) geodetic coordinates of corner reflector point is converted into the earth's core rectangular space coordinate:

Wherein, (X_P,Y_P,Z_P) be corner reflector point the earth's core rectangular space coordinate, N is radius of curvature in prime vertical, (B, L, H) point Not Wei corner reflector point geodetic coordinates latitude, longitude and geodetic height, e be reference ellipsoid the first eccentricity；

S2.2 calculates the corresponding satellite imagery moment t of main image；Specifically, calculating the satellite imagery moment according to following formula (2) T:

T=t₀+Δt·y (2)

Wherein, t₀For the starting imaging time of main image, that is, the imaging time of the first row image, Δ t are imaging time interval, y The orientation coordinate for being corner reflector in radar fix system；

S2.3, the second order polynomial being fitted between satellite position vectors and the satellite imagery moment t, the second order polynomial is such as Shown in formula (3):

Wherein, (X_S,Y_S,Z_S) it is satellite position vectors to be fitted,Second order polynomial system is fitted for position Matrix number；

S2.4, the second order polynomial being fitted between satellite velocity vector and the satellite imagery moment t, the second order polynomial is such as Shown in formula (4):

Wherein,For satellite velocity vector to be fitted, t is the satellite imagery moment,For speed It is fitted second order polynomial coefficient matrix；

S2.5 constructs the tight imaging model of SAR, shown in the tight imaging model of SAR such as formula (5):

Wherein, Vx_S·(X_S-X_P)+Vy_S·(Y_S-Y_P)+Vz_S·(Z_S-Z_P)=- f_d·λ·(R₀+ Δ Rx)/2 be Doppler side Journey,For range equation, (X_S,Y_S,Z_S) it is that satellite position is sweared Amount, (X_P,Y_P,Z_P) be corner reflector point the earth's core rectangular space coordinate,For satellite velocity vector, x is that angle is anti- The distance of emitter point is to coordinate, R₀For perigee oblique distance, Δ R is range resolution, f_dFor doppler centroid, λ is radar Wavelength；

S2.6, the tight imaging model of SAR constructed using iterative method process of solution S2.5 solve the initial pixel of corner reflector point Coordinate；

S2.7 repeats step S2.1-S2.6, calculates the initial pixel coordinate of all corner reflector points, obtain all corner reflectors Set { the P of the initial pixel coordinate of point_init}。

3. according to the method described in claim 2, it is characterized in that, step S2.6 specifically includes following sub-step:

S2.6.1, the initial orientation of the Doppler equation iterative calculation corner reflector point based on formula (5) is to coordinate y；

S2.6.2 calculates the initial distance of corner reflector point to coordinate x according to the range equation of formula (5).

4. according to the method described in claim 3, it is characterized in that, the step S2.6.1 includes following sub-step:

S2.6.1.1, the pixel coordinate initial value of set angle reflector points, wherein the pixel coordinate initial value set is SAR main image Central point (x_center,y_center)；

S2.6.1.2, according to the orientation y of corner reflector point_centerCalculate satellite imagery time t_y, and calculated by polynomial fitting Corresponding satellite sensor position vector (X_S,Y_S,Z_S) and velocity vector

S2.6.1.3 utilizes satellite sensor position vector (X_S,Y_S,Z_S), velocity vectorWith corner reflector point Set vector (X_P,Y_P,Z_P), doppler values f is calculated according to Doppler equation_d, while doppler changing rate f is calculated_d', In conjunction with Doppler measurement f_dop, the time variable Δ t for needing to adjust is calculated according to the following formula:

Δ t=- (f_d-f_dop)/f_d'；

S2.6.1.4, using time variable Δ t to satellite imagery time t_yIt modifies, obtains new imaging moment t_y+1=t_y+ Δt；

S2.6.1.5 utilizes t_y+1The iterative calculation for carrying out step S2.6.1.2-S2.6.1.4, until Δ t absolute value is both less than threshold Value ε, iteration terminate, and export main image imaging moment t_y；

S2.6.1.6, by the main image imaging moment t resolved_y, the orientation coordinate y:y=of corner reflector point is calculated according to the following formula (t_y-t₀)/Δt_y, wherein t₀For initial time, Δ t_yFor orientation time interval.

5. according to the method described in claim 3, it is characterized in that, the step S2.6.2 includes following sub-step:

S2.6.2.1, the initial orientation of the corner reflector point obtained from step S2.6.1 are defended to coordinate y according to the calculating of following formula Star sensor position vector (X_S,Y_S,Z_S):

Wherein, t is the current calculating main image corresponding satellite imagery moment, and t0 is the starting imaging time of main image, that is, first The imaging time of row image, Δ t are imaging time interval, and y is orientation coordinate of the corner reflector point in radar fix system, (X_S,Y_S,Z_S) it is satellite position vectors to be fitted,Second order polynomial coefficient matrix is fitted for position；

S2.6.2.2, by range equationIn conjunction with corner reflection Position vector (the X of device point_P,Y_P,Z_P), perigee oblique distance R₀, range resolution Δ R, calculate corner reflection point initial distance to Coordinate x.

6. the method according to claim 1, wherein step S3 specifically includes following sub-step:

S3.1, with the initial pixel coordinate (x of the corner reflector point obtained by step S2₀,y₀) centered on, obtain window size For the complex matrix of W1 × W1；

S3.2 calculates the W₁×W₁Complex matrix amplitude, obtain the first amplitude value matrix；

Specifically, calculating the W according to following formula₁×W₁Complex matrix amplitude:

Wherein, Re (i, j) indicates W₁×W₁Complex matrix pixel value (i, j) real part, Im (i, j) indicate W₁×W₁Plural number The imaginary part of the pixel value (i, j) of matrix；

S3.3 finds out the peak swing value in the first amplitude value matrix, and calculates the peak swing value point in master Pixel coordinate in image

S3.4, with the pixel coordinateCentered on, acquisition window size is W2 × W₂Complex matrix；

S3.5 calculates the W₂×W₂Complex matrix amplitude, obtain the second amplitude value matrix；

Specifically, calculating the W according to following formula₂×W₂Complex matrix amplitude:

Wherein, Re (i, j) indicates W₂×W₂Complex matrix pixel value (i, j) real part, Im (i, j) indicate W₂×W₂Plural number The imaginary part of the pixel value (i, j) of matrix；

S3.6 carries out N to the amplitude in the second amplitude value matrix₁Times over-sampling, obtains inverse-transform matrix, and calculate institute The amplitude for stating each pixel in inverse-transform matrix obtains third amplitude value matrix；

S3.7 finds out the peak swing value in the third amplitude value matrix, and calculates the peak swing value point in master Pixel-level coordinate (x in image₁,y₁)；

S3.8 repeats step S3.1-S3.7 to all corner reflector points, obtains pixel of all corner reflector points in main image Grade coordinate set { (P_pix)}。

7. according to the method described in claim 6, it is characterized in that, W2 < W1.

8. according to the method described in claim 6, it is characterized in that, step S3.6 specifically includes following sub-step:

S3.6.1 carries out Fast Fourier Transform (FFT) to the amplitude in the second amplitude value matrix, goes to complex number space；

S3.6.2 carries out N to the amplitude in the second amplitude value matrix for going to complex number space₁Times over-sampling；

S3.6.3 calculates the amplitude of each pixel in the inverse-transform matrix, obtains third amplitude value matrix；

Specifically, calculating the amplitude of the inverse-transform matrix according to following formula:

Wherein, Re (i, j) indicates that the real part of the pixel value (i, j) of the inverse-transform matrix, Im (i, j) indicate the inverse transformation square The imaginary part of the pixel value (i, j) of battle array.

9. the method according to claim 1, wherein step S4 specifically includes following sub-step:

S4.1 is based on the tight imaging pattern of SAR and interferometric phase equation, determines the Three-dimensional Gravity modeling based on visual direction amount Orthogonal Decomposition Type, wherein the Three-dimension Reconstruction Model are as follows:

Wherein, P is coordinate of ground point, S₁For satellite position vectors, r₁For visual direction amount,Indicate the corresponding satellite position of main image To the distance of ground target,For visual direction amount r₁The unit visual direction amount obtained after normalizing；

S4.2 calculates the satellite position vectors S₁, calculation formula isWherein, (X_S,Y_S,Z_S) it is satellite position vectors to be fitted,It is fitted second order polynomial coefficient matrix for position, t is The main image corresponding satellite imagery moment；

S4.3, the distance of the corresponding satellite position of calculating main image to ground targetWith unit visual direction amount

S4.4, the corresponding satellite position vectors S of main image obtained based on step S4.2₁, step S4.3 obtainsWithIt pushes away Export the Three-dimension Reconstruction Model are as follows:

Wherein, P is coordinate of ground point, S₁For satellite position vectors, r₁For visual direction amount, r_s1Indicate the corresponding satellite of main image to ground The distance of appearance punctuate, For transition matrix,For unit visual direction amount, r_vFor unit visual direction amountIn V axis On component, r_nFor unit visual direction amountComponent on N axis, r_wFor unit visual direction amountComponent on W axis, b are interference base Line vector, v are carrier aircraft flight velocity vector,For the unit vector of flying speed, λ is radar wavelength,Absolutely to interfere phase Position, Q is antenna operating mode, | b_pv|, b long for the baseline vector mould perpendicular to satellite velocities direction_vIt is baseline vector in satellite The projection value of directional velocity, f_dopFor doppler centroid, B is that interference baseline mould is long；

S4.5 rewrites the Three-dimension Reconstruction Model of step S4.4, obtains the observational equation of Three-dimension Reconstruction Model:

In formula, F_x(X)、F_y(X)、F_zIt (X) is respectively the function of X, Y, Z-direction,It is to be optimized Interference system parameter, b_lFor baseline length, b_αFor baseline angle, r₀For nearly oblique distance,For interferometric phase, f_dFor Doppler's frequency Rate；S_1x,S_1y,S_1zFor the corresponding coordinate of the satellite position of main image,For main image to the oblique distance of ground corner reflector, u₁₁,u₁₂, u₁₃,u₂₁,u₂₂,u₂₃,u₃₁,u₃₂,u₃₃For spin matrixCalculating variable；r_vFor unit visual direction amountComponent on V axis, r_n For unit visual direction amountComponent on N axis, r_wFor unit visual direction amountComponent on W axis；P_x,P_y,P_zFor the ground of corner reflector Heart rectangular space coordinate；

S4.6 calculates the partial derivative of the interference system parameter in the observational equation of the Three-dimension Reconstruction Model；

S4.7, based on the partial derivative of the obtained observational equation of step S4.5 and S4.6 the interference system parameter derived, utilization is several A corner reflector point carries out Nonlinear least squares fitting to interference system parameter；

S4.8 resolves the interference system parameter based on the principle of least square；

S4.9, the interference system parameter obtained using step S4.8, interferes the corner reflector picpointed coordinate in step S3 Parameter measurement verification, rejects the biggish corner reflector point of error, and obtain model optimization parameter；

S4.10, the model optimization parameter obtained using step S4.9.5, including baseline length b_l, baseline angle b_α, nearly oblique distance r₀、 Interferometric phaseWith Doppler frequency f_d, corner reflector point is recalculated according to the tight imaging model of SAR that abovementioned steps S2 is constructed Image picpointed coordinate, final corner reflector point coordinate be (x₂,y₂), and all coordinate sets are denoted as { P_{pix_f}}。

10. a kind of computer-readable medium, there is computer program in medium, the computer program is wanted for executing such as right Seek method described in any one of 1-9.