CN112508227B - A rapid calculation method for regional target visible windows in complex conical fields of view of remote sensing satellites - Google Patents

Abstract

The invention discloses a regional target visible window rapid calculation method of a remote sensing satellite complex conical view field, wherein the calculation efficiency is improved by 10 times of magnitude in comparison with the conventional tracking propagation calculation method in terms of visibility judgment, and is improved by 10 times of magnitude in terms of global visible window calculation in comparison with the conventional method. Meanwhile, the relative error of the calculation precision compared with the STK is within 0.1 percent, and the calculation precision meets the engineering practical requirement; according to the method, satellite orbit information is obtained through Lagrange interpolation, the calculation process is independent from an orbit model, and any orbit model can be met, so that the algorithm has wider applicability; the boundary of the regional target is described as a large arc, so that the size of the regional target has no influence on the method. Compared with the traditional algorithm, the method has higher calculation efficiency and better precision under the large-area target.

Description

A rapid calculation method for regional target visible windows in complex conical fields of view of remote sensing satellites

Technical field

The invention belongs to the technical field of remote sensing satellite imaging, and specifically relates to a method for rapid calculation of regional target visible windows in complex conical fields of view of remote sensing satellites.

Background technique

Judging from existing research, the visibility issue first focused on the visibility of satellites to point targets. The most general method is called the propagation algorithm, which is to judge one by one by traversing discrete time steps. This method is accurate in calculation, but has low computational efficiency and is often used as a comparison algorithm in the literature. Follow-up research will mainly focus on improving the problem of excessive calculation amount. Lawton et al. proposed for the first time the visual function method to judge the visibility of satellite-satellite and satellite-point targets, which greatly improved the calculation efficiency. However, it is only suitable for small eccentricity orbits. AliI et al. used a large arc approximate orbit in a The sub-satellite point trajectory under periodicity has a simplified model and high calculation efficiency but is only suitable for low-orbit satellites. Mai Y et al. considered the visible window calculation of the satellite's ability to perform attitude maneuvers. Han C et al. used Hermit in the calculation of the visual function. Interpolation method, Zhang Jinxiu et al. proposed a coarse search method for satellites with variable step sizes at low latitudes and high latitudes.

However, the above studies are all methods for calculating the visible window of satellites for point targets. For regional targets, there are few existing studies, and because the type of field of view has a greater impact on the visibility of regional targets, the original methods cannot be directly promoted. Some existing research is still based on an extension of point targets, mainly sampling the interior of the region or the boundary of the region to approximate the solution, but this method will lead to a significant increase in the amount of calculation. For example, Song Zhiming et al. did not consider the satellite field of view, directly discretized the regional target boundaries into point targets and then calculated the union of their visible windows. Wang Rongfeng et al. also considered the coverage information during the calculation process. Some studies consider how to increase the coarse search range during the calculation process. E Zhibo et al. proposed using the dichotomy method for precise search to increase the initial search step size. It can be seen that there is currently no more general and efficient algorithm for the visibility calculation of regional targets, and the types of fields of view involved in the above research are only conical and rectangular fields of view. Therefore, it is very necessary to propose an accurate and efficient method suitable for complex conical fields of view.

Contents of the invention

To address this problem, the present invention proposes a semi-analytical rapid calculation method for regional target visible windows based on complex conical fields of view, which can quickly determine accurate regional target visible windows.

A method for calculating the visible window of a regional target in a complex cone field of view of a remote sensing satellite, including the following steps:

Step 1. Establish a mathematical description of the complex cone field of view: The projections of the circle corresponding to the inner half-cone angle and the circle corresponding to the outer half-cone angle of the complex cone field of view on the earth's surface are two ordinary arcs. The minimum clock angle and The projection of the two sides corresponding to the maximum clock angle on the earth's surface is two great arcs; among them, the center of the great arc coincides with the center of the sphere, and the circumference is on the earth's surface; an ordinary arc is the circumference of the circle on the earth minus the great arc. Any arc outside;

Step 2. Establish a mathematical description of the regional target: the regional target refers to a closed convex polygon on the earth's surface, and each side of the convex polygon is a great arc on the earth's sphere;

Step 3. Determine whether there is an overlapping area between the visible area A of the complex cone field of view projected on the earth's surface at each time and the area B where the regional target is located. Specifically:

(81) Determine whether area A is within area B:

(811) Traverse each edge of area B and calculate the limit point of area A relative to the edge; where the limit point refers to the point on the boundary of area A that is most likely to be outside area B; each boundary edge of area B is Corresponds to a set of limit points;

