CN109993719B - Multi-rail splicing imaging optimization method for area coverage - Google Patents

Multi-rail splicing imaging optimization method for area coverage Download PDF

Info

Publication number
CN109993719B
CN109993719B CN201910214391.XA CN201910214391A CN109993719B CN 109993719 B CN109993719 B CN 109993719B CN 201910214391 A CN201910214391 A CN 201910214391A CN 109993719 B CN109993719 B CN 109993719B
Authority
CN
China
Prior art keywords
yaw angle
angle
polygon
roll
area
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910214391.XA
Other languages
Chinese (zh)
Other versions
CN109993719A (en
Inventor
沈欣
张过
李仕学
刘钰霖
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Wuhan University WHU
Original Assignee
Wuhan University WHU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Wuhan University WHU filed Critical Wuhan University WHU
Priority to CN201910214391.XA priority Critical patent/CN109993719B/en
Publication of CN109993719A publication Critical patent/CN109993719A/en
Application granted granted Critical
Publication of CN109993719B publication Critical patent/CN109993719B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/50Image enhancement or restoration using two or more images, e.g. averaging or subtraction
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10032Satellite or aerial image; Remote sensing
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20212Image combination
    • G06T2207/20221Image fusion; Image merging

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Image Analysis (AREA)

Abstract

The application provides a multi-track splicing imaging optimization method facing area coverage, and the method comprises the following steps of firstly, constructing n side swing angle scheme sets P formed by combining m side swing angles, and correspondingly generating a set F formed by n coverage rates, wherein m is the number of tracks, and each side swing angle value is randomly generated between the upper limit and the lower limit; calculating the coverage rate of each side sway angle scheme in the P by adopting a polygon logic calculation method, and recording the maximum value of the coverage rate and the corresponding side sway angle scheme; adopting a PSODE algorithm, carrying out evolution operations such as crossing and variation on each side swing angle scheme in the scheme P, and then reserving the scheme with better coverage rate to obtain a corresponding new element set P'; and after the last operation is repeatedly executed for k times, taking the yaw angle scheme with the maximum coverage rate in the P' as an optimal solution. The invention overcomes the defects of large calculation time consumption and low efficiency of the traditional grid point method.

Description

Multi-rail splicing imaging optimization method for area coverage
Technical Field
The invention relates to the technical field of satellite remote sensing, in particular to a multi-orbit splicing imaging optimization method facing to area coverage.
Background
Because the satellite orbit is low in height, the load field of view is limited, the imaging width is small, the side swing imaging is needed, and the target imaging in a large-range area can be realized in a multi-strip splicing mode. The sidesway angle during each satellite transit observation is required to be optimized within a given task duration, and the requirement of area coverage is met to the maximum extent.
In the existing regional imaging multi-strip yaw angle optimization research, two common ideas are a grid point method and a strip method. The grid point method is characterized in that a regional target is converted into a point target set, and a yaw scheme is formed by optimizing the coverage quantity of grid points, so that the modeling process can be effectively simplified, the task preprocessing workload is large, and the optimization effect depends on the density of the grid points. The strip method firstly decomposes a target area into a plurality of imaging strips parallel to a satellite motion trail, then uses a search algorithm to solve the optimal yaw angle of each strip by taking the maximum coverage area as a target. For the different stripe dividing modes, two types can be summarized: 1) based on a fixed-width stripe division method, according to the flight radial direction of a satellite and the width of a satellite-borne sensor, a region is divided into parallel stripes with fixed width, if the whole division is controlled by adopting parameters of direction and offset, and a group of imaging stripes with overlapped edges are used for covering a region target, the stripe division method is simple to realize, but the change of the imaging width under the condition of different side swing angles is not considered, and the model precision is limited; 2) based on the method of regional dynamic division, in each imaging time window of a satellite to a regional target, a group of imaging strips with equal intervals are constructed according to the satellite sidesway angle to serve as candidate strips, and the corresponding sidesway angle is searched in a certain step length.
Aiming at the regional imaging task requirement and the defects of the existing method, an optimization method for solving a sidesway angle optimization model is provided.
Disclosure of Invention
The invention provides a region coverage rate rapid calculation method based on vector polygon logical operation and a PSODE algorithm based on region coverage rate calculation and a sidesway angle optimization model aiming at the region coverage rate calculation and the sidesway angle optimization model on the basis of constructing a region imaging side swing angle optimization model, and provides a satellite sidesway angle optimization method based on the PSODE algorithm and facing to a region imaging task.
An area coverage-oriented multi-rail splicing imaging optimization method comprises the following steps:
step S1, randomly selecting 1 yaw angle between the upper limit and the lower limit of the yaw angle of each orbit arc segment aiming at m orbit arc segments of the satellite to form a yaw angle combination; wherein m is a positive integer;
step S2, repeating step S1 n times to obtain n yaw angle combinations, and recording the n yaw angle combinations as a set P, wherein n is a positive integer;
step S3, determining the coverage rate corresponding to each element in the set P, and recording the coverage rate corresponding to each element in the set P as a set F;
step S4, taking a yaw angle combination corresponding to the maximum coverage rate in the set F as an optimal yaw angle combination;
step S5, updating, mutating and crossing each element in the set P to obtain a set P' consisting of corresponding new elements;
step S6, obtaining a set F 'according to the coverage rate corresponding to each element in the set P';
step S7, comparing each element of the set F 'with the corresponding element in the set F, and if each element of the set F' is larger than the corresponding element in the set F, replacing the corresponding element in the set F with the element of the set F ', and replacing the corresponding element in the set P with the element of the set P';
step S8, comparing the largest element in the set F obtained in step S7 with the coverage rate corresponding to the optimal roll angle combination in step S4, and if the largest element in the set F obtained in step S7 is greater than the coverage rate corresponding to the optimal roll angle combination in step S4, taking the element in the set P corresponding to the largest element in the set F obtained in step S7 as the optimal roll angle combination to be optimized;
step S9, replacing the set P and the set F in step S5 with the set P and the set F obtained in step S7, and repeating steps S5 to S8K times to obtain a final optimal yaw angle combination, where K is an evolutionary algebra and K is an integer greater than or equal to 1.
Further, determining the coverage corresponding to each element in the set P includes:
determining a rectangular strip corresponding to each yaw angle according to each yaw angle in a yaw angle combination, a central projection imaging geometric relation and the position and the speed of a satellite, and obtaining vertex coordinates of four vertexes of each rectangular strip;
solving all the rectangular strips according to the vertex coordinates of all the rectangular strips to obtain a covering polygon;
intersecting the coverage polygon with a target area to obtain an effective coverage polygon;
and dividing the area of the effective coverage polygon by the area of the target area to obtain the coverage rate.
Further, the rectangular strip is formed on the ground surface through four intersection points of characteristic rays from the imaging center S corresponding to each yaw angle to four corner points of the area array sensor and the ground surface.
Further, the rectangular strip is formed by scanning the characteristic light rays from the imaging center S corresponding to each side swing angle to two corner points of the linear array sensor on the ground surface for a period of time.
Further, the vertex coordinates of each vertex of each rectangular strip are determined according to the following steps:
in an inertial coordinate system, the earth is taken as a standard rotation ellipsoid, and the equation of the standard rotation ellipsoid is a first formula:
Figure BDA0002001557820000031
wherein, (x, y, z) is a point on the surface of the standard ellipsoid of revolution, and a and b are respectively the length of the major semi-axis and the minor semi-axis of the standard ellipsoid of revolution;
the equation of the characteristic light ray where the vertex is located in the geocentric inertial coordinate system is a second formula:
Figure BDA0002001557820000032
wherein (x)0,y0,z0) The coordinate of the imaging center S under the geocentric inertial coordinate system is shown, and the (A, B and C) are the direction numbers of the characteristic light rays where the vertex is located under the geocentric inertial coordinate system;
and combining the first formula and the second formula to obtain the vertex coordinates of each vertex of each rectangular strip.
Further, the area of the effective coverage polygon is obtained according to the following steps:
dividing the effective coverage polygon into a plurality of triangles, solving the area of each triangle one by adopting a vector product method, and summing the areas of all the triangles to obtain the area of the effective coverage polygon; wherein, the calculation formula of the area of the effective coverage polygon is as follows:
Figure BDA0002001557820000041
wherein x isi、yiIs the plane coordinate of the ith vertex, q is the number of the vertices of the effective coverage polygon, and S is the area of the effective coverage polygon, wherein i and q are positive integers.
Further, the unit of the vertex coordinates of the four vertices of each rectangular strip is mm;
rounding the vertex coordinates of all the rectangular strips;
and solving all the rectangular strips according to the rounded vertex coordinates of all the rectangular strips to obtain a covering polygon.
Further, the intersection rule of intersection includes an Unlike intersection rule and a like intersection rule.
Further, the upper and lower limits of the yaw angle of each track arc segment are determined according to the following steps:
obtaining the track distance d between each vertex of the rectangular strip and the satellite point according to the position relation between the satellite point track of each orbit arc segment of the satellite and the rectangular strip;
determining a yaw angle beta required when a satellite main optical axis points to each vertex for imaging according to the satellite orbit height h and the camera field angle IFOV:
β=atan(d/h)+IFOV
sequencing the sidesway angles beta corresponding to the vertexes, wherein the maximum value of the sidesway angle is the maximum sidesway angle u of the track arc section, and the minimum value of the sidesway angle is the minimum sidesway angle l of the track arc section;
let the range of the satellite yaw angle be (-Roll)max,Rollmax);
If l is<-RollmaxAnd u>Rollmax, the upper and lower limits of the side swing angle of each track arc section are respectively-Rollmax and Rollmax;
if l < -Rollmax and u < Rollmax, the upper and lower limits of the yaw angle of each track arc section are respectively-Rollmax and u;
if l > -RollmaxAnd u < RollmaxThe upper and lower limits of the yaw angle of each track arc section are l and u respectively;
if l > -RollmaxAnd u>RollmaxThe upper and lower limits of the yaw angle of each track arc segment are l, Rollmax
Further, the roll angle optimization model for obtaining the final optimal roll angle combination in step S9 is as follows:
Figure BDA0002001557820000051
wherein, P ═ { x ═ x1,x2,…,xi,…xn},xi∈[li,ui]In the formula xiFor a yaw angle combination, Scov (x)i) To cover the area of the polygon effectively, SobjIs the area of the target region, n is the number of elements in the set P, li、uiAnd i and n are positive integers, which are the upper and lower limits of the yaw angle of the ith track arc section.
Compared with the prior art, the method has the following advantages:
according to the invention, the satellite imaging process is decomposed into a plurality of orbit arc segment imaging processes, and a side swing angle is randomly selected on each orbit arc segment in consideration of the change of the imaging width under the condition of different side swing angles, different side swing angles form different rectangular strips, and the widths of the different rectangular strips are different, so that the plurality of rectangular strips with different widths in the application are different from the strips which only use a single type and have fixed widths in the prior art. The method overcomes the defect that the traditional method based on the division of the strips with fixed width does not take the change of the imaging width under different side swing angle conditions into consideration, further adopts a PSODE algorithm to continuously update, cross and mutate n side swing angle combinations to generate new side swing angle combinations, obtains the side swing angle combinations meeting the imaging conditions from the new side swing angle combinations, and improves the accuracy of determining the optimal side swing angle combinations according to the final optimal side swing angle combinations meeting the imaging conditions.
Since the vector polygon logical operation method needs to be implemented by integer numbers, vertex coordinates of the covered polygons need to be converted into integer numbers, and the calculation precision depends on the unit of floating point numbers and rounding errors of the integer numbers converted from the floating point numbers. In order to ensure the accuracy of the solution of the covering polygon, the invention converts the vertex coordinate unit of the covering polygon into mm and then carries out rounding, and ensures that the coordinate accuracy loss caused by rounding is less than 1 mm.
Drawings
FIG. 1 is a diagram of the imaging range of a characteristic light and area array sensor of the present invention on a single track arc segment;
FIG. 2 is a schematic view of a polygon covered by the area array sensor of the present invention;
FIG. 3 is a schematic diagram of Vatti algorithm polygon reconstruction and scan line dynamic update according to the present invention;
FIG. 4 is a schematic diagram of a polygon intersection rule according to the present invention;
FIG. 5 is a graph of the upper and lower bounds of the roll angle of the present invention;
FIG. 6 is a diagram of the coverage rate evolution process with iterative algebra according to the present invention;
FIG. 7 is a diagram illustrating the coverage effect obtained from step S1 to step S4 according to the present invention;
FIG. 8 is a diagram illustrating the coverage effect obtained from step S5 to step S9 according to the present invention;
FIG. 9 is a flow chart of the steps of the present invention.
Detailed Description
The technical scheme of the invention is explained in detail in the following by combining the drawings and the embodiment.
In order to overcome the defects of large calculation time consumption and low efficiency of the traditional grid point method, the invention adopts a coverage rate calculation method based on vector polygon logical operation, which comprises the following steps: the method comprises the steps of firstly calculating the imaging strip boundary corresponding to each orbit arc segment of the satellite according to the imaging geometric relation, then carrying out logic operation on a coverage area formed by a plurality of strips and a target area to obtain the boundary of an effective coverage area, and obtaining the coverage rate. And finally, optimizing the satellite yaw angle by taking the obtained coverage rate as an objective function based on a PSODE algorithm.
As shown in fig. 9, a multi-track stitching imaging optimization method facing area coverage includes the following steps:
step S1, randomly selecting 1 yaw angle between the upper limit and the lower limit of the yaw angle of each orbit arc segment aiming at m orbit arc segments of the satellite to form a yaw angle combination; wherein m is a positive integer;
wherein, the upper and lower limits of the sidesway angle of each track arc section are determined according to the following steps:
as shown in fig. 5, the distance d between each vertex of the rectangular strip and the track of the satellite point is obtained according to the position relationship between the track of the satellite point under the satellite and the rectangular strip of each orbital arc segment of the satellite;
determining a yaw angle beta required when a satellite main optical axis points to each vertex for imaging according to the satellite orbit height h and the camera field angle IFOV:
β=atan(d/h)+IFOV
sequencing the sidesway angles beta corresponding to the vertexes, wherein the maximum value of the sidesway angle is the maximum sidesway angle u of the track arc section, and the minimum value of the sidesway angle is the minimum sidesway angle l of the track arc section;
let the range of the satellite yaw angle be (-Roll max, Roll)max);
If l is<-RollmaxAnd u>RollmaxThe upper and lower limits of the yaw angle of each track arc segment are respectively-Rollmax,Rollmax
If l is<-RollmaxAnd u < RollmaxThe upper limit and the lower limit of the side swing angle of each track arc section are-Rollmax, u respectively;
if l is-Rollmax and u is less than Rollmax, the upper and lower limits of the yaw angle of each track arc section are l and u respectively;
if l > -RollmaxAnd u>RollmaxThe upper and lower limits of the yaw angle of each track arc segment are respectivelyl,Rollmax
Step S2, repeating step S1 n times to obtain n yaw angle combinations, and recording the n yaw angle combinations as a set P, wherein n is a positive integer;
since the roll angles in the elements (i.e., the roll angle combinations) in the set P are all randomly generated, the case where two elements in the set P are the same may occur, but the case where two elements in the set P are the same has a very small occurrence probability, so that the case can be ignored here. Even if two elements in the set P are identical, there is no influence on the present embodiment.
Step S3, determining the coverage rate corresponding to each element in the set P, and recording the coverage rate corresponding to each element in the set P as a set F;
determining the coverage rate corresponding to each element in the set P, including:
and step 31, determining a rectangular strip corresponding to each yaw angle according to each yaw angle in one yaw angle combination, the central projection imaging geometric relationship and the position and the speed of the satellite, and obtaining the vertex coordinates of four vertexes of each rectangular strip.
The method for determining the rectangular strip includes, but is not limited to, the following two methods:
mode 1:
the rectangular strips are formed on the ground surface through four intersection points of characteristic rays from the imaging center S corresponding to each side swing angle to four corner points of the area array sensor and the ground surface.
Mode 2:
the rectangular strip is formed by scanning characteristic light rays from an imaging center S corresponding to each side swing angle to two angular points of the linear array sensor on the ground surface for a period of time.
Wherein the vertex coordinates of each vertex of each rectangular strip are determined according to the following steps:
in an inertial coordinate system, the earth is taken as a standard rotation ellipsoid, and the equation of the standard rotation ellipsoid is a first formula:
Figure BDA0002001557820000071
wherein, (x, y, z) is a point on the surface of the standard ellipsoid of revolution, and a and b are respectively the length of the major semi-axis and the minor semi-axis of the standard ellipsoid of revolution;
the equation of the characteristic light ray where the vertex is located in the geocentric inertial coordinate system is a second formula:
Figure BDA0002001557820000081
wherein (x)0,y0,z0) The coordinate of the imaging center S under the geocentric inertial coordinate system is shown, and the (A, B and C) are the direction numbers of the characteristic light rays where the vertex is located under the geocentric inertial coordinate system;
and combining the first formula and the second formula to obtain the vertex coordinates of each vertex of each rectangular strip.
As shown in fig. 1, according to the central projection geometry, the satellite position, the velocity and the yaw angle, the intersection points of the characteristic rays from the imaging center S to the four angular points of the area array sensor and the ground surface are determined (for the line array sensor, the intersection points of the characteristic rays from the imaging center S to the two angular points of the line array sensor and the ground surface are determined), which are P respectively1、P2、P3、P4(ii) a Obtaining a rectangular strip, P, formed on the surface of the earth for each orbital arc of the satellite1、P2、P3、P4Forming a rectangle on the ground surface, and then obtaining the vertex P of each rectangular strip1、P2、P3、P4The coordinates of (a). As shown in fig. 2, the coverage area of the area array sensor in one step is an approximately rectangular area formed by P1 and P2 corresponding to the start point of the step and P3 and P4 corresponding to the end point of the step.
The characteristic light rays are four homonymous light rays corresponding to the outermost side of the instantaneous field angle of the sensor, and determine the imaging range of a rectangular strip (an approximate rectangle consisting of P1, P2, P3 and P4 in FIG. 1) formed on the ground surface by the sensor at each orbital arc segment of the satellite according to a collinear condition.
The coordinates of P1, P2, P3 and P4 can be obtained by the intersection point of the characteristic ray and the earth surface. In an inertial coordinate system, the earth is taken as a standard rotation ellipsoid, and the equation is as follows:
Figure BDA0002001557820000082
in the above formula, a and b are respectively the length of the major semi-axis and the minor semi-axis of the standard ellipsoid of revolution. According to the WGS84(World Geodetic System 1984, WGS84) earth model, a 6378.137km and b 6356.752 km.
The equation of the characteristic ray in the inertial coordinate system is as follows:
Figure BDA0002001557820000083
in the formula (2), (x)0,y0,z0) The coordinate of the imaging center S in the geocentric inertial coordinate system, and the (A, B and C) are the direction numbers of the characteristic light rays in the geocentric inertial coordinate system. The coordinates of the imaging center S can be obtained by orbit recursion (assuming that the imaging center is located at the center of mass of the satellite), and the direction number of the characteristic light needs to be obtained by coordinate change.
The coordinate conversion process of the direction number of the characteristic light rays is as follows:
in the sensor coordinate system, the direction vector of the characteristic light ray is expressed as:
Figure BDA0002001557820000091
wherein, VL、VRDividing the light into direction vectors of characteristic light rays on two sides, wherein IFOV is a camera field angle;
converting the direction vector under the sensor coordinate system into the direction vector under the geocentric inertial coordinate system as follows:
Figure BDA0002001557820000092
wherein R isTIs a whole-star attitude angle matrix,
Figure BDA0002001557820000093
Figure BDA0002001557820000094
is a roll angle, omega is a pitch angle, and kappa is a yaw angle; rREAs determined by the state of motion of the satellites,
Figure BDA0002001557820000095
wherein the content of the first and second substances,
Figure BDA0002001557820000096
Y2=Z2∧X2p (t) and v (t) are a position vector and a velocity vector at time t, respectively.
The coordinates of P1, P2, P3 and P4 in the geocentric inertial system can be solved through the joint type (1) and (2), and the coordinates of the P1, the P2, the P3 and the P4 in the geostationary coordinate system can be obtained through coordinate transformation.
Step 32, solving all the rectangular strips according to the vertex coordinates of all the rectangular strips to obtain a covering polygon;
in this embodiment, a vector polygon logical operation mode is adopted, and a coverage polygon is obtained by solving all rectangular strips according to vertex coordinates of all rectangular strips.
The unit of the vertex coordinates of the four vertexes of each rectangular strip is mm;
rounding the vertex coordinates of all the rectangular strips;
and solving all the rectangular strips according to the rounded vertex coordinates of all the rectangular strips to obtain a covering polygon.
Since the vector polygon logical operation method needs to be implemented by integer numbers, vertex coordinates of the covered polygons need to be converted into integer numbers, and the calculation precision depends on the unit of floating point numbers and rounding errors of the integer numbers converted from the floating point numbers. In order to ensure the accuracy of solving the coverage polygon, in this embodiment, the vertex coordinate unit of the coverage polygon after gaussian projection is converted into mm, and then rounding operation is performed, so as to ensure that the coordinate accuracy loss caused by the rounding operation is less than 1 mm.
Step 33, intersecting the coverage polygon with the target area to obtain an effective coverage polygon;
the rule for solving and obtaining the coverage polygon for all rectangular strips is as follows:
(1)(LC∪LS)or(LS∪LC)=LI;
(2)(RC∪RS)or(RS∪RC)=RI;
(3)(LS∪RC)or(LC∪RS)=MN;
(4)(RS∪LC)or(RC∪LS)=MX;
the rule for intersecting the coverage polygon with the target area is as follows:
the rule for Unlike edges is:
(1)(LC∩LS)or(LS∩LC)=LI;
(2)(RC∩RS)or(RS∩RC)=RI;
(3)(LS∩RC)or(LC∩RS)=MX;
(4)(RS∩LC)or(RC∩LS)=MN;
the rule for intersecting Like edges is as follows:
(5)(LC∩RC)or(RC∩LC)=LI and RI;
(6)(LS∩RS)or(RS∩LS)=LI and RI;
l, R represents the left and right boundaries of the polygon, S, C represents the subject polygon and the clip polygon respectively, LC, RC, LS and RS represent the left boundary of the clip polygon, the right boundary of the clip polygon, the left boundary of the subject polygon and the left boundary of the subject polygon respectively, MN and MX represent the local minimum value and the local maximum value respectively, and LI and RI represent the middle points of the left boundary and the right boundary of the polygon after operation respectively.
The invention adopts Vatti algorithm to realize the logic operation of complex polygon. The manner in which the left and right boundaries of the polygon are determined is as follows:
the local maximum and minimum values (vertical coordinates) of the polygon vertex coordinates redefine the polygons and divide each polygon into left and right boundary pairs. Each boundary starts at a local minimum point and ends at a local maximum point; the definition of the left and right boundaries depends on whether the boundaries are located on the left or right side of the interior of the polygon. In the polygon shown in fig. 3, the local maximum points are P6 and P2, and the local minimum points are P0 and P4, and the polygon is described by two boundary pairs:
boundary pair 1: left boundary (P0P8, P8P7, P7P6), right boundary (P0P1, P1P 2);
boundary pair 2: left boundary (P4P5, P5P6), right boundary (P4P3, P3P 2).
In the intersection calculation shown in fig. 4, A, B, C, D four intersection points are determined by the intersection rule, where point a is LI (the middle point of the left boundary of the result polygon), point B is MN (the local minimum point of the result polygon), point C is RI (the middle point of the left boundary of the result polygon), and point D is MN (the local maximum point of the result polygon).
And (3) adopting a Scan-beams scanning area (figure 3) consisting of two scanning lines scanned in the horizontal direction to the obtained polygon, and carrying out intersection and union logic operation in the area according to intersection and union rules to obtain left and right boundary pairs of the result polygon. And the scanning area is updated from the bottommost end of the polygon along the direction from the vertical direction to the top end, and the left and right boundary pairs of the polygon after the logic operation are continuously updated, so that a new polygon is finally obtained.
And step 34, dividing the area of the effective coverage polygon by the area of the target area to obtain the coverage rate.
The area of the effective coverage polygon is obtained according to the following steps:
dividing the effective coverage polygon into a plurality of triangles, solving the area of each triangle one by adopting a vector product method, and summing the areas of all the triangles to obtain the area of the effective coverage polygon; wherein, the calculation formula of the area of the effective coverage polygon is as follows:
Figure BDA0002001557820000111
wherein x isi、yiIs the ithAnd the plane coordinates of the vertexes, q is the number of the vertexes of the effective coverage polygon, and S is the area of the effective coverage polygon, wherein i and q are positive integers.
Step S4, taking a yaw angle combination corresponding to the maximum coverage rate in the set F as an optimal yaw angle combination;
the roll angle optimization model for obtaining the final optimal roll angle combination in step S9 is as follows:
Figure BDA0002001557820000121
wherein, P ═ { x1, x2, …, xi, … xn }, xi is equal to [ li, ui }, and xi is equal to [ li, ui ],]in the formula xiFor a yaw angle combination, Scov (x)i) To cover the area of the polygon effectively, SobjIs the area of the target region, n is the number of elements in the set P, li、uiAnd i and n are positive integers, which are the upper and lower limits of the yaw angle of the ith track arc section.
As shown in fig. 7 and 8, the present invention decomposes the process of satellite imaging into a plurality of orbital arc imaging processes, and considers the variation of imaging width under different side swing angles, a side swing angle is randomly selected on each orbital arc, and different side swing angles form different rectangular strips, and the widths of different rectangular strips are different (as visible from left to right in fig. 7, the width of a first rectangular strip is different from that of a second rectangular strip), so that the plurality of rectangular strips with different widths in the present application is different from the prior art which only uses a single type of strip with fixed width. The method overcomes the defect that the traditional method based on the division of the strips with fixed width does not take the change of the imaging width under different yaw angles into consideration, further adopts a PSODE algorithm to continuously update, cross and mutate a plurality of yaw angle combinations to generate new yaw angle combinations, obtains the yaw angle combinations meeting the imaging conditions from the new yaw angle combinations, and improves the accuracy of determining the optimal yaw angle combinations according to the final optimal yaw angle combinations meeting the imaging conditions.
Particle Swarm Optimization (PSO) and Differential Evolution (DE) are both population-based intelligent algorithms, and a Particle Swarm-Differential Evolution (PSODE) algorithm is provided by combining the advantages of the two methods.
Performing multiple iterations by using a PSODE algorithm to determine a final optimal yaw angle combination, wherein the steps specifically comprise:
step S5, updating, mutating and crossing each element in the set P to obtain a set P' consisting of corresponding new elements;
step S6, obtaining a set F 'according to the coverage rate corresponding to each element in the set P';
step S7, comparing each element of the set F 'with the corresponding element in the set F, and if each element of the set F' is larger than the corresponding element in the set F, replacing the corresponding element in the set F with the element of the set F ', and replacing the corresponding element in the set P with the element of the set P';
step S8, comparing the largest element in the set F obtained in step S7 with the coverage rate corresponding to the optimal roll angle combination in step S4, and if the largest element in the set F obtained in step S7 is greater than the coverage rate corresponding to the optimal roll angle combination in step S4, taking the element in the set P corresponding to the largest element in the set F obtained in step S7 as the optimal roll angle combination to be optimized;
and step S9, replacing the set P and the set F obtained in the step S5 with the set P and the set F obtained in the step S7, and repeatedly executing the steps S5 to S8 for K times to obtain the final optimal yaw angle combination, wherein K is an integer greater than or equal to 1.
TABLE 1 satellite orbits and sensor parameters
Figure BDA0002001557820000131
An example is proposed, wherein a Chinese area is taken as a target area, the orbit of a satellite and related parameters of a sensor are shown in table 1, and the satellite completes the facing area imaging task through multiple transit. The simulation starting time is coordinated Universal Time (UTC) in 2017, 5, month, 1, 4:00:00, and the terminating time is UTC in 2017, 5, month, 10, 4:00: 00.
Step a1, randomly selecting 1 yaw angle between the upper limit and the lower limit of the yaw angle of each orbit arc segment aiming at 30 orbit arc segments of the satellite; the respective selected yaw angles for the 30 orbital arc segments are shown in table 2.
Step a2, repeating step Sa1 n times to obtain n yaw angle combinations, and recording the n yaw angle combinations as a set P; because of the large number, all combinations of yaw angles are not listed here.
Step a3, determining the coverage rate corresponding to each element in the set P, and recording the coverage rate corresponding to each element in the set P as a set F; namely, determining the coverage rate corresponding to each of the plurality of individuals (elements in the set P), wherein the coverage rate obtained according to each individual is called the individual fitness value of the individual;
TABLE 2 yaw angle of each arc segment
Numbering Side swing angle (°) Numbering Side swing angle (°)
x1 1 x16 2
x2 2 x17 6
x3 5 x18 4
x4 6 x19 5
x5 4 x20 4
x6 7 x21 5
x7 1 x22 6
x8 2 x23 4
x9 3 x24 1
x10 5 x25 5
x11 6 x26 6
x12 4 x27 4
x13 1 x28 2
x14 5 x29 4
x15 1 x30 2
Step a31, calculating the vertex coordinates of the rectangular strip for the crossed track arc with the number x1 in the time period from 14:24:33 in 5/1/2017 to 14:27:31 in 5/1/2017.
Taking the calculation method of the vertex coordinates of the point P1 as an example:
first, t ═ 2101.603, a ═ 0.624, B ═ 0.348, and C ═ 0.727 were determined.
The following equations are simultaneously set forth:
Figure BDA0002001557820000151
Figure BDA0002001557820000152
the P1 point Cartesian coordinates (three-dimensional coordinates x, y and z) are (-4893.443, 4136.054 and 6311.774) and the unit is meter, and the coordinates after conversion into longitude and latitude are (133.474 and 56.2042).
Similarly, the coordinates P1, P2, P3, and P4 corresponding to the four corner points of the covering strip are sequentially obtained as follows: (133.474, 56.2042),(130.29, 55.7156),(135.211, 44.0479),(137.706, 44.4612).
Step a32, calculating a coverage rate for 30 arcs from 1/5/2017 to 10/5/2017 by using the vector polygon logical operation method provided herein.
Step a32, obtaining the vertex coordinates of all rectangular strips according to the yaw angle of the table, and finding out the total area of the rectangular strips corresponding to 30 orbit arc segments as: 9.6008*1018mm2
And then, performing parallel operation and intersection on 30 rectangular strips to obtain the effective coverage area of the polygon as follows: 5.1136*1018mm2
As shown in fig. 7, the area of the effective coverage polygon was divided by the area of the chinese region to obtain a coverage rate of 53.2625%.
And obtaining the individual adaptation values of the rest individuals according to the steps.
Step a4, using the individual with the largest individual fitness value corresponding to each individual in the set P as the overall extreme value of the initial generation population, where the overall extreme value refers to the yaw angle combination with the largest coverage rate obtained correspondingly from all the yaw angle combinations.
A5, updating, mutating and crossing each element (namely individual) in the set P to obtain a set P' consisting of corresponding new elements;
step a6, obtaining a set F 'according to the set P';
step a7, comparing each element of the set F 'with the corresponding element in the set F, if each element of the set F' is larger than the corresponding element in the set F, replacing the corresponding element in the set F with the element of the set F ', and replacing the corresponding element in the set P with the element of the set P';
for example, the position and the speed of the first individual in the set P are updated, varied and crossed according to a PSODE algorithm to obtain a first individual in the set P ', the individual fitness value obtained by the first individual in the set P is compared with the individual fitness value of the first individual in the set P ', and if the individual fitness value obtained by the first individual in the set P is greater than or equal to the individual fitness value of the first individual in the set P ', the first individual in the set P is retained; and if the individual fitness value obtained by the first individual in the set P is smaller than that of the first individual in the set P ', taking the first individual in the set P' as the first individual in the set P.
TABLE 3 optimization results of the yaw angles of the arc segments
Numbering Side swing angle (°) Numbering Side swing angle (°)
x1 18.2035 x16 -31.3346
x2 -32.9375 x17 18.5529
x3 -22.6816 x18 -0.4457
x4 30.4829 x19 -6.3483
x5 28.7086 x20 -22.5342
x6 24.4881 x21 -20.2448
x7 -12.7036 x22 -33.4622
x8 30.8017 x23 -34.6931
x9 11.214 x24 -10.7434
x10 15.4382 x25 -20.5764
x11 -33.0486 x26 11.0523
x12 -14.7236 x27 4.7268
x13 -34.6931 x28 27.587
x14 34.0358 x29 -16.2205
x15 2.0449 x30 -19.6252
A step a8, comparing the largest element in the set F obtained in the step a7 with the coverage rate corresponding to the whole extreme value of the initial generation population in the step a4, and if the largest element in the set F obtained in the step a7 is larger than the coverage rate corresponding to the whole extreme value of the initial generation population in the step a4, taking the element in the set P corresponding to the largest element in the set F obtained in the step 7 as the second generation whole extreme value;
step a9, replacing the set P and the set F in step a5 with the set P and the set F obtained in step a7, and repeating 500 times steps a5 to a8 to obtain the final overall extremum, i.e. the final overall extremum is the optimal roll angle combination, such as the roll angle combinations shown in table 3. As shown in fig. 8, the coverage map is obtained by combining the optimal roll angles after 500 iterations.
Because the PSODE algorithm is adopted for iteration process, more data are provided, and the optimal yaw angle combination selected by each generation and the corresponding change of the coverage rate are better embodied by using the graph shown in FIG. 6. In fig. 6, the abscissa is an evolution algebra, and the ordinate is a coverage rate obtained by the corresponding algebra, and it can be seen from fig. 6 that the coverage rate gradually increases with the increase of the number of iterations.
The coverage was 93.8754% according to the roll angle of table 3. 93.8754% the coverage increased by nearly 1-fold compared to 53.2625%.
The specific embodiments described herein are merely illustrative of the spirit of the invention. Various modifications or additions may be made to the described embodiments or alternatives may be employed by those skilled in the art without departing from the spirit or ambit of the invention as defined in the appended claims.
For the system embodiment, since it is basically similar to the method embodiment, the description is simple, and for the relevant points, refer to the partial description of the method embodiment.
The embodiments in the present specification are described in a progressive manner, each embodiment focuses on differences from other embodiments, and the same and similar parts among the embodiments are referred to each other.
The area coverage-oriented multi-track splicing imaging optimization method provided by the application is described in detail, a specific example is applied in the method to explain the principle and the implementation mode of the application, and the description of the embodiment is only used for helping to understand the method and the core idea of the application; meanwhile, for a person skilled in the art, according to the idea of the present application, there may be variations in the specific embodiments and the application scope, and in summary, the content of the present specification should not be construed as a limitation to the present application.

