A calculation method for ground target access considering terrain occlusion

Abstract

The present invention discloses a method for calculating access to ground targets considering terrain occlusion, which is particularly applicable to the field of calculating access to ground targets. Aiming at the problem that traditional calculations of access to ground targets only establish the spatio-temporal relationship between satellites and ground points, and information acquisition may fail under certain specific satellite observation postures, terrain occlusion constraints are considered, and a method and implementation process for calculating access to ground targets considering terrain occlusion are designed. This method calculates the occlusion information of the target point by constructing the spatial relationship between the target point and the occluder; after the traditional target access calculation (satellite observing the target point) is completed, the azimuth and elevation angles of the ground point observing the satellite are calculated, and by comparing with the occlusion information, the information occluded by the mountain is filtered out, thereby improving the success rate of data acquisition. This method has the characteristics of high calculation result accuracy, good practicability, wide adaptability, etc.

Description

Technical Field

[0001] The present invention relates to a method for calculating access to ground targets considering terrain occlusion, and is particularly applicable to the field of calculating access to ground targets. Background Art

[0002] With the rapid development of the space industry, satellite remote sensing imaging has become an important means for humans to obtain ground scene information. Obtaining ground scene information through satellite imaging has the advantages of being unrestricted by airspace, having a wide observation range, and high timeliness of obtaining information. In order to quickly and accurately obtain the required ground scene information, it is necessary to establish a spatio-temporal relationship between the satellite and the ground point. Through spatio-temporal relationship analysis, the optimal time and attitude for the satellite to observe the ground point are calculated, that is, the satellite access calculation technology for ground targets.

[0003] Traditional satellite access calculation technology for ground targets only establishes a spatio-temporal relationship between the satellite and the ground point, that is, only calculates the optimal access time and attitude of the satellite to the ground point, without analyzing the occlusion situation around the ground point. In some specific satellite observation attitudes, the required scene information may not be obtained (due to mountain occlusion), resulting in mission failure.

[0004] To solve this problem, the present invention proposes a method for calculating access to ground targets considering terrain occlusion. In addition to considering the time and attitude of the satellite observing the ground point, this technology adds the ability to analyze the occlusion azimuth and occlusion angle of the ground point, thereby solving the problem of being unable to obtain the specified target scene information due to terrain occlusion, avoiding waste of satellite resources, and improving the mission success rate. Summary of the Invention

[0005] The problem to be solved by the present invention is to avoid the problem of information acquisition failure caused by mountain occlusion in the above background art, and provide a method for calculating access to ground targets considering terrain occlusion. In addition to considering the time and attitude of the satellite observing the ground point, this technology adds the ability to analyze the occlusion azimuth and occlusion angle of the ground point, thereby solving the problem of being unable to obtain the specified target scene information due to terrain occlusion, avoiding waste of satellite resources, and improving the mission success rate.

[0006] The technical problem to be solved by the present invention is achieved by the following technical solutions:

[0007] A method for calculating access to ground targets considering terrain occlusion includes the following steps:

[0008] (1) Obtain the geographical information of the target point and the geographical information of the surrounding occluders of the target point, and pre-calculate the surrounding occlusion angle information of the target point; wherein, the geographical information includes geographical latitude, geographical longitude, and geographical altitude, and the occlusion angle information includes occlusion azimuth angle and occlusion pitch angle;

[0009] (2) Obtain satellite ephemeris information, calculate the attitude angle of the satellite's observation target according to the space transformation formula, and filter out the time periods that meet the satellite attitude angle ability constraints according to the satellite attitude angle ability constraints, and save the corresponding satellite attitude angle information;

[0010] (3) Calculate the azimuth and elevation angle information of the satellite relative to the target point, and filter out the occluded time periods according to the pre-calculated occlusion angle information, and screen out the time periods that are not occluded by the occluder;

[0011] (4) Further filter out the visible time periods and information that are not occluded according to the satellite attitude ability constraints.

[0012] Among them, the specific calculation of the occlusion angle information in step (1) is as follows:

[0013] (101) Convert the ground point and the occluder from geodetic coordinates to earth-fixed coordinates respectively;

[0014] (102) Convert the occluder from earth-fixed coordinates to the ground point's local horizon coordinate system matrix, and obtain the position r of the occluder in the ground point's local horizon coordinate system zxdp ;

[0015] (103) Decompose r zxdp to obtain the occlusion azimuth angle a and the occlusion elevation angle b, then:

[0016]

[0017]

