A method for extracting the boundary of the hazardous area caused by an asteroid hitting the Earth

Abstract

The present invention discloses a method for extracting the boundaries of a hazard zone in the event of an asteroid impact on Earth, relating to the field of planetary defense. The method comprises the following steps: S1. obtaining the longitude, latitude, and hazard radius data for all airburst and / or impact points; S2. evenly dividing the longitude and latitude into multiple intervals according to a preset resolution to obtain corresponding map grids, and setting the initial values of the map grid state variables to 0; S3. calculating the state variables in all map grids; S4. calculating the preliminary boundaries of the hazard zone; and S5. calculating the final boundaries of the hazard zone. The present invention provides a method for extracting the boundaries of a hazard zone in the event of an asteroid impact on Earth. The method uses a large number of airburst / impact points and hazard radii obtained from an asteroid impact hazard assessment engineering model as input, first extracting the preliminary boundaries of the hazard zone using the map grid state variables, and then obtaining the final boundaries using a concave hull algorithm. In this process, parallel computing is used to handle the massive computational effort resulting from the massive input.

Description

Technical Field

[0001] The present invention relates to the field of planetary defense in the direction of asteroid impact hazard assessment, and in particular to a method for extracting the boundary or envelope of an asteroid impact hazard area. Background Art

[0002] Asteroid impacts pose a real threat to human survival and development. Planetary defense has been a hot research area internationally for nearly 20 years and represents a major security need in my country. Asteroid impact hazard assessment is a crucial component and research direction of planetary defense, providing the basis for defense decision-making and ground-based civil defense. Extracting the boundaries of the asteroid impact hazard zone helps visualize the hazard range in geographic information systems, thereby better supporting planetary defense.

[0003] The hazards of asteroids hitting the Earth are characterized by low probability, high hazard, and randomness, so the hazard area needs to include any possible hazard point. In hazard assessment, engineering models are usually used to obtain any airburst / impact point and hazard radius within the input parameter range of impact position and kinematic parameters, asteroid characteristics, entry into the atmosphere and impact surface model. Each airburst / impact point and hazard radius determines the local hazard area, and the hazard assessment will provide a huge number (tens of millions) of airburst / impact point and hazard radius combinations. For example Figure 10 As shown ( Figure 10 A in the middle is the boundary of the hazard area. P i For the i An airburst\impact point, R i For the i Airblast\impact point hazard radius R, P bj For the j The boundary A of the hazard area is the n massive airbursts\impact points mentioned above. P i ( i =1, 2…n) and hazard radius R i ( i =1, 2…n) combination determines the boundary of the global area. However, in practical applications, there are two difficulties in extracting the boundary of the hazard area:

[0004] First, the number of airburst / impact points is huge, and shortening the calculation time is difficult;

[0005] Second, the boundary of the hazard area has a complex shape, and the extracted boundary needs to have high fidelity.

[0006] However, there is currently no work in China to extract the boundaries of the hazardous areas caused by asteroid impacts on the Earth. Summary of the Invention

[0007] An object of the present invention is to solve at least the above problems and / or disadvantages and to provide at least the advantages which will be described hereinafter.

[0008] In order to achieve these objectives and other advantages of the present invention, a method for extracting the boundary of an asteroid impact hazard zone is provided, comprising:

[0009] S1. Obtain the longitude, latitude, and hazard radius data of all airburst and / or impact points;

[0010] S2. Divide the longitude and latitude into multiple intervals evenly according to the preset resolution to obtain corresponding map grids, and set the initial values of the map grid state variables to 0;

[0011] S3. For any input data, calculate the index values of the longitude interval and latitude interval in the map grid to obtain the corresponding local grid interval;

[0012] In the local grid interval, it is judged whether each grid center falls within the hazard radius of the input data. If it does, the corresponding grid state variable is set to 1;

[0013] S4. Based on the judgment of whether the state variables of adjacent grids are equal, determining whether the corresponding grid is a preliminary boundary value of the hazard area, traversing the longitude grid lines to obtain a corresponding preliminary boundary data set;

[0014] S5. Taking the preliminary boundary dataset as input, the concave hull algorithm is used to obtain the final boundary dataset of the hazard area, thereby completing the boundary extraction of the hazard area.

[0015] Preferably, in S2, the multiple intervals refer to evenly dividing the longitude within the full range of [-180°, 180°] into ni-1 intervals, and evenly dividing the latitude within the full range of [-90°, 90°] into nj-1 intervals, where ni-1 is the number of longitude intervals and nj-1 is the number of latitude intervals;

[0016] The state variables of each grid are mgv i,j To characterize.

[0017] Preferably, in S3, the index value is obtained in the following manner:

[0018] S30, for i Input data, based on the corresponding airburst and / or impact point P i Longitude L i ,latitude B iand hazard radius R i , calculate the corresponding longitude interval [ L i - dL , L i + dL ]、latitude range[ B i - dB , B i + dB ];

[0019] in, dL, dB They are characterized by the following formulas:

[0020] ,

[0021] In the above formula, dL is half the maximum length of the longitude interval, dB is half of the maximum length of the latitude interval, r e is the average radius of the Earth;

[0022] S31. Based on the longitude interval and latitude interval calculated in S30, calculate the index values corresponding to the longitude interval and latitude interval in the map grid.

[0023] Preferably, in S5, after taking the preliminary boundary dataset as input, the method further includes:

[0024] S50, dividing the input data into m parts, performing the concave hull algorithm on each part of the data to obtain a sub-hazard boundary;

[0025] S51, gather all sub-hazard boundaries and execute the concave hull algorithm again to obtain the final boundary data set;

[0026] Each data point in the data set represents the longitude and latitude information of the corresponding boundary point, and all data points are arranged in an orderly manner in a counterclockwise or clockwise direction on the map so that they can be directly input into the geographic information system for boundary drawing in the later stage.

[0027] The present invention has at least the following beneficial effects: the present invention can be regarded as a post-processing process of the asteroid impact hazard assessment engineering calculation, providing support for illustrating the hazard range in the geographic information system. In actual implementation, the massive airburst / impact points and hazard radius obtained from the asteroid impact hazard assessment engineering model are used as input, the preliminary boundary of the hazard area is first extracted using the map grid state variables, and then the final boundary is obtained using the concave hull algorithm. In this process, parallel computing is used to handle the massive calculation caused by the massive input.

[0028] Other advantages, objectives and features of the present invention will be reflected in part from the following description and will be understood by those skilled in the art through study and practice of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0029] Figure 1 A flow chart of the method for extracting the boundary of the hazard zone of an asteroid impact on the Earth;