Claims (10)

1. An area coverage-oriented multi-rail splicing imaging optimization method is characterized by comprising the following steps:
step S1, randomly selecting 1 yaw angle between the upper limit and the lower limit of the yaw angle of each orbit arc segment aiming at m orbit arc segments of the satellite to form a yaw angle combination; wherein m is a positive integer;
step S2, repeating step S1 n times to obtain n yaw angle combinations, and recording the n yaw angle combinations as a set P, wherein n is a positive integer;
step S3, determining the coverage rate corresponding to each element in the set P, and recording the coverage rate corresponding to each element in the set P as a set F;
step S4, taking a yaw angle combination corresponding to the maximum coverage rate in the set F as an optimal yaw angle combination;
step S5, updating, mutating and crossing the speed and position of each element in the set P by adopting a PSODE algorithm to obtain a set P' consisting of corresponding new elements;
step S6, obtaining a set F 'according to the coverage rate corresponding to each element in the set P';
step S7, comparing each element of the set F 'with the corresponding element in the set F, and if each element of the set F' is larger than the corresponding element in the set F, replacing the corresponding element in the set F with the element of the set F ', and replacing the corresponding element in the set P with the element of the set P';
step S8, comparing the largest element in the set F obtained in step S7 with the coverage rate corresponding to the optimal roll angle combination in step S4, and if the largest element in the set F obtained in step S7 is greater than the coverage rate corresponding to the optimal roll angle combination in step S4, taking the element in the set P corresponding to the largest element in the set F obtained in step S7 as the optimal roll angle combination to be optimized;
step S9, replacing the set P and the set F in step S5 with the set P and the set F obtained in step S7, and repeating steps S5 to S8K times to obtain a final optimal yaw angle combination, where K is an evolutionary algebra and K is an integer greater than or equal to 1;
the PSODE algorithm is adopted to continuously update, cross and mutate a plurality of side sway angle combinations to generate new side sway angle combinations, the side sway angle combinations meeting the imaging conditions are obtained from the new side sway angle combinations, and the accuracy of determining the optimal side sway angle combinations is improved according to the final optimal side sway angle combinations meeting the imaging conditions.
2. The method of claim 1, wherein determining the coverage corresponding to each element in the set P comprises:
determining a rectangular strip corresponding to each yaw angle according to each yaw angle in a yaw angle combination, a central projection imaging geometric relation and the position and the speed of a satellite, and obtaining vertex coordinates of four vertexes of each rectangular strip;
solving all the rectangular strips according to the vertex coordinates of all the rectangular strips to obtain a covering polygon;
intersecting the coverage polygon with a target area to obtain an effective coverage polygon;
and dividing the area of the effective coverage polygon by the area of the target area to obtain the coverage rate.
3. The method of claim 2, wherein the rectangular strips are formed on the ground surface through four intersections of feature rays from the imaging center S corresponding to each yaw angle to four corner points of the area array sensor and the ground surface.
4. The method of claim 2, wherein the rectangular strips are formed by scanning feature rays from the imaging center S corresponding to each yaw angle to two corner points of the linear array sensor over a period of time on the ground surface.
5. The method of claim 2, wherein the vertex coordinates of each vertex of each rectangular strip are determined according to the following steps:
in an inertial coordinate system, the earth is taken as a standard rotation ellipsoid, and the equation of the standard rotation ellipsoid is a first formula:
Figure FDA0003269326550000021
wherein, (x, y, z) is a point on the surface of the standard ellipsoid of revolution, and a and b are respectively the length of the major semi-axis and the minor semi-axis of the standard ellipsoid of revolution;
the equation of the characteristic light ray where the vertex is located in the geocentric inertial coordinate system is a second formula:
Figure FDA0003269326550000022
wherein (x)0,y0,z0) The coordinate of the imaging center S under the geocentric inertial coordinate system is shown, and the (A, B and C) are the direction numbers of the characteristic light rays where the vertex is located under the geocentric inertial coordinate system;
and combining the first formula and the second formula to obtain the vertex coordinates of each vertex of each rectangular strip.
6. The method of claim 2, wherein the area of the effective coverage polygon is obtained by:
dividing the effective coverage polygon into a plurality of triangles, solving the area of each triangle one by adopting a vector product method, and summing the areas of all the triangles to obtain the area of the effective coverage polygon; wherein, the calculation formula of the area of the effective coverage polygon is as follows:
Figure FDA0003269326550000031
wherein x isi、yiIs the plane coordinate of the ith vertex, q is the number of the vertices of the effective coverage polygon, and S is the area of the effective coverage polygon, wherein i and q are positive integers.
7. The method of claim 2, wherein the vertex coordinates of the four vertices of each rectangular strip are in mm;
rounding the vertex coordinates of all the rectangular strips;
and solving all the rectangular strips according to the rounded vertex coordinates of all the rectangular strips to obtain a covering polygon.
8. The method of claim 2, wherein the intersection rules of intersection include Unlike intersection rules and like intersection rules.
9. The method of claim 2, wherein the upper and lower limits of the yaw angle for each orbital arc are determined according to the following steps:
obtaining the track distance d between each vertex of the rectangular strip and the satellite point according to the position relation between the satellite point track of each orbit arc segment of the satellite and the rectangular strip;
determining a yaw angle beta required when a satellite main optical axis points to each vertex for imaging according to the satellite orbit height h and the camera field angle IFOV:
β=atan(d/h)+IFOV
sequencing the sidesway angles beta corresponding to the vertexes, wherein the maximum value of the sidesway angle is the maximum sidesway angle u of the track arc section, and the minimum value of the sidesway angle is the minimum sidesway angle l of the track arc section;
let the range of the satellite yaw angle be (-Roll)max,Rollmax);
If l is<-RollmaxAnd u>RollmaxThe upper and lower limits of the yaw angle of each track arc segment are respectively-Rollmax,Rollmax
If l is<-RollmaxAnd u < RollmaxThe upper and lower limits of the yaw angle of each track arc segment are respectively-Rollmax,u;
If l > -RollmaxAnd u < RollmaxThe upper and lower limits of the yaw angle of each track arc section are l and u respectively;
if l > -RollmaxAnd u>RollmaxThe upper and lower limits of the yaw angle of each track arc segment are l, Rollmax
10. The method according to claim 1, wherein the roll angle optimization model for obtaining the final optimal roll angle combination in step S9 is as follows:
Figure FDA0003269326550000041
wherein, P ═ { x ═ x1,x2,…,xi,…xn},xi∈[li,ui]In the formula xiFor a yaw angle combination, Scov (x)i) To cover the area of the polygon effectively, SobjIs the area of the target region, n is the number of elements in the set P, li、uiAnd i and n are positive integers, which are the upper and lower limits of the yaw angle of the ith track arc section.
CN201910214391.XA 2019-03-20 2019-03-20 Multi-rail splicing imaging optimization method for area coverage Active CN109993719B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910214391.XA CN109993719B (en) 2019-03-20 2019-03-20 Multi-rail splicing imaging optimization method for area coverage

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910214391.XA CN109993719B (en) 2019-03-20 2019-03-20 Multi-rail splicing imaging optimization method for area coverage