[0018] In the formula, arctan represents the arctangent, with the unit of radian; when x>0 and y>0, a = a; when x>0 and y<0, a = 180 - a; when x<0 and y<0, a = 270 - a; when x<0 and y>0, a = 360 - a, arcsin represents the arcsine, with the unit of radian, x represents the x component of r zxdp the y component of r zxdp the z component of r zxdp the modulus of r zxdp is len.

[0019] Among them, the calculation method of the attitude angle of the satellite's observation target in step (2) is as follows:

[0020] (201) Calculate the vector of the ground point in the orbital coordinate system

[0021] Record the position of the ground point in the J2000.0 inertial coordinate system as r1, and the position of the satellite in the J2000.0 inertial coordinate system as r0, then the position r of the ground point in the satellite orbital coordinate system orb is:

[0022] r orb = C(r1 - r0)

[0023] (202) Calculate the satellite's pointing attitude towards the target

[0024] Set the angle of the satellite rotating around the X-axis as the roll angle a, the angle of rotating around the Y-axis as the pitch angle b, the Z-axis pointing to the center of the earth, the X-axis as the satellite flight direction, O as the center of the coordinate system, and P as the target point; and when the satellite rotates, the rotation around the Z-axis is not considered, and the rotation sequence is roll first and then pitch;

[0025] Then, the calculation formula for the roll angle a is:

[0026] a = arctan(Y / Z)

[0027] That is, the arctangent of the ratio of the component Y to the component Z in the position r of the ground point in the satellite orbit coordinate system orb ;

[0028] The calculation formula for the pitch angle b is:

[0029] b = arcsin(X / |OP|)

[0030] That is, the arcsine of the ratio of the component X in the position r of the ground point in the satellite orbit coordinate system orb to the modulus value of the OP vector.

[0031] Among them, step (4) is specifically:

[0032] Set the satellite's side-sway angle range as: the maximum side-sway angle MaxCbj, the minimum side-sway angle MinCbj, set the pitch angle range as: the maximum pitch angle MaxFyj, the minimum pitch angle MinFyj, the azimuth angle of the shelter as ZbwFwj, and the shelter angle of the shelter as ZbwZbj;

[0033] Set in the time period information that meets the satellite attitude angle ability constraints, the side-sway angle is denoted as Cbj, the pitch angle is denoted as Fyj, the azimuth angle is denoted as Fwj, and the elevation angle is denoted as Yj;

[0034] Then, further screen out the visible time periods that are not blocked and the information meets: Cbj >= MinCbj, and Cbj <= MaxCbj, and Fyj >= MinFyj, and Fyj <= MaxFyj, fabs(Fwj - ZbwFwj) < 0.1, and Yj < ZbwZbj; where fabs(Fwj - ZbwFwj) < 0.1 means that the difference between Fwj and ZbwFwj does not exceed 0.1 degrees.

[0035] The present invention has the following advantages compared with the background technology:

[0036] 1. The present invention uses a standard space coordinate system conversion model and parameters, and the calculation results have high accuracy;

[0037] 2. The present invention can calculate the target occlusion information in advance and save it, so as to calculate the occlusion information in advance and directly use it as known data when performing real-time analysis of satellite access to ground targets, which has good engineering practicability;

[0038] 3. The present invention can be used in the fields of optical satellite access calculation to ground targets, SAR satellite access calculation to ground targets, electronic satellite access calculation to ground targets, etc., and has the characteristics of wide adaptability to access. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Figure 1 is the overall flowchart of the method and implementation for calculating satellite access to ground targets considering terrain occlusion of the present invention.

[0040] Figure 2 is the flowchart of calculating ground occlusion information of the present invention.

[0041] Figure 3 is the flowchart of calculating the yaw angle and pitch angle of the satellite to the ground target of the present invention.

[0042] Figure 4 is the schematic diagram of the satellite pointing to the target attitude angle of the present invention.

[0043] Figure 5 is the flowchart of calculating the azimuth and pitch of the satellite relative to the ground point of the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0044] Next, in combination with Figure 1 , Figure 2 , Figure 3 , Figure 4 , Figure 5 the present invention will be further described.

[0045] In combination with Figure 1 , a method for calculating satellite access to ground targets considering terrain occlusion includes the following steps:

[0046] (1) Obtain the geographical information of the target point (geographical latitude, geographical longitude, geographical altitude), obtain the geographical information of the surrounding occluders of the target point (geographical latitude, geographical longitude, geographical altitude), and pre-calculate the surrounding occlusion angle information of the target point (occlusion azimuth angle, occlusion pitch angle).