[0030] Figure 2 This is a flow chart of the map grid state variable calculation method based on openMP parallelism;

[0031] Figure 3 Flowchart of the method for calculating preliminary boundaries of the hazard area;

[0032] Figure 4 This is a flow chart of the method for calculating the final boundary of the hazard area based on region decomposition and concave hull algorithm;

[0033] Figure 5 This is an example of a map grid variable cloud map;

[0034] Figure 6 This is an example of the preliminary boundary of the hazard area;

[0035] Figure 7 for Figure 6 A partial enlarged view of

[0036] Figure 8 This is an example diagram of the final boundary of the hazard area;

[0037] Figure 9 for Figure 8 A partial enlarged view of

[0038] Figure 10 Schematic diagram of the boundaries of the area where an asteroid could hit the Earth. DETAILED DESCRIPTION

[0039] The present invention will be described in further detail below in conjunction with the accompanying drawings so that those skilled in the art can implement the invention with reference to the description.

[0040] like Figure 1 As shown, the present invention provides a method for extracting the boundary of an asteroid impact hazard zone on Earth, comprising:

[0041] S1, read in the longitude, latitude and hazard radius data of all airburst / impact points (a total of n pieces of data), i Data mid-air explosion\impact point P i The longitude, latitude and hazard radius are respectively denoted as L i 、B i 、 R i .

[0042] S2. According to the preset resolution, the full range of longitude [-180°, 180°] is evenly divided into ni-1 intervals, and the full range of latitude [-90°, 90°] is evenly divided into nj-1 intervals, and the map grid state variables are set. mgv i,j =0( i =1、2…ni, j =1, 2…nj).

[0043] In S2, the map grid is divided as follows:

[0044] S20, set the map grid length resolution to d , based on the average radius of the Earth r e Get the angular resolution of the grid δ ,and δ It is characterized by the following formula:

[0045]

[0046] S21, according to angular resolution δ, The entire longitude range [-180°, 180°] is evenly divided into ni-1 intervals, and the entire latitude range [-90°, 90°] is evenly divided into nj-1 intervals. The calculation formulas for ni and nj are:

[0047]

[0048]

[0049] In the above formula, int() is the rounding function;

[0050] S22. Calculate the grid point data separately, that is, calculate the longitude discrete value separately by the following formula Lc i and latitude Bc j :

[0051] ,

[0052] ,

[0053] S23. Set the initial value of the state variable mgv i,j =0( i =1、2…ni,j =1, 2…nj).