Publications (2)

Publication Number Publication Date
CN109993719A CN109993719A (en) 2019-07-09
CN109993719B true CN109993719B (en) 2022-01-18

Family

ID=67130723

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910214391.XA Active CN109993719B (en) 2019-03-20 2019-03-20 Multi-rail splicing imaging optimization method for area coverage

Country Status (1)

Country Link
CN (1) CN109993719B (en)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111693994B (en) * 2020-04-30 2022-06-03 中国科学院空天信息创新研究院 Airborne synthetic aperture radar route laying method, device, equipment and storage medium
CN111950877B (en) * 2020-07-31 2023-11-14 上海卫星工程研究所 Multi-star formation collaborative region imaging autonomous task planning method and system

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102479289A (en) * 2010-11-30 2012-05-30 中国人民解放军国防科学技术大学 Regional division method for satellite observation
CN103294925A (en) * 2013-06-20 2013-09-11 中国科学院遥感与数字地球研究所 Method for realizing multi-satellite combined imaging
CN109063938A (en) * 2018-10-30 2018-12-21 浙江工商大学 Air Quality Forecast method based on PSODE-BP neural network

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101614813B (en) * 2009-07-23 2011-08-17 航天东方红卫星有限公司 Revisiting orbit determining method of all-weather coverage satellite
CN108446254A (en) * 2018-02-08 2018-08-24 中国科学院遥感与数字地球研究所 Collecting method and device
CN108845976B (en) * 2018-06-25 2019-05-10 湖南国科轩宇信息科技有限公司 Large-scale area observation scheduling method and system under multi satellites joint imaging

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102479289A (en) * 2010-11-30 2012-05-30 中国人民解放军国防科学技术大学 Regional division method for satellite observation
CN103294925A (en) * 2013-06-20 2013-09-11 中国科学院遥感与数字地球研究所 Method for realizing multi-satellite combined imaging
CN109063938A (en) * 2018-10-30 2018-12-21 浙江工商大学 Air Quality Forecast method based on PSODE-BP neural network

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
一种基于改进PSO算法的高时间分辨率遥感卫星星座优化设计方法;沈欣等;《武汉大学学报》;20181231;第43卷(第12期);全文 *
一种改进的差分进化算法及其应用;夏莘媛等;《河北工业大学学报》;20150228;第44卷(第1期);全文 *
光学遥感卫星轨道设计若干关键技术研究;沈欣;《中国博士学位论文全文数据库 工程科技Ⅱ辑》;20130615;全文 *
面向区域成像任务的环月卫星侧摆角优化方法;李仕学等;《武汉大学学报》;20180614;摘要部分、第1-4部分,图1-9,表1-4 *