[0047] (2) Obtain the satellite ephemeris information (satellite position information and velocity information in the J2000 coordinate system), calculate the attitude angle of the satellite observing the target according to the space conversion formula, and screen the time periods that meet the satellite attitude angle ability constraints according to the satellite attitude angle ability constraints, and save the corresponding satellite attitude angle information.

[0048] (3) Calculate the azimuth and elevation angle information of the satellite relative to the target point, and perform occlusion period screening (screen the periods not occluded by the occluder) according to the pre-calculated occlusion angle information (the azimuth and elevation angle information of the occluder relative to the target point).

[0049] (4) Further screen out the visible periods and information that are not occluded according to the satellite attitude ability constraints (slew angle range, elevation angle range).

[0050] The above steps are further explained as follows:

[0051] Let the satellite ephemeris data at a specific moment be SATECI_0(X0, Y0, Z0, VX0, VY0, VZ0), where X0, Y0, Z0 represent the position in the J2000 coordinate system; VX0, VY0, VZ0 represent the velocity in the J2000 coordinate system.

[0052] Among them, combined with Figure 2 , the method for calculating the occlusion angle information in step (1) is as follows:

[0053] (101) Convert the geodetic coordinates of the ground point to the geocentric-fixed coordinate system

[0054] The calculation method for converting the geodetic coordinates (longitude, latitude, altitude) of the ground point to the position coordinates in the geocentric-fixed coordinate system is as follows:

[0055]

[0056] Among them, (B, L, H) are the geodetic latitude, geodetic longitude, and geodetic height; Re = 6378137m is the equatorial radius of the earth, and c is the eccentricity of the earth's meridian.

[0057] (e c ) 2 = 2f - f 2

[0058] f is the geometric flattening of the reference ellipsoid, and the calculation formula is:

[0059] f = 1 / 298.257223563

[0060] (102) Convert the geodetic coordinates of the occluder to the geocentric-fixed coordinate system

[0061] The calculation formula is referred to (101).

[0062] (103) Calculate the matrix from the geocentric-fixed coordinate system to the local horizon coordinate system

[0063] Let the position vector in the geocentric-fixed coordinate system be The position vector in the local horizon coordinate system is The position coordinates of the measuring station in the geocentric-fixed coordinate system are If the geodetic longitude and latitude of the survey station are L and B, then:

[0064]

[0065] Where:

[0066]

[0067] M is the matrix sought.

[0068] (104) The vector of the obstacle in the local horizon coordinate system of the station

[0069] The position of the obstacle in the geocentric fixed coordinate system is r zd , the position of the ground point in the geocentric fixed coordinate system is r mb , then the position r zxdp of the obstacle in the local horizon coordinate system of the ground point station is:

