Summary of the invention
Purpose of the present invention: overcome the deficiencies in the prior art, provide a kind of spaceborne double antenna SAR of comprehensive Parameter of Overall Design to interfere scaler to lay method, the method makes overall design and interferes calibration closely to connect, can and interfere calibration to carry out l-G simulation test as unified integral body with overall design, SAR imaging processing, interference processing, both improved interference SAR system overall design quality, and given from overall angle again and interfere calibration that important references is provided.
Technical solution of the present invention: a kind of spaceborne double antenna SAR of comprehensive Parameter of Overall Design interferes scaler to lay method, to analyze borne SAR to the coverage condition on ground according to population parameter, the method of recycling sensitivity matrix analysis draws the ground calibration device and lays rule, and obtain the latitude and longitude information of each scaler by coordinates transformation method, can analyze for each latitude zone, the whole world, realize the interference calibration analysis for global range.
Principle of the present invention is: the ultimate principle that the spaceborne double antenna SAR of given first interferes calibration.Figure 1 shows that the principle of interference figure of spaceborne double antenna InSAR, h
InSARThe earth's core distance for the terrain object place can obtain following expression according to geometric relationship among the figure:
Wherein, f () is the height conversion relation; X is interferometric parameter, comprises satellite altitude H, oblique distance r
1, base length B, baseline inclination angle theta
bAnd interferometric phase
In the interferometry process, above-mentioned parameter can affect the precision of digital elevation model (DEM) data all with error, usually will utilize ground control point (GCP:Grond Control Point) that these parameters are calibrated, and namely interferes calibration.
Figure 2 shows that the InSAR system interferes the schematic diagram of calibration, establishes the parameter error in the data acquisition
The interferometric parameter that then obtains is expressed as
For certain impact point in the mapping band, the linearization error model of its height reconstruction error is
In the mapping band, lay L scaler, and the elevation of known each scaler is h
i, i=1,2 ..., L can obtain the poor of the reconstruction elevation of each scaling point and actual elevation, and matrix equation can be expressed as
Δ=F·ΔX+M (4)
Wherein,
Be the vertical error of i scaling point, i=1,2 ..., L; M is the linearization error matrix;
Be sensitivity matrix, the concrete form of sensitivity matrix is suc as formula shown in (5),
The interferometric parameter error can be obtained by the finding the solution of system of linear equations of formula (4), is expressed as: Δ X=F
+Δ.Wherein, F
+Generalized inverse for sensitivity matrix F.Can realize thus the interference calibration of spaceborne double antenna InSAR system.Fig. 2 has provided concrete realization flow.
In the interferometry of reality, the position of ground calibration device can't accurately be known, all with certain error, and interferometric parameter is also with certain stochastic error, therefore observation data Δ band is made an uproar, when finding the solution overdetermined equation (4), error can pass to by certain weight interferometric parameter X, the degree of propagation of error depends on the conditional number of sensitivity matrix F, and the situation that lays of the conditional number of sensitivity matrix and ground calibration device has much relations, must carry out scaler by certain constraint condition for this reason and lay, just can guarantee the calibration precision of interferometric parameter.Namely meeting under certain geometrical constraint condition, making conditional number C (F) minimum of sensitivity matrix F.C (F) computing method are suc as formula (6),
C(F)=||F||
2·||F
+||
2 (6)
Wherein || F||
2With || F
+||
2Difference representing matrix F and F
+2 norms.
The present invention introduces population parameter and interferes calibration, can carry out the analysis that scaler lays rule for each latitude area of global range, will based on the wave beam covering analyzing of satellite orbit parameter and unified based on the interference calibration of sensitivity matrix analysis in same set of coordinate system, be the key point that realizes the method.Used coordinate system should reflect the physical process of spaceborne double antenna InSAR work in overall design analysis and the digital emulation, therefore needs comparatively complete coordinate system framework, so that interfere sensitivity matrix also can represent in this coordinate frame.
Fig. 3 has provided the space geometry graph of a relation of spaceborne double antenna interference SAR, has clearly reflected the inner link between ellipsoid earth model, satellite orbit, major-minor SAR antenna and the ground scene, can provide being defined as follows of each coordinate system accordingly:
(1) earth inertial coordinates system E
O: true origin is at the earth centre of sphere, and X-axis in the plane, is pointed to the first point of Aries under the line, and Z axis points to the positive arctic along the axis of rotation of the earth, and Y-axis in the plane, meets right-hand rule relation with X, Z axis under the line.
(2) body-fixed coordinate system E
G: true origin is at the earth centre of sphere, and in the plane, by first branch of Greenwich meridian, Z axis points to the positive arctic along the axis of rotation of the earth to X-axis under the line, and Y-axis in the plane, meets right-hand rule relation with X, Z axis under the line.
(3) satellite orbit coordinate system E
V: true origin is the earth centre of sphere (being positioned at the satellite orbit plane), and X-axis is in the satellite orbit plane, and forward points to pericenter, Z axis is perpendicular to the satellite orbit plane, forward points to the angular-momentum vector direction of satellite, and Y-axis is in the satellite orbit plane, and direction is determined by the right-hand rule.
(4) satellite platform coordinate system E
R: true origin is at centroid of satellite, and X-axis is in the satellite orbit plane, and forward points to the design heading of satellite, and Z axis is perpendicular to the satellite orbit plane, and forward points to the angular-momentum vector direction of satellite, and Y-axis is in the satellite orbit plane, and direction is determined by the right-hand rule.
(5) satellite celestial body coordinate system E
E: true origin is at centroid of satellite, and X-axis is along satellite celestial body y direction, and forward points to the Live Flying direction of satellite, and Y-axis and Z axis are respectively along two other principal axis of inertia direction of satellite celestial body.
(6) main antenna coordinate system E
A: true origin is at the antenna phase center point, and the X-axis forward points to the Live Flying direction of satellite, and Y-axis is along antenna boresight, and forward points to the ground direction of bowl.
(7) satellite synchronization scene coordinate system E
SS: as shown in Figure 4, true origin is the ground aiming point of main antenna beam center, and XOY plane and earth's spheroid are tangential on the O point, the X-axis positive dirction is pointed to the satellite flight direction, the normal that Z axis is ordered at O along earth's spheroid, forward deviate from the ground direction of bowl, and Y direction is determined by the right-hand rule.
(8) earth synchronization scenarios coordinate system E
SE: as shown in Figure 4, true origin is the ground aiming point of main antenna beam center, and XOY plane and earth's spheroid are tangential on the O point, the X-axis forward points to the due east direction, the Y-axis positive dirction is pointed to direct north, and the normal that Z axis is ordered at O along earth's spheroid, forward deviate from the ground direction of bowl.As we know from the figure, earth synchronization scenarios coordinate system E
SeRotate to an angle around Z axis and can obtain satellite synchronization scene coordinate system Ess.
Performing step of the present invention is as follows:
The first step is found the solution satellite Covering time t
0:
Can set up following system of equations according to geometric relationship shown in Figure 5,
δ=arcsin[sin(ω+φ)sini] (9)
λ=Ω+arctan[tan(ω+φ)cosi]-[G
0+ω
e(t
0-t
p)] (10)
Λ
P-λ=arccos[(cosψ-sinδsinΦ
P)/cosδcosΦ
P] (11)
Wherein τ is the moment of satellite process pericenter, and a represents the major semi-axis of elliptical orbit, and e represents orbital eccentricity, ω is argument of perigee, and i is orbit inclination, and Ω represents the red footpath of ascending node, above six parameter general name satellite orbit six key elements will have overall design to provide.μ represents gravitational field gravitational constant, ω
eBe earth spin rotating speed, ω
e=7.2921158 * 10
-5Rad/s, G
0Be t
pThe sidereal hour angle of moment Greenwich, Λ
P, Φ
PRepresent respectively geocentric longitude and geocentric latitude that terrain object P is ordered, above-mentioned parameter is known quantity.δ, λ represent respectively t
0Geocentric latitude and the geocentric latitude of moment substar, ψ is substar N
SAnd the earth's core angle between the impact point P, E represent inclined to one side pericenter angle, and φ represents true pericenter angle, and these parameters are unknown quantity, can find the solution by above-mentioned system of equations.Formula (7) is for describing the Kepler's equation of satellite position and time relationship, and formula (8) represents the relation at very near heart angle and partially near heart angle, and formula (9-11) is set up according to the spherical trigonometry relation, specifically can be referring to Fig. 5.O represents the earth centre of sphere among the figure, and N represents the earth's axis arctic, and G represents Greenwich meridian and equatorial node, and D represents the perigee, S
mBe satellite position, P is t
0The ground aiming point at moment antenna beam center, N
SBe substar.
Utilize numerical computation method can find the solution above-mentioned system of equations, obtain satellite Covering time t
0
Substar N
SAnd the earth's core angle ψ computing method between the impact point P are as follows:
R wherein
eRepresent local earth radius, θ
AExpression antenna beam centre visual angle.
Second step calculates the main antenna wave beam at the coverage condition on ground, obtains the longitude and latitude of the ground aiming point of beam center, near-end wave beam line and far-end wave beam line, and calculates mapping swath width, specifically is calculated as:
If the main antenna centre visual angle is θ
A, distance is Δ θ to beam angle
Ar,
(2.1) the longitude Λ of calculating main antenna beam center aiming point
0With latitude Φ
0
At first utilize following formula Calculation of Satellite pole of orbit radius vector r:
Can calculate t thus
0Satellite is at earth inertial coordinates system E constantly
OIn position (x
Os, y
Os, z
Os), shown in (14):
A wherein
OVBe satellite orbit coordinate system E
VTo earth inertial coordinates system E
OTransition matrix.
Then set up main antenna coordinate system E
AMiddle any point (x
a, y
a, z
a) at body-fixed coordinate system E
GIn coordinate.Because the Y-axis of antenna coordinate system overlaps with antenna boresight, so the main antenna aiming point is at coordinate system E
AIn coordinate be (0, y, 0), aiming point is at body-fixed coordinate system E
GIn coordinate (x
Go, y
Go, z
Go) be:
(x wherein
e, y
e, z
e) be that antenna phase center is with respect to satellite celestial body coordinate system E
EPosition vector, (x
Os, y
Os, z
Os) for satellite at earth inertial coordinates system E
OIn position vector, A
GOEarth inertial coordinates system E
OTo body-fixed coordinate system E
GTransition matrix, A
OVSatellite orbit coordinate system E
VTo earth inertial coordinates system E
OTransition matrix, A
VRSatellite platform coordinate system E
RTo satellite orbit coordinate system E
VTransition matrix, A
RESatellite celestial body coordinate system E
ETo satellite platform coordinate system E
RSatellite platform coordinate system E
RTransition matrix, A
EATo be θ by the main antenna centre visual angle
AThe antenna coordinate of determining is tied to satellite celestial body coordinate system E
ETransition matrix.With (x
Go, y
gO, z
Go) substitution ellipsoid model of globe equation:
E wherein
aAnd E
bBe respectively major semi-axis and the minor semi-axis of earth's spheroid.Solve an equation and to try to achieve y (get less in two solutions, cast out larger solution), can obtain aiming point at body-fixed coordinate system E simultaneously
GIn coordinate (x
Go, y
Go, z
Go), aiming point longitude Λ is then arranged
0For:
Aiming point latitude Φ
0For:
(2.2) longitude and latitude (Λ of the ground aiming point of calculating main antenna near-end wave beam line
N, Φ
N), the longitude and latitude (Λ of the ground aiming point of far-end wave beam line
F, Φ
F);
The longitude of near-end wave beam line ground aiming point is:
The latitude of near-end wave beam line ground aiming point is:
Wherein
That the ground aiming point of main antenna near-end wave beam line is at body-fixed coordinate system E
GIn coordinate vector, the computation process of its computation process and main antenna beam center aiming point is similar, difference is for calculating transition matrix A
EAAntenna look angle, use main antenna centre visual angle θ in the step (2.1)
A, use the near-end view angle theta herein
AN, θ as shown in Figure 6
AN=θ
A-Δ θ
Ar/ 2.
The longitude of far-end wave beam line ground aiming point is:
The latitude of far-end wave beam line ground aiming point is:
Wherein
That the ground aiming point of main antenna far-end wave beam line is at body-fixed coordinate system E
GIn coordinate vector, the computation process of its computation process and main antenna beam center aiming point is similar, difference is for calculating transition matrix A
EAAntenna look angle, use main antenna centre visual angle θ in the step (2.1)
A, use the far-end view angle theta herein
AF, θ as shown in Figure 6
AF=θ
A+ Δ θ
Ar/ 2.
(2.3) cross the width W that the constantly calculation of longitude ﹠ latitude mapping of the ground aiming point of corresponding antenna proximal end wave beam line and far-end wave beam line of top is with according to satellite:
W=ψ
NF·R
e (23)
ψ wherein
NFCorresponding the earth's core angle between the ground aiming point for antenna proximal end wave beam line and far-end wave beam line.
Can obtain ψ according to spherical trigonometry relation among Fig. 5
NFComputing formula:
ψ
NF=arccos[sin(Φ
N)sin(Φ
F)+cos(Φ
N)cos(Φ
F)cos(Λ
F-Λ
N)] (24)
The 3rd step, utilize mapping swath width W that second step tries to achieve and satellite system parameter for InSAR system made sensitivity matrix, described satellite system parameter comprises satellite altitude H, base length B, baseline inclination angle theta
bAccording to the constraint condition that scaler is laid, to try to achieve scaler and lay parameter, the described parameter that lays comprises distance d between scaler number L, first scaler and the substar
1And the distance interval delta d between the scaler, it is as follows that scaler lays the calculation of parameter step:
(3.1) set up the sensitivity matrix F of spaceborne double antenna interference SAR:
Wherein
The ground elevation f at expression the i scaler place is to the susceptibility of interferometric parameter X, and described interferometric parameter X comprises satellite altitude H, main antenna oblique distance r
1, base length B, baseline inclination angle theta
bAnd interferometric phase
(illustrate:
In follow-up formula, be not used for calculating, only be used for
The expression of this symbol.)
Wherein
Re is local earth radius, and h is the height that scaler place terrain object arrives the earth ellipsoid face, and θ is main antenna visual angle corresponding to scaler place, and Δ r is that the major-minor antenna is poor to the oblique distance of terrain object.
(3.2) be located at the interior edge of mapping band distance to uniformly-spaced laying L scaler, change the distance d between first scaler and the substar
1And the distance interval delta d between each scaler, adopt numerical analysis method to calculate (L, d
1, Δ d), under the condition that satisfies formula (31), make conditional number C (F) minimum of sensitivity matrix F.
D wherein
NExpression mapping band near-end distance, D
FExpression mapping band distal end distance, Δ d
MinHave at least in the expression (26-30) one can not carry out at the scaling point place that linearization represents minimally apart from the interval.
The 4th goes on foot, and calculates the longitude and latitude of each scaler;
(4.1) according to the 3rd distance d that goes on foot between scaler number L, first scaler and the substar of obtaining
1And the distance interval delta d between the scaler, calculate each scaler at satellite synchronization scene coordinate system E
SSIn position vector;
If the distance of L scaler is respectively d
i, i=1 ..., L, i.e. i scaler shown in Fig. 7 and substar N
SBetween arc length, then each scaler is at satellite synchronization scene coordinate system E
SSIn position vector be V
SS=(0, d
i-d
0, 0), i=1 ..., L, wherein d
0Be the ground aiming point of main antenna beam center and the distance between the substar, i.e. P shown in Fig. 7 and N
SBetween arc length, computation process is as follows:
d
0=ψ
0·R
e (32)
ψ wherein
0Be the earth's core angle between beam center ground aiming point and the substar, can do following calculating according to the spherical trigonometry relation:
ψ
0=arccos[sin(Φ
0)sin(δ)+cos(Φ
0)cosδcos(Λ
0-λ)] (33)
(4.2) with scaler at satellite synchronization scene coordinate system E
SSIn position vector be transformed into earth synchronization scenarios coordinate system E
SEIn,
V
SE=A
ESV
SS (34)
A wherein
ESFor by satellite synchronization scene coordinate system E
SSTo earth synchronization scenarios coordinate system E
SETransition matrix.
(4.3) each scaler position coordinates is transformed into body-fixed coordinate system E
GIn, and then calculate the longitude and latitude of each scaler, process is as follows,
A
GSBe corresponding transition matrix, the longitude of i scaler and latitude are respectively suc as formula (36) and (37):
The present invention's advantage compared with prior art is:
(1) the present invention combines the overall system design parameter with the sensitivity matrix analysis, makes the overall design of spaceborne double antenna interference SAR and interferes calibration to combine and study, can the sophisticated systems performance evaluation, and the optimization system design.
(2) the present invention is based on satellite orbit parameter and SAR systematic parameter to interfering calibration to be analyzed, adopted ellipsoid earth model and perfect coordinate-system, can carry out calibration analysis for each latitude zone, the whole world, and obtain the concrete latitude and longitude information of scaler.
(3) the present invention carries out satellite beams covering and calibration analysis under perfect coordinate system framework, overall design, interfere that processing, interference are calibrated etc. can be carried out comprehensive simulating and verify.
Embodiment
The below utilizes concrete spaceborne double antenna InSAR parameter to test to verify that the spaceborne double antenna InSAR of comprehensive population parameter interferes scaler to lay the validity of method.
Table 1 has provided the orbit parameter of satellite, and table 2 has provided main SAR systematic parameter.
Table 1 satellite orbit parameter
Satellite orbit parameter |
Numerical value |
Semi-major axis of orbit a/km |
6886.856 |
Orbital eccentricity e |
0.000726 |
Orbit inclination i/ ° |
97.63898 |
Argument of pericentre ω/° |
209.6589 |
The red footpath Ω of ascending node/° |
218.8837 |
Time of pericenter passage τ/s |
0.0 |
Table 2 SAR systematic parameter
Systematic parameter |
Numerical value |
Signal wavelength lambda/m |
0.03 |
The main antenna view angle theta
A/°
|
28.4 |
Beam angle/° |
7.3 |
Base length B/m |
200.0 |
The baseline inclination angle theta
b/°
|
45.0 |
As shown in Figure 8, the present invention's step when the spaceborne double antenna InSAR of analysis interferes scaler to lay rule is as follows:
Step (1): zone to be analyzed is equator, and establishing latitude is north latitude 0.25 degree, can be 1347.26 seconds in the hope of the satellite Covering time according to formula (7)-(11).
Step (2.1): according to the satellite Covering time, orbital elements, SAR systematic parameter are tried to achieve the aiming point longitude and latitude of main antenna beam center on ground according to formula (13-18) and are (0.248088 ° ,-64.158359 °);
Step (2.2): distance is to beam angle Δ θ
Ar=7.3 °, near-end antenna look angle θ
AN=24.75 °, far-end antenna look angle θ
AF=32.05 °, calculate the intersection point longitude and latitude on antenna beam near-end and far-end and ground according to formula (19-22), be respectively:
(Λ
N,Φ
N)=(-64.366509°,0.285488°)
(Λ
F,Φ
F)=(-63.933913°,0.207757°)
Step (2.3): with the above results substitution formula (23) and (24), calculating mapping swath width W is 48.9275Km.
Step 3: the susceptibility of at first analyzing spaceborne double antenna interference SAR according to formula (26-30), provide its sensitivity curve, shown in Fig. 9 a-Fig. 9 e, select accordingly main interferometric parameter and carry out calibration analysis, can find out that by the sensitivity curve that Fig. 9 a-Fig. 9 e provides the interference survey is high very little to the satellite altitude susceptibility, main interferometric parameter is main antenna oblique distance, base length, baseline inclination angle and interferometric phase.
To different scaler numbers (here from L=2), utilize formula (25-30) to set up sensitivity matrix, be optimized calculating according to formula (31), try to achieve minimal condition corresponding to each L value and count C
L(F), L=2 ..., 30 and corresponding d
1With Δ d value.Figure 10 and Figure 11 have provided respectively the sensitivity matrix conditional number with Δ d and d
1Change curve, can draw accordingly in the situation that the scaler number is certain, scaler distributes overstepping the bounds of propriety loose along mapping band, the sensitivity matrix conditional number is less, namely near the low coverage end, each scaler is evenly distributed in the whole mapping band first scaler as far as possible as far as possible.
Choose C
L(F), L=2 ..., the minimum value C in 30
Min(F), its corresponding scaler number is L
Opt, L
OptCorresponding allocation optimum is d
1, optWith Δ d
Opt, the scaler that obtains thus optimum lays and is configured to (L
Opt, d
1, opt, Δ d
Opt).Figure 12 has provided the change curve of sensitivity matrix minimal condition number with the scaler number, and can draw the scaler number by figure is that 8 o'clock sensitivity matrix conditional numbers are minimum, so L
Opt=8.This moment d
1Be 134.474Km, Δ d is 6115.9m.Can get allocation optimum is:
L
opt=8
d
1,opt=134.474Km
Δd
opt=6115.9m
Step (4.1): the distance that can be calculated main antenna beam center ground aiming point place by population parameter and coordinate system translation operation is 158.0162Km, i.e. d
0=158.0162km, then the position vector of each scaling point in satellite synchronization scene coordinate system is as shown in table 3.
Table 3 scaling point is position vector in satellite synchronization scene coordinate system
The scaling point sequence number |
Position vector/m |
1 |
(0,-23542.10,0) |
2 |
(0,-16552.47,0) |
3 |
(0,-9562.84,0) |
4 |
(0,-2573.21,0) |
5 |
(0,4416.41,0) |
6 |
(0,11406.04,0) |
7 |
(0,18395.67,0) |
8 |
(0,25385.30,0) |
Step (4.2) and (4.3): the angle of Calculation of Satellite velocity and due east direction, and then try to achieve the transition matrix A that arrives earth synchronization scenarios coordinate system in the satellite synchronization scene coordinate system
ES, recycling formula (34-37) is calculated the latitude and longitude information of each scaler, and is as shown in table 4.
Table 4 scaler latitude and longitude information
The scaling point sequence number |
(longitude/°, latitude/°) |
1 |
(-64.366509,0.255488) |
2 |
(-64.304710,0.248669) |
3 |
(-64.242910,0.241851) |
4 |
(-64.181111,0.235032) |
5 |
(-64.119311,0.228213) |
6 |
(-64.057512,0.221394) |
7 |
(-63.995712,0.214576) |
8 |
(-63.933913,0.207757) |
In a word, the present invention combines the design of the population parameters such as satellite orbit parameter and radar system parameter with interference calibrating method based on the sensitivity matrix analysis, provided spaceborne double antenna interference SAR scaler and laid method.Carry out the satellite covering analyzing based on population parameter, give and interfere calibration that necessary data message is provided, and then set up spaceborne double antenna interference SAR sensitivity matrix, lay scheme according to interfering the high constraint condition analysis to the sensitivity matrix conditional number of survey to obtain optimum scaler.The present invention combines overall design and interference calibration, and can in overall design, take into full account and interfere calibration to the improvement of system performance, and for interfering the Project Realization of calibrating that reference is provided.By above-mentioned test, can obtain optimum scaler number, lay the interval and specifically lay the position, validity of the present invention has been described.
The non-elaborated part of the present invention belongs to techniques well known.