(812) For each edge of area B, determine whether its corresponding limit points are all inside the edge. If there is a limit point that is not inside the edge, exit, return "No", and enter step (82);

(813) If all the limit points on area A are inside each corresponding side of area B, then area A is within area B at that moment;

(82) Determine whether there is at least one intersection point on the boundary between area A and area B. If so, then the two areas have an intersection point; if not, exit, return "No", and enter step (83);

(83) Determine whether area B is within area A:

Determine in turn whether each boundary point of area target B is within the projection area A of the complex cone field of view. If all points are satisfied, then area B is within area A; otherwise, if one point is not satisfied, exit and return " no";

For a certain moment, as long as one of the above three conditions (81), (82), and (83) is satisfied, the field of view at that moment is visible to the regional target; if none of the above three conditions are satisfied, then the field of view at that moment is Invisible to area targets;

Step 4. Based on the positional relationship between area A and area B at each time calculated in step 3, determine the visible window of the complex cone field of view to area target B.

Preferably, in step (811), the method for calculating the limit point is:

First, calculate the Earth's half-cone angle corresponding to the inner and outer cone half-angles of the complex cone field of view. and/> Where D represents the distance of the satellite from the center of the earth, R represents the radius of the earth, θ represents the half-cone angle of the satellite's field of view,/> Indicates the desired half-cone angle of the earth’s field of view to the satellite:/>

Secondly, calculate the "limit points": take up to 6 "limit points" in one judgment, among which points 1, 2, 3, and 4 are the 4 boundary points of area A, while points 5 and 6 are based on the area To calculate the target boundary, the specific steps are as follows:

Calculate the unit vectors n _x and n _y of the unitized speed and reverse angular velocity, and obtain the unit direction vector corresponding to each "clock side" through the definition of the clock angle. Take the "clock side" where points 1 and 3 are located as an example. The expression of its unit direction vector v _1,3 is as follows:

Represents the minimum clock angle;

Using the rotation method, the position vector of the center of the circle where the complex cone field of view projection area is located is rotated to the v _{1, 3} directions respectively. angle,/> The position vector of point 1 is obtained in the same way as the position vector of point 3. Use this method to obtain the position vectors of points 2 and 4;

Calculate the 5th and 6th limit points that may exist on the boundary line of complex conical field of view projection area A relative to area B: among them, when the difference between the maximum and minimum clock half angles is less than 180 degrees, the fifth limit point may exist; when it is greater than 180 degrees When the degree is reached, there are the 5th and 6th limit points;

For the possible 5th and 6th limit points on the boundary of area A, use the normal vector of the area target boundary as its direction vector, and use the calculation method with points 1, 2, 3, 4 to calculate the position vector of the possible limit points; then according to The following formula determines whether the possible limit point is in area A:

If it satisfies the above formula, it proves that the possible limit point is within the field of view, that is, the limit point exists, otherwise the limit point does not exist;

Among them; p ₅ is the position vector of the possible limit point, the norm function represents the unitization, and the expression of the f function is as follows:

Represents the maximum clock angle; v _1/2 = f (norm (v _1,3 + v _2,4 )) represents the midpoint vector of area A; v _2,4 represents the unit of the "clock side" where point 2,4 is located direction vector.

Preferably, the rotation method is as follows:

Find the vector p' after rotating α to the unit vector v'. p' is represented by the following formula: p'=cos(α)p+sin(α)v ^* , where v ^* =norm(v'-(v'· p)v'), norm(·) represents the unitization of the vector; the unit vector v' represents the unit direction vector of a great arc or an ordinary arc.

Preferably, in step (812), the method for judging whether the corresponding limit points are all inside the side is:

Let v represent the normal vector of each side of area B, and let p represent the position vector of a certain point. If p v ≥ 0, then the vector p is on the same side as the normal vector of the side, otherwise, it is outside the normal vector of the side.

Preferably, in step (82), when determining whether area A and area B have an intersection, the following method is used to determine the maximum and small clock angle projections of the large arc where the target boundary of each area is located and the complex cone field of view. Whether the arc has an intersection point, specifically:

The intersection point Q of the above two large arcs is calculated according to the following formula Q=±norm(v×v'), and it is judged whether Q is inside the large arc d ₁ d _2. At this time, it is necessary to judge whether (d ₁ ×Q)· (Q×d ₂ ) ≥ 0 and Q · (d ₁ + d ₂ ) ≥ 0. If satisfied, it proves that Q is inside the great arc, that is, the above two great arcs have an intersection.

Preferably, in step (82), when determining whether there is an intersection between area A and area B, the following method is used to determine whether there is an intersection between the large arc at the target boundary of each area and the ordinary arc of the conic projection inside and outside the complex cone field of view. ,Specifically:

Suppose the intersection point of the ordinary arc and the great arc is Q', and let the angle between the unit vectors v and v' be β; where v and v' are respectively the unit normal vector of the plane where the great arc is located and the location of the ordinary arc. The unit normal vector of the plane;

If β+θ=90°, then the ordinary arc is tangent to the great arc; θ represents the half-cone angle of the satellite field of view corresponding to the ordinary arc;

If β+θ>90°, then the ordinary arc intersects with the great arc;

If β+θ<90°, there is no intersection between the ordinary arc and the great arc;

Among them, β is calculated as follows β=arccos(|v·v'|);

When it is determined that there is an intersection, the vector needs to be calculated: assuming that the angle between the intersection Q' and the current projection plane is γ, then cos(γ)=cos(θ)/sin(β); then the direction vector in the constructed plane is v _// =norm(v'-(v'·v)v) and the unit vector v in the vertical plane _⊥ =norm(v'×v);

Then the intersection point Q' is calculated according to the following formula: Q'=cos(γ)v _// ±sin(γ)v _⊥ ;

Finally, it is necessary to determine whether the intersection point is within the above-mentioned large arc and ordinary arc. Project Q' onto the ordinary circle plane first, and then judge whether there is an intersection point between the circles.

Preferably, the specific method of step 4 is:

First, based on the discrete step length, determine the visibility of the remote sensing satellite to the regional target at that moment; then arrange the visibility at each step length in chronological order and determine the visible start and end time ranges of each visible window; and then The precise moment of each visible window is obtained through binary search, and finally the information of each precise visible window is summarized and output.

Preferably, the satellite orbit information used in the binary search is obtained by Lagrangian interpolation.

The invention has the following beneficial effects:

For the first time, a calculation method for the visible window of regional targets in complex conical fields of view is proposed;

The computational efficiency is 10^3 times higher than the traditional tracking propagation calculation method in terms of visibility judgment, and 10^5 times faster than the traditional method in global visible window calculation. At the same time, the relative error of calculation accuracy is within 0.1% compared with STK, which meets the actual requirements of the project;

The method proposed by the present invention obtains satellite orbit information through Lagrangian interpolation. The calculation process is independent of the orbit model and can satisfy any orbit model. Therefore, the applicability of this algorithm is wider;

Since this method describes the boundary of the regional target as a large arc, the size of the regional target has no impact on this method. Compared with traditional algorithms, this method has higher computational efficiency and better accuracy under large-area targets.

Description of the drawings

Figure 1 is a schematic diagram of complex cone field of view imaging;

Figure 2 is a schematic diagram of a great circle and an ordinary circle;

Figure 3 is a schematic diagram of finding the intersection point between a large arc and an ordinary arc;

Figure 4 shows the half-cone angle of the earth’s field of view to the satellite;

Figure 5 is a schematic diagram of the meaning of "limit point";

Figure 6 is a schematic diagram for solving the "limit point" of a complex cone field of view;

Figure 7 shows the unit vector of the complex cone field of view;

Figure 8 is the overall program algorithm flow chart.

Detailed ways

The present invention will be described in detail below with reference to the accompanying drawings and examples.

The implementation steps of the present invention are as follows:

Step 1. Establish a mathematical description of the complex cone field of view

The imaging schematic diagram of a complex cone field of view is shown in Figure 1. The shape and size of the field of view are represented by four parameters: inner cone half angle, outer cone half angle, minimum clock angle, and maximum clock angle. What needs to be made clear is that under the assumption of a spherical earth, the projections of the circle corresponding to the inner half-cone angle and the outer half-cone angle of the complex cone field of view on the earth's surface are two ordinary arcs, and the minimum and maximum The projections of the two sides corresponding to the clock angle on the earth's surface are two large arcs.

Among them, the great circle mentioned in the previous paragraph refers to a circle whose center coincides with the center of the sphere and whose circumference is on the surface of the sphere. The great circle is also called a geodesic, because the distance between two points on the sphere along the geodesic is the shortest. An ordinary circle refers to any circle whose circumference is on the sphere. For the convenience of distinction, unless otherwise specified, the ordinary circles described later in the present invention do not include great circles.

Step 2. Establish a mathematical description of the regional target

The regional target refers to a closed convex polygon on the earth's surface. Each side of the convex polygon is a great arc on the earth's sphere, that is, a geodesic. By convention, the initial description method of regional target boundary points is described in the geodetic coordinate system (longitude, latitude, altitude).

Step 3. Description of the visibility problem of regional targets from the complex conical field of view of the satellite

The visible window of a satellite to a regional target refers to the time range within which the field of view of the payload it carries can observe the regional target. According to the mathematical description of steps 1 and 2, the visibility problem of the satellite's complex conical field of view to regional targets can be expressed as: whether the visible area projected by the satellite's complex cone field of view on the earth's surface overlaps with the regional target.

Step 4. In order to reduce the calculation amount of the problem of "determining whether a region is inside another region", an idea of calculating the "limit point" is proposed.

For the complex-shaped projection field of view area A and area target B, the traditional method to determine whether area A is inside area B is: sequentially follow the boundaries of area B, traverse all points in area A, and determine whether they are all in the area. inside each edge of B. This processing method is computationally intensive. Therefore, the present invention proposes the concept of "limit point", whose essential meaning is to find the point that is most likely not to satisfy the situation that "area A is located inside area B". This point is called a "limit point". By judging the "limit point" "The positional relationship with area B is used to determine whether the positional relationship between the two areas satisfies the situation of "area A is located inside area B". Its meaning is shown in Figure 5. That is, to determine whether the area A shown in the figure is on a certain side in the direction of the arrow, you only need to verify one point shown in the figure. This point is the "limit point" described in the present invention. The principle is similar for the case of a three-dimensional sphere. It should be noted that the location of the "limit point" is related to the direction of the boundary of area B (i.e., the direction of the arrow in Figure 5). Therefore, the "limit point" of area A corresponding to each edge of area target B should be calculated independently.

Step 5. To calculate the overlapping problem of two areas on the sphere, establish a basic algorithm model of vector geometry.

By transforming the problem into an overlapping problem of two areas on a spherical surface, the use of spherical space vector calculations will significantly reduce the amount of calculation. To this end, the present invention proposes the following five vector geometry basic algorithms for rapid calculation of relevant required information:

(51) Find the vector of a spherical vector rotated along a certain direction;

This algorithm is designed to solve the "limit point" of the complex cone field of view relative to the target boundary of each area. As shown in Figure 2, if you find the vector p' after rotating α to v' (unit vector), then p' can be expressed by the following formula: p'=cos(α)p+sin(α)v ^* , where v ^* =norm(v'-(v'·p)v'), norm(·) represents the normalization of the vector.

(52) Determine whether a point on the sphere is on a certain side of the great circle;

This algorithm is designed to determine the relationship between each point and the target boundary of each area. as shown in picture 2. If it is judged that p is on one side of the great circle 2. If p·v≥0, then p is on the "same side" of the normal vector of the great circle 2, otherwise, it is on the "outside of the normal vector" of the great circle

(53) Determine whether a point on the spherical surface is on a certain side of an ordinary circle;

This algorithm is designed to determine the relationship between the target boundary points in each area and the ordinary circle of the inner and outer conic projections of the complex conical field of view. as shown in picture 2. To determine whether p is on a certain side of an ordinary circle. Similarly, calculate the size of p·v'. If p·v≥cos(θ), it means that p is on the "same side of the normal vector" of the ordinary circle, otherwise, it is on the "outside of the normal vector" of the circle.

(54) Determine whether two large arcs have an intersection;

This algorithm is designed to determine whether there is an intersection point between the target boundary of each area and the great circle of the largest and small clock angle projection of the complex cone field of view. The general idea of determining whether two great arcs have an intersection can be summarized as the following two points: first, calculate the intersection of two great circles; second, determine whether the point is on two great arcs.

Taking Figure 2 as an example, first, the intersection point Q of the big circle 2 and the big circle 1 can be calculated according to the following formula Q = ±norm (v×v'). Since there are two intersection points, there are two positive and negative situations in the formula. Then, determine whether Q is inside the great arc d ₁ d _2. At this time, you need to determine (d ₁ ×Q)·(Q×d ₂ )≥0 and Q·(d ₁ +d ₂ )≥0. If the above conditions are met Two conditions, it is proved that Q is inside the great arc d ₁ d ₂ .

(55) Determine whether a large arc and an ordinary arc have intersection points.

This algorithm is designed to determine whether there is an intersection point between the target boundary of each area and the ordinary circle of cone projection inside and outside the complex cone field of view. The calculation process is the same as algorithm (54). The following is a detailed description using the side view perpendicular to the v and v' directions in Figure 2. As shown in Figure 3. Among them, the intersection point of the ordinary circle and the great circle 2 is Q', and the angle between v and v' is β.

Similarly, the intersection point Q' is required first, but it is necessary to first determine whether the intersection point exists. This is because there are three possibilities for the intersection point of an ordinary circle and a great circle.

If β+θ=90°, then the ordinary circle is tangent to the great circle;

If β+θ>90°, then the ordinary circle intersects the great circle.

If β+θ<90°, then there is no intersection between the ordinary circle and the great circle.

Among them, β is calculated as follows β = arccos (|v·v'|), which is the two cases considering the positive and negative normal vector of the great circle.

When it is determined that there is an intersection, its vector needs to be calculated. Assume that the angle between the intersection point Q' and the current projection plane is γ, then cos(γ)=cos(θ)/sin(β). At this time, the direction vector in the constructed plane is v _// =norm(v'-(v'·v)v) and the unit vector v _⊥ =norm(v'×v) in the vertical plane.

Then the intersection point Q' can be calculated according to the following formula: Q'=cos(γ)v _// ±sin(γ)v _⊥

Finally, it is necessary to determine whether the intersection point is within two arcs. You only need to project Q' onto the ordinary circular plane first. The subsequent ideas are consistent with algorithm (54) and will not be repeated here.

Step 6. Algorithm preprocessing:

(61) Calculation of normal vector of area target boundary

In order to facilitate subsequent calculations, the boundary points of the regional targets are sorted counterclockwise, and their coordinates are converted to the ground-fixed system. The definition can be easily transformed and normalized to obtain DΝ ₁ , DΝ ₂ ,. ..,DΝ _n (DistrictNode).

The obtained coordinates DΝ ₁ , DΝ ₂ ,..., DΝ _n are cross-multiplied with each other in order to obtain the normal vector of each boundary. In order to obtain the geodesic formed by the boundary points of two adjacent areas, it is necessary to cross-multiply the boundary point vectors to obtain its normal vector, the expression is:

And do normalization processing:

The regional boundary normal vectors are obtained as NV ₁ , NV ₂ ,..., NV _n .

(62) Calculate the half-cone angle of the inner and outer cones of the earth corresponding to the satellite field of view.

In the process of finding the "limit point", it is necessary to use the half cone angle of the earth's field of view to the satellite. For general earth observation satellites, its eccentricity is close to zero. Therefore, during the entire satellite operation, the half-cone angle of the Earth's field of view remains almost unchanged. In order to reduce the amount of calculation, during the preprocessing process of the regional target, this algorithm calculates the half cone angle in advance.

The half-cone angle of the earth's field of view to the satellite can be represented by Figure 4. Among them, D represents the distance of the satellite from the center of the earth, R represents the radius of the earth, and θ represents the half-cone angle of the satellite's field of view. Represents the desired half-cone angle of the earth’s field of view to the satellite, and its calculation formula is shown in Figure 4.

For complex conical fields of view, it is necessary to calculate the Earth's half-cone angle corresponding to the inner and outer cone half-angles. and/>

Step 7. Calculate the "limit point" of the complex conic projection field of view

For a complex conical field of view, due to the complex boundaries of its projection area, in order to facilitate subsequent unified processing, the present invention combines the shape of its projection area and selects at most 6 "limit points" in one judgment, among which points 1, 2, 3 and 4 are the four boundary points of the complex cone field of view projection area A, while points 5 and 6 are calculated based on the situation of the regional target boundary (which may not exist). The schematic diagram of its "limit point" is shown in Figure 6.

(71) First, calculate the Earth’s half-cone angle corresponding to the inner and outer cone angles. and/> For details on the algorithm, see Step 6 of "Calculate the half-cone angle of the inner and outer cones of the earth's corresponding satellite field of view";

(72) Secondly, calculate the “limit points” 1, 2, 3, 4, 5, 6 of the complex cone field of view projection area in Figure 7.

Calculate the unit vectors n _x and n _y of the unitized velocity and reverse angular velocity as the unit vectors of the plane, as shown in Figure 7.

After obtaining the unit vector, the unit direction vector corresponding to each "clock side" is obtained through the definition of the clock angle. Taking the "clock side" where points 1 and 3 are located in Figure 7 as an example, the expression of the unit direction vector v _1,3 as follows

Rotate the position vector of the center of the circle where the complex cone field of view projection area is located to the v _{1 and 3} directions respectively. horn, The position vector of 1 point and the position vector of 3 points can be obtained by using the angle. For the rotation method, please refer to the algorithm (51) in step 5, "Finding the vector of a vector on the spherical surface after rotating the angle along a certain direction." According to the above method, the position vectors of points 1, 2, 3, and 4 can be obtained.

For the possible 5 and 6 points, the calculation method is the same, except that the direction vector is the normal vector of the regional target boundary. It should be noted that since the projected area of the complex conical field of view is not a complete ring, the points obtained by rotating using the above method may not exist in the effective visible area of the complex conical field of view, so its existence needs to be determined.

First, we need to calculate the midpoint vector v _1/2 of the projection area. The calculation formula is as follows:

v _1/2 =f(norm(v _1,3 +v _2,4 ))

The norm function represents unitization, and the expression of the f function is as follows:

If the position vector of the 5 points obtained according to the original algorithm is p ₅ , and if it satisfies the following relationship, it proves that the point is within the field of view, that is, the point exists, otherwise the point does not exist.

Step 8. Determine whether there is an overlapping area between the two areas at a certain time based on the visibility judgment analysis algorithm.

One of the foundations of solving the visible window is to be able to judge the visibility of regional targets at each moment. According to the basic algorithm of vector geometry that has been proposed, the present invention proposes to divide the visible situation into three situations for the mathematical model of this problem, and judge whether it is satisfied for each situation in turn. If one situation is satisfied, it can be judged that the moment is visible. If the visible area of the complex cone field of view projected on the earth's surface is called area A, and the area target is called area B, then the above three situations are: (a) A is within B, (b) A and B have boundaries The intersection point, (c) B is within A. The algorithm flow of these three cases is as follows:

(81) Determine whether it belongs to case (a), that is, whether the conical field of view is located in the regional target:

(811) Traverse each edge of area target B and calculate the "limit point" of area A relative to this edge. For details on how to calculate the “limit point”, see step 7;

(812) For each edge of the regional target, determine whether all "limit points" are inside the edge. For the specific method of determining whether the point is inside the edge, please refer to the algorithm (52) in step 5 "Judge whether a point on the sphere is on the edge." One side of the great circle". If there is a limit point that is not inside a certain edge, exit, return "No", and enter step (82).

(813) If after traversing each edge, all "limit points" are inside each edge, then return "Yes".

Among them, the "limit points" of the complex cone field of view are the special points searched for in case (a). The purpose is to obtain the relative relationship between the projection area of the complex cone field of view and the target area by proving the positional relationship between these points and the target area. positional relationship to reduce the amount of calculation. In addition to the four points on the boundary of the complex cone field of view projection area, the "limit points" can be divided into two categories according to whether the difference between the maximum and minimum clock half angles is greater than 180 degrees, as shown in Figure 7. When it is less than 180 degrees, it is necessary to judge whether there is point 5 based on the situation of the regional target boundary; when it is greater than 180 degrees, it is necessary to judge whether there are points 5 and 6 based on the situation of the regional target boundary.

(82) Determine whether it belongs to case (b), that is, determine whether there is an intersection between the complex conical field of view projection area and the target area:

Determine in turn whether each boundary of the regional target intersects with the boundary of the complex cone field of view in the earth's projection area. If there is an intersection, then case (b) is met and "yes" is returned; if not, "no" is returned. For specific methods to determine whether there is an intersection, please refer to the algorithm (54) "To determine whether two large arcs have an intersection" and algorithm (55) "To determine whether a large arc and an ordinary arc have an intersection" in step 5.

(83) Determine whether it belongs to case (c), that is, whether the target area is located within the complex cone field of view projection area:

Determine in turn whether the boundary points of each area target are within the projection area of the complex cone field of view. If all points are satisfied, then case (c) is satisfied and "Yes" is returned; if one point is not satisfied, "No" is returned. . For details of the specific algorithm, please refer to the algorithm (52) of step 5 "Judge whether a point on the spherical surface is on a certain side of a great circle" and (53) "Judge whether a point on a spherical surface is on a certain side of an ordinary circle".

Among them, being located inside the complex conic field of view projection area means being on the "different side of the normal vector" of the ordinary circle of the inner cone projection, and being on the "same side of the normal vector" of the ordinary circle of the outer cone projection and other projection boundaries.

Step 9. Semi-analytically determine the visible window of the complex cone field of view

Through the visibility judgment analysis algorithm proposed above, a semi-analytic method for final determination of the visible window is obtained. First, determine the visibility at that moment based on the discrete step length; then arrange the visibility at each step length in chronological order and determine the visible start and end time ranges of each visible window; and then obtain each visible window through binary search. The fine moment of the visible window, in which the satellite orbit information used in the binary search is obtained by Lagrangian interpolation; finally, the precise visible window information is summarized and output. The overall operation process can be explained according to the flow chart shown in Figure 8.

To sum up, the above are only preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention shall be included in the protection scope of the present invention.

Claims (8)

1. A method for calculating the visible window of a regional target in a complex conical field of view of a remote sensing satellite, which is characterized by including the following steps: Step 1. Establish a mathematical description of the complex cone field of view: The projections of the circle corresponding to the inner half-cone angle and the circle corresponding to the outer half-cone angle of the complex cone field of view on the earth's surface are two ordinary arcs. The minimum clock angle and The projection of the two sides corresponding to the maximum clock angle on the earth's surface is two great arcs; among them, the center of the great arc coincides with the center of the sphere, and the circumference is on the earth's surface; an ordinary arc is the circumference of the circle on the earth minus the great arc. Any arc outside; Step 2. Establish a mathematical description of the regional target: the regional target refers to a closed convex polygon on the earth's surface, and each side of the convex polygon is a great arc on the earth's sphere; Step 3. Determine whether there is an overlapping area between the visible area A of the complex cone field of view projected on the earth's surface at each moment and the area B where the regional target is located. Specifically: (81) Determine whether area A is within area B: (811) Traverse each edge of area B and calculate the limit point of area A relative to the edge; where the limit point refers to the point on the boundary of area A that is most likely to be outside area B; each boundary edge of area B is Corresponds to a set of limit points; (812) For each edge of area B, determine whether its corresponding limit points are all inside the edge. If there is a limit point that is not inside the edge, exit, return "No", and enter step (82); (813) If all the limit points on area A are inside each corresponding side of area B, then area A is within area B at that moment; (82) Determine whether there is at least one intersection point on the boundary between area A and area B. If so, then the two areas have an intersection point; if not, exit, return "No", and enter step (83); (83) Determine whether area B is within area A: Determine in turn whether each boundary point of area target B is within the projection area A of the complex cone field of view. If all points are satisfied, then area B is within area A; otherwise, if one point is not satisfied, exit and return " no"; For a certain moment, as long as one of the above three conditions (81), (82), and (83) is satisfied, the field of view at that moment is visible to the regional target; if none of the above three conditions are satisfied, then the field of view at that moment is Invisible to area targets; Step 4. Based on the positional relationship between area A and area B at each time calculated in step 3, determine the visible window of the complex cone field of view to area target B. 2. A method for calculating the regional target visible window of a complex conical field of view of a remote sensing satellite as claimed in claim 1, characterized in that in the step (811), the method for calculating the limit point is: First, calculate the Earth's half-cone angle corresponding to the inner and outer cone half-angles of the complex cone field of view. and/> Where D represents the distance of the satellite from the center of the earth, R represents the radius of the earth, θ represents the half-cone angle of the satellite's field of view,/> Represents the desired half-cone angle of the earth’s field of view to the satellite: Secondly, calculate the "limit points": take up to 6 "limit points" in one judgment, among which points 1, 2, 3, and 4 are the 4 boundary points of area A, while points 5 and 6 are based on the area To calculate the target boundary, the specific steps are as follows: Calculate the unit vectors n _x and n _y of the unitized speed and reverse angular velocity, and obtain the unit direction vector corresponding to each "clock side" through the definition of the clock angle. Taking the "clock side" where points 1 and 3 are located as an example, The expression of its unit direction vector v _1,3 is as follows: Represents the minimum clock angle; Using the rotation method, the position vector of the center of the circle where the complex cone field of view projection area is located is rotated to the v _{1, 3} directions respectively. angle,/> The position vector of point 1 is obtained in the same way as the position vector of point 3. Use this method to obtain the position vectors of points 2 and 4; Calculate the 5th and 6th limit points that may exist on the boundary line of complex conical field of view projection area A relative to area B: among them, when the difference between the maximum and minimum clock half angles is less than 180 degrees, the fifth limit point may exist; when it is greater than 180 degrees When the degree is reached, there are the 5th and 6th limit points; For the possible 5th and 6th limit points on the boundary of area A, use the normal vector of the area target boundary as its direction vector, and use the calculation method with points 1, 2, 3, 4 to calculate the position vector of the possible limit points; then according to The following formula determines whether the possible limit point is in area A: If it satisfies the above formula, it proves that the possible limit point is within the field of view, that is, the limit point exists, otherwise the limit point does not exist; Among them; p ₅ is the position vector of the possible limit point, the norm function represents the unitization, and the expression of the f function is as follows: Represents the maximum clock angle; v _1/2 = f (norm (v _1,3 + v _2,4 )) represents the midpoint vector of area A; v _2,4 represents the unit of the "clock side" where points 2 and 4 are located direction vector. 3. A method for calculating the regional target visible window of a complex conical field of view of a remote sensing satellite as claimed in claim 2, characterized in that the rotation method is as follows: Find the vector p' after rotating α to the unit vector v'. p' is represented by the following formula: p'=cos(α)p+sin(α)v ^* , where v ^* =norm(v'-(v'· p)v'), norm(·) represents the unitization of the vector; the unit vector v' represents the unit direction vector of a great arc or an ordinary arc. 4. A method for calculating the regional target visible window of a complex conical field of view of a remote sensing satellite as claimed in claim 1, characterized in that in the step (812), it is determined whether the corresponding limit points are all on the inside of the side. The method is: Let v represent the normal vector of each side of area B, and let p represent the position vector of a certain point. If p·v≥0, then the vector p is on the same side as the normal vector of the side, otherwise, it is outside the normal vector of the side. 5. A method for calculating the regional target visible window of a complex conical field of view of a remote sensing satellite as claimed in claim 1, characterized in that in step (82), when judging whether there is an intersection between area A and area B, the method is as follows: The method determines whether there is an intersection point between the large arc where the target boundary of each area is located and the large arc projected by the maximum and small clock angles of the complex cone field of view, specifically as follows: The intersection point Q of the above two large arcs is calculated according to the following formula Q=±norm(v×v'), and it is judged whether Q is inside the large arc d ₁ d _2. At this time, it is necessary to judge whether (d ₁ ×Q)· (Q×d ₂ ) ≥ 0 and Q · (d ₁ + d ₂ ) ≥ 0. If satisfied, it proves that Q is inside the great arc, that is, the above two great arcs have an intersection. 6. A method for calculating the regional target visible window of a complex conical field of view of a remote sensing satellite as claimed in claim 1, characterized in that in the step (82), when judging whether there is an intersection between area A and area B, the method is as follows: The method determines whether there is an intersection point between the large arc at the target boundary in each area and the ordinary arc of the cone projection inside and outside the complex cone field of view, specifically as follows: Suppose the intersection point of the ordinary arc and the great arc is Q', and let the angle between the unit vectors v and v' be β; where v and v' are respectively the unit normal vector of the plane where the great arc is located and the location of the ordinary arc. The unit normal vector of the plane; If β+θ=90°, then the ordinary arc is tangent to the great arc; θ represents the half-cone angle of the satellite field of view corresponding to the ordinary arc; If β+θ>90°, then the ordinary arc intersects with the great arc; If β+θ<90°, there is no intersection between the ordinary arc and the great arc; Among them, β is calculated as follows β=arccos(|v·v'|); When it is determined that there is an intersection, the vector needs to be calculated: assuming that the angle between the intersection Q' and the current projection plane is γ, then cos(γ)=cos(θ)/sin(β); then the direction vector in the constructed plane is v _// =norm(v'-(v'·v)v) and the unit vector v in the vertical plane _⊥ =norm(v′×v); Then the intersection point Q' is calculated according to the following formula: Q'=cos(γ)v _// ±sin(γ)v _⊥ ; Finally, it is necessary to determine whether the intersection point is within the above-mentioned large arc and ordinary arc. First, project Q' onto the ordinary circle plane, and then judge whether there is an intersection point between the circles. 7. A method for calculating the regional target visible window of a complex conical field of view of a remote sensing satellite as claimed in claim 1, characterized in that the specific method of step 4 is: First, based on the discrete step length, determine the visibility of the remote sensing satellite to the regional target at that moment; then arrange the visibility at each step length in chronological order and determine the visible start and end time ranges of each visible window; and then The precise moment of each visible window is obtained through binary search, and finally the information of each precise visible window is summarized and output. 8. A method for calculating the regional target visible window of a complex conical field of view of a remote sensing satellite as claimed in claim 7, wherein the satellite orbit information used in the binary search is obtained by Lagrangian interpolation.