[0054] S3. Calculate the state variables in all map grids mgv i,j ( i =1、2…ni, j =1, 2…nj), if the grid ( i,j ) is within the local hazard radius, then mgv i,j This step, on the one hand, only performs calculations within the local grid domain of each input airburst / impact point, avoiding ineffective calculations of searching the global grid; on the other hand, it uses multi-core parallel computing to shorten the calculation time;

[0055] like Figure 2 As shown, get the map grid state variable mgv i, The steps are as follows:

[0056] S30, for i Input data, i.e. airburst and / or impact point P i Longitude L i ,latitude B i and hazard radius R i , calculate the local longitude interval [ L i - dL , L i + dL ], and the local latitude range [ B i - dB , B i + dB ],in, dL 、 dB It is characterized by the following formula:

[0057] ,

[0058] Then use the following formula to calculate the index of the longitude interval in the map grid [ i1 , i2 ], and the index of the latitude interval in the map grid [ j1 , j2 ].

[0059] ,

[0060] ,

[0061] S31, within the local map grid formed by the local longitude interval and latitude interval, that is, for each grid ( iL , jL )( i ∈[ i1 , i2 ]、 j ∈[ j1 , j2 ]), computational grid ( iL , jL ) and the airburst and / or impact point P i distance D ,like D ≤ R i , then the grid center is within the hazard range, then let mgv iL,jL = 1. Distance D The calculation formula is:

[0062]

[0063]

[0064]

[0065]

[0066]

[0067]

[0068]

[0069] in,( x , y , z ) is the grid ( iL , jL )'s center in the geocentric coordinate system, ( x i , y i , z i ) is the airburst and / or impact point P i The coordinates in the geocentric coordinate system are Lc iL and Bc jL That is the grid calculated in step S22 ( iL , jL )coordinate, r e is the average radius of the Earth, δ The grid angular resolution set for the S20 step.

[0070] S32. Traverse all input data and execute S30~S31. During this process, use shared memory parallel programming (OpenMP) technology to accelerate the loop. Specifically, add some parallel instructions to the algorithm, such as $ompparallel, $omp do, etc.

[0071] S4. Calculate the preliminary boundary of the hazard area. mgv i,j The value of is used to obtain the preliminary boundary of the hazard area.

[0072] like Figure 3 As shown, the steps for calculating the preliminary boundary of the hazard area are:

[0073] S40, for iL Longitude grid lines, find the boundary points. The characteristic of the boundary points is that the adjacent grid state variables are not equal, then it is considered that ( iL , jL) The center of the map grid is the preliminary hazard boundary, which is recorded.

[0074] S41, traverse all longitude grid lines ( i ∈[1,ni]), execute step S40 to record all map grids belonging to the preliminary boundary of the hazard.

[0075] S5. Calculate the final boundary of the hazard zone. Using the calculated preliminary boundary as input, the concave hull algorithm is used to obtain the final boundary. This step utilizes multi-core parallel computing to shorten computation time.

[0076] like Figure 4 As shown, the steps to calculate the final boundary of the hazard area are:

[0077] S50 , taking the preliminary boundary of the hazard area as input, dividing the input data into m parts (corresponding to m CPU cores).

[0078] S51. Execute the concave hull algorithm on each piece of data to obtain the sub-hazard boundary.

[0079] S52. Gather all sub-hazard boundaries together and execute the concave hull algorithm again to obtain the final boundary.

[0080] Example:

[0081] To more clearly illustrate the technical solution of the present invention, a method for extracting the boundary of an asteroid impact hazard zone is described through an example. This example also uses the method to extract the hazard boundary based on the calculation results of the hazard assessment engineering tabletop exercise conducted at the 6th International Academy of Astronautics (IAA) Planetary Defense Conference.

[0082] The specific implementation process is as follows Figure 1 As shown, including:

[0083] S1. Read in hazard assessment result data

[0084] These data are given in text or binary data form, and the corresponding hazard assessment result data table is obtained (as shown in Table 1):

[0085] Table 1

[0086]

[0087] Table 1 shows 9 rows of data on the hazard results, each row represents a possible hazard result. Each column of data in Table 1 represents a different variable, including the longitude of the airburst\impact point L i ,latitude B i and the hazard radius R i In this embodiment, there are 10 million rows of data.

[0088] S2. Divide the latitude and longitude grid

[0089] Set the grid length resolution d If the grid angle is 2.5 km, the grid angular resolution is 0.0225°. At this interval, ni and nj points are evenly distributed within the longitude and latitude of the earth, respectively, with ni=16013 and nj=8007.

[0090] S3, such as Figure 2 As shown, calculate the state variables in all map grids mgv i,j ( i =1、2…ni, j =1, 2…nj), we get Figure 5 The results shown, where Figure 5 middle mgv i,j The values are represented by colors (see the legend on the right), and the extent of the hazard zone (green) can be seen from the figure. The calculation process shows that on a 10-core CPU, the parallelization based on OpenMP reduces the time by 5 times compared to the serial calculation.

[0091] S4. Use Figure 3The process shown in the figure calculates the preliminary boundary of the hazard area and obtains the following Figure 6 The results shown. Figure 6 The hazard zone boundary can be clearly seen, but Figure 7 The enlarged diagram shows that these boundary points are redundant, and the data are not all arranged counterclockwise or clockwise, and cannot be used for geographic information system diagrams.

[0092] S5. Use Figure 4 The process shown in the figure calculates the final boundary of the hazard area and obtains Figure 8 The results shown in the figure show that all calculation processes in this example utilize 10 CPU cores of Intel i9-9900K. For input data of tens of millions, high-fidelity hazard boundaries can be obtained within 5 minutes. Figure 8 and Figure 6 , we can see that the final boundary of the hazard area can reflect the characteristics of the preliminary boundary, that is, it has high fidelity; at the same time, check Figure 9 , there is no redundancy in boundary data. Figure 8 The data shown are ordered data, that is, all data points are arranged in a counterclockwise or clockwise order. Figure 10 Through P bj Give an example of an ordered arrangement so that it can be directly input into a geographic information system (i.e. various navigation maps on the market) for illustration.

[0093] The above solution is only an illustration of a preferred embodiment, but is not limited thereto. When implementing the present invention, appropriate replacements and / or modifications can be made according to user needs.

[0094] Although the embodiments of the present invention have been disclosed above, they are not limited to the applications listed in the specification and exemplary embodiments. They can be applied to a variety of fields suitable for the present invention. Further modifications will be readily apparent to those skilled in the art. Therefore, the present invention is not limited to the specific details and illustrations shown and described herein without departing from the general concept defined by the claims and their equivalents.

Claims

1. A method for extracting the boundary of an asteroid impact hazard zone, characterized in that: include: S1. Obtain the longitude, latitude, and hazard radius data of all airburst and / or impact points; S2. Divide the longitude and latitude into multiple intervals evenly according to the preset resolution to obtain corresponding map grids, and set the initial values of the map grid state variables to 0; S3. For any input data, calculate the index values of the longitude interval and latitude interval in the map grid to obtain the corresponding local grid interval; In the local grid interval, it is judged whether each grid center falls within the hazard radius of the input data. If it does, the corresponding grid state variable is set to 1; S4. Based on the judgment of whether the state variables of adjacent grids are equal, determining whether the corresponding grid is a preliminary boundary value of the hazard area, traversing the longitude grid lines to obtain a corresponding preliminary boundary data set; S5. Taking the preliminary boundary dataset as input, the concave hull algorithm is used to obtain the final boundary dataset of the hazard area, thereby completing the boundary extraction of the hazard area.

2. The method for extracting the boundary of the asteroid impact hazard zone according to claim 1, wherein: In S2, the multiple intervals refer to evenly dividing the longitude within the full range of [-180°, 180°] into ni-1 intervals, and evenly dividing the latitude within the full range of [-90°, 90°] into nj-1 intervals, where ni-1 is the number of longitude intervals and nj-1 is the number of latitude intervals; The state variables of each grid are mgv i,j To characterize.

3. The method for extracting the boundary of the asteroid impact hazard zone according to claim 1, wherein: In S3, the index value is obtained as follows: S30, for i Input data, based on the corresponding airburst and / or impact point P i Longitude L i ,latitude B i and hazard radius R i , calculate the corresponding longitude interval [ L i - dL , L i +dL ]、latitude range[ B i - dB , B i +dB ]; in, dL, dB They are characterized by the following formulas: ， In the above formula, dL is half the maximum length of the longitude interval, dB is half the maximum length of the latitude interval, r e is the average radius of the Earth; S31 . Based on the longitude interval and latitude interval calculated in S30 , calculate the corresponding index values of the longitude interval and latitude interval in the map grid.

4. The method for extracting the boundary of the asteroid impact hazard zone according to claim 1, wherein: In S5, after taking the preliminary boundary dataset as input, it also includes: S50, dividing the input data into m parts, performing the concave hull algorithm on each part of the data to obtain a sub-hazard boundary; S51, gather all sub-hazard boundaries and execute the concave hull algorithm again to obtain the final boundary data set; Each data point in the data set represents the longitude and latitude information of the corresponding boundary point, and all data points are arranged in an orderly manner in a counterclockwise or clockwise direction on the map so that they can be directly input into the geographic information system for boundary drawing in the later stage.

Images