Also Published As

Publication number Publication date
CN109993719A (en) 2019-07-09

Similar Documents

Publication Publication Date Title
CN109828607B (en) Unmanned aerial vehicle path planning method and system for irregular obstacles
CN103869820B (en) A kind of rover ground navigation planning control method
CN102243074B (en) Method for simulating geometric distortion of aerial remote sensing image based on ray tracing technology
CN106780712B (en) Three-dimensional point cloud generation method combining laser scanning and image matching
US8610708B2 (en) Method and apparatus for three-dimensional image reconstruction
CN109993719B (en) Multi-rail splicing imaging optimization method for area coverage
CN106597434B (en) It is a kind of based on pushing away the quick Satellite Targets decomposition method and system for sweeping track
CN111949922B (en) Method and system suitable for on-board rapid calculation of multi-time window of ground detection task
CN111666661B (en) Method and system for planning imaging multi-strip splicing task in single track of agile satellite
CN102346922A (en) Space remote sensing load imaging geometric distortion three-dimensional visualization simulation method
Holenstein et al. Watertight surface reconstruction of caves from 3D laser data
CN115564926B (en) Three-dimensional patch model construction method based on image building structure learning
Majeed et al. A multi-objective coverage path planning algorithm for UAVs to cover spatially distributed regions in urban environments
CN111336994A (en) Remote sensing satellite coverage analysis method based on combination of graphics and numerical calculation
CN110866015B (en) Moving target moving range recording method based on local grid
CN104764443A (en) Optical remote sensing satellite rigorous imaging geometrical model building method
Stucker et al. ResDepth: Learned residual stereo reconstruction
CN102519436A (en) Chang&#39;e-1 (CE-1) stereo camera and laser altimeter data combined adjustment method
CN114721436A (en) Automatic air route planning method for unmanned aerial vehicle-mounted hyperspectral imaging system
CN112149911A (en) Hypersensitive short satellite same-orbit multipoint target in-motion imaging task planning method
CN103778610B (en) A kind of spaceborne line array sensor hangs down the geometry preprocess method of rail sweeping image
KR100965838B1 (en) An Implicit Geometric Regularization of Building Polygon Using LiDAR Data
Shang et al. Topology-based UAV path planning for multi-view stereo 3D reconstruction of complex structures
CN111856747B (en) Optimal transmission-based reflector design method
Rau et al. Lod generation for 3d polyhedral building model

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant