CN109993719B - Multi-rail splicing imaging optimization method for area coverage - Google Patents
Multi-rail splicing imaging optimization method for area coverage Download PDFInfo
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 57
- 238000003384 imaging method Methods 0.000 title claims abstract description 45
- 238000005457 optimization Methods 0.000 title claims abstract description 21
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 15
- 238000004364 calculation method Methods 0.000 claims abstract description 14
- 239000013598 vector Substances 0.000 claims description 15
- 230000003287 optical effect Effects 0.000 claims description 3
- 238000012163 sequencing technique Methods 0.000 claims description 3
- 230000007547 defect Effects 0.000 abstract description 5
- 230000008569 process Effects 0.000 description 8
- 230000008859 change Effects 0.000 description 6
- 238000010586 diagram Methods 0.000 description 6
- 230000000694 effects Effects 0.000 description 3
- 238000006243 chemical reaction Methods 0.000 description 2
- 239000002245 particle Substances 0.000 description 2
- 230000006978 adaptation Effects 0.000 description 1
- 238000007792 addition Methods 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 238000010845 search algorithm Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20212—Image combination
- G06T2207/20221—Image 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
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:
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:
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:
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:
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:
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:
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:
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:
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:
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:
wherein R isTIs a whole-star attitude angle matrix,
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:
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:
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
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 | |
2 |
| |
2 | x17 | 6 |
| |
5 | x18 | 4 |
| x4 | 6 | |
5 |
| x5 | 4 | x20 | 4 |
| x6 | 7 | |
5 |
| x7 | 1 | x22 | 6 |
| |
2 | x23 | 4 |
| x9 | 3 | x24 | 1 |
| |
5 | |
5 |
| x11 | 6 | x26 | 6 |
| x12 | 4 | x27 | 4 |
| x13 | 1 | |
2 |
| |
5 | x29 | 4 |
| x15 | 1 | |
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:
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:
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:
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:
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:
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.
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)
| 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)
| 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)
| 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 |
-
2019
- 2019-03-20 CN CN201910214391.XA patent/CN109993719B/en active Active
Patent Citations (3)
| 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)
| 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 | |
| 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 | |
| CN102496181A (en) | True-orthophotomap making method oriented to large-scale production | |
| Stucker et al. | ResDepth: Learned residual stereo reconstruction | |
| CN102519436A (en) | Chang'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 | |
| CN113566793A (en) | True orthoimage generation method and device based on unmanned aerial vehicle oblique image | |
| 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 |
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 |