[0070] r zxdp = (M)[r zd - r mb

[0071] (105) Calculation of obstacle information

[0072] Decompose r zxdp to obtain the azimuth angle and elevation angle. Denote the azimuth angle as a and the elevation angle as b, then:

[0073]

[0074]

[0075] In the above formulas, arctan represents the calculation of the arctangent, with the unit of radian; when x > 0 and y > 0, a = a; when x > 0 and y < 0, a = 180 - a; when x < 0 and y < 0, a = 270 - a; when x < 0 and y > 0, a = 360 - a, arcsin represents the calculation of the arcsine, with the unit of radian, x represents the x - component of r zxdp , y represents the y - component of r zxdp , z represents the z - component of r zxdp , and len represents the modulus of r zxdp .

[0076] The calculation method of the attitude angle of the satellite observing the target in step (2) is as follows:

[0077] Combined with Figure 3 , construct the spatio - temporal relationship between the satellite and the ground point. Assume that the input conditions include the geodetic coordinates (longitude, latitude, altitude) of the ground point and the position - velocity vector of the satellite in the J2000 inertial coordinate system.

[0078] (201) Geodetic Coordinates of Ground Points Transformed into Earth-Fixed Coordinates

[0079] See formula (101) for the calculation formula.

[0080] (202) Earth-Fixed Coordinates Transformed into J2000 Inertial Coordinates

[0081] Define (PR) as the precession matrix, (NR) as the nutation matrix, (ER) as the Earth rotation matrix, and (EP) as the polar motion matrix. They are respectively expressed by the following formulas:

[0082] (EP) = R y (-x p )R x (-y p )

[0083] (ER) = R z (S G )

[0084] (NR) = R x (-Δε)R y (Δθ)R z (-Δμ)

[0085] (PR) = R z (-z A )R y (θ A )R z (-ζ A )

[0086] In each formula, x p , y p are polar motion components, S G is Greenwich sidereal time, and other quantities are precession and nutation values.

[0087] For the position vectors r J2000 and r EG in the J2000.0 inertial coordinate system and the Earth-fixed coordinate system, there is:

[0088] r J2000 = (PR) - (NR) - (ER) - (EP) - r EG

[0089] (203) Calculate the Transformation Matrix from J2000 Inertial Coordinate System to Orbital Coordinate System

[0090] Denote the transformation matrix from the J2000.0 inertial coordinate system to the satellite orbital coordinate system as C. According to the definition of the satellite orbital coordinate system, the elements of the transformation matrix C are:

[0091]

[0092]

[0093] C(1,i) = C(2,j) × C(3,k)

[0094] where: i = 1, 2, 3 correspond to the three components of each row vector in the transformation matrix C. C(1,i), C(2,j), and C(3,k) are the row vectors of the first row, second row, and third row in matrix C respectively; represents the position vector in the J2000 coordinate system; represents the velocity vector in the J2000 coordinate system.

[0095] (204) Calculate the vector of the ground point in the orbital coordinate system

[0096] If the position of the ground point in the J2000.0 inertial coordinate system is r1 and the position of the satellite in the J2000.0 inertial coordinate system is r0, then the position r of the ground point in the satellite orbital coordinate system orb is

[0097] r orb = C(r1 - r0)

[0098] (205) Calculate the attitude angles of the satellite pointing to the target

[0099] Combined with Figure 4 , it is stipulated that the angle of the satellite rotating around the X-axis is the roll angle a, and the angle of rotating around the Y-axis is the pitch angle b. Assuming the satellite rotates in the order of roll first and then pitch (without considering rotation around the Z-axis), then as Figure 3 shown, a is the roll angle, b is the pitch angle, the Z-axis points to the center of the earth, the X-axis is the flight direction of the satellite, O is the coordinate system center, and P is the target point.

[0100] The calculation formula for angle a is:

[0101] a = arctan(Y / Z)

[0102] That is, the arctangent of the ratio of the Y component to the Z component in the position r orb in the satellite orbital coordinate system of the ground point.

[0103] The calculation formula for angle b is:

[0104] b = arcsin(X / |OP|)

[0105] That is, the arcsine of the ratio of the X component to the modulus of the OP vector in the position r orb in the satellite orbital coordinate system of the ground point.

[0106] The calculation method of the azimuth and elevation angle information of the satellite relative to the target point in step (3) is as follows:

[0107] Combined with Figure 5 , assuming that the input conditions include the target point position parameters (longitude, latitude, altitude) and the satellite J2000 inertial coordinate system, the calculation process of the satellite's azimuth and elevation relative to the ground point is as follows:

[0108] (301) Convert the geodetic coordinates of the ground point to the coordinates in the earth-fixed system

[0109] The calculation formula is shown in (101).

[0110] (302) Calculate the matrix for converting the earth-fixed system to the local horizon coordinate system

[0111] The calculation formula is shown in (103).

[0112] (303) Convert the position and velocity in the J2000 inertial system to the coordinates in the earth-fixed system

[0113] For the position vectors r J2000 and r EG in the J2000.0 inertial coordinate system and the earth-fixed system, there is:

[0114] r EG =(EP)(ER)(NR)(PR)r J2000

[0115] The element descriptions are shown in (202).

[0116] (304) The vector of the satellite in the local horizon coordinate system

[0117] The calculation formula is shown in (104).

[0118] (305) Calculate the azimuth and elevation information

[0119] The calculation formula is shown in (105).

[0120] In step (4), according to the satellite attitude ability constraints (yaw angle range, pitch angle range), the unobstructed visible time periods and information are further screened out. The basis is:

[0121] Let the yaw angle range be: MaxCbj (maximum yaw angle), MinCbj (minimum yaw angle), and the pitch angle range be: MaxFyj (maximum pitch angle), MinFyj (minimum pitch angle). The azimuth angle of the obstacle is ZbwFwj, and the obstacle angle is ZbwZbj.

[0122] Let the yaw angle be denoted as Cbj, the pitch angle be denoted as Fyj, the azimuth angle be denoted as Fwj, and the elevation angle be denoted as Yj in the returned time period information.

[0123] The returned period information satisfies: Cbj >= MinCbj, and Cbj <= MaxCbj, and Fyj >= MinFyj, and Fyj <= MaxFyj, fabs(Fwj - ZbwFwj) < 0.1, and Yj < ZbwZbj. Where fabs(Fwj - ZbwFwj) < 0.1 means the difference between Fwj and ZbwFwj does not exceed 0.1 degree.

[0124] Complete the calculation of the ground target access considering terrain occlusion.

Claims

1. A method for calculating access to ground targets considering terrain occlusion, characterized in that, It includes the following steps: (1) Obtain the geographical information of the target point and the geographical information of the surrounding obstacles of the target point, and pre-calculate the surrounding occlusion angle information of the target point; wherein, the geographical information includes geographical latitude, geographical longitude and geographical altitude, and the occlusion angle information includes occlusion azimuth and occlusion elevation angle; (2) Obtain the satellite ephemeris information, calculate the attitude angle of the satellite observing the target according to the space conversion formula, and screen the time periods that meet the satellite attitude angle ability constraints according to the satellite attitude angle ability constraints, and save the corresponding satellite attitude angle information; (3) Calculate the azimuth and elevation angle information of the satellite relative to the target point, and perform occlusion time period screening according to the pre-calculated occlusion angle information to screen out the time periods that are not occluded by the obstacles; (4) Further screen out the visible time periods and information that are not occluded according to the satellite attitude ability constraints; Among them, the specific calculation of the occlusion angle information in step (1) is: (101) Convert the ground point and the obstacle from geodetic coordinates to earth-fixed system coordinates respectively; (102) Convert the geodetic coordinate system to the local-level coordinate system matrix of the ground point, and obtain the position r of the obstacle in the local-level coordinate system of the ground point zxdp ; (103) For r zxdp Perform disassembly to obtain the occlusion azimuth angle a and the occlusion elevation angle b, then: where arctan represents the calculation of the arctangent, with the unit being radians; when x > 0 and y > 0, a = a; when x > 0 and y < 0, a = 180 - a; when x < 0 and y < 0, a = 270 - a; when x < 0 and y > 0, a = 360 - a, arcsin represents the calculation of the arcsine, with the unit being radians, x represents the x-component of r zxdp , y represents the y-component of r zxdp , and z represents the z-component of r zxdp , len represents the modulus of r zxdp .

2. The method for calculating access to ground targets considering terrain occlusion according to claim 1, characterized in that, The calculation method of the attitude angle of the satellite observing the target in step (2) is: (201) Calculate the vector of the ground point in the orbital coordinate system Denote the position of the ground point in the J2000.0 inertial coordinate system as r1, and the position of the satellite in the J2000.0 inertial coordinate system as r0. Then the position r of the ground point in the satellite orbit coordinate system is: orb as follows: r orb ＝C(r1-r0) (202) Calculate the satellite pointing target attitude Set the angle of the satellite rotating around the X-axis as the roll angle a, the angle of rotating around the Y-axis as the pitch angle b, the Z-axis points to the earth's center, the X-axis is the satellite flight direction, O is the coordinate system center, and P is the target point; and the satellite rotation does not consider rotation around the Z-axis, and the rotation order is roll first and then pitch; Then, the calculation formula of the roll angle a is: a = arctan(Y / Z) That is, the position r of a ground point in the satellite orbit coordinate system orb The arctangent of the ratio of component Y to component Z in The calculation formula of the pitch angle b is: b = arcsin(X / |OP|) That is, the position r of a ground point in the satellite orbit coordinate system orb The arcsine of the ratio of the component X to the magnitude of the OP vector in 3. The method for calculating access to ground targets considering terrain occlusion according to claim 1, characterized in that, Step (4) is specifically: Set the satellite side-sway angle range as: maximum side-sway angle MaxCbj, minimum side-sway angle MinCbj, set the pitch angle range as: maximum pitch angle MaxFyj, minimum pitch angle MinFyj, the occlusion azimuth of the shelter is ZbwFwj, and the occlusion angle of the shelter is ZbwZbj; Set in the time period information that meets the satellite attitude angle ability constraints, the side-sway angle is denoted as Cbj, the pitch angle is denoted as Fyj, the azimuth angle is denoted as Fwj, and the elevation angle is denoted as Yj; Then, the visible time periods and information that are not occluded are further screened to satisfy: Cbj >= MinCbj, and Cbj <= MaxCbj, and Fyj >= MinFyj, and Fyj <= MaxFyj, fabs(Fwj - ZbwFwj) < 0.1, and Yj < ZbwZbj; where fabs(Fwj - ZbwFwj) < 0.1 means the difference between Fwj and ZbwFwj does not exceed 0.1 degree.

Images