Method for extracting boundary of hazard area of impact of asteroid on earth

Abstract

The invention discloses an asteroid impact earth hazard area boundary extraction method, and relates to the field of planet defense, and the method comprises the steps: S1, obtaining the longitude, latitude and hazard radius data of all air explosion and / or impact points; s2, according to a preset resolution ratio, uniformly dividing the longitude and the latitude into a plurality of intervals so as to obtain corresponding map grids, and setting initial values of state variables of the map grids to be 0; s3, calculating state variables in all map grids; s4, calculating a preliminary boundary of the hazard area; and S5, calculating the final boundary of the hazard area. The invention provides an asteroid impact earth hazard area boundary extraction method, which comprises the following steps of: taking a huge amount of air blast impact points and hazard radiuses obtained by an asteroid impact earth hazard evaluation engineering model as input, firstly extracting an initial boundary of a hazard area by utilizing a map grid state variable, and then obtaining a final boundary by utilizing a concave envelope algorithm; in the process, parallel computing is used for processing huge amount calculation brought by huge amount input.

Description

Technical Field

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

[0002] An asteroid impact on the Earth is a real threat to the survival and development of humanity. Planetary defense is a hot research field that has emerged internationally in the past 20 years and is also a major security requirement in our country. The assessment of the hazard of an asteroid impact on the Earth is an important part and research direction of planetary defense and is the basis for making defense decisions and civil defense on the ground. Extracting the boundary of the hazard area of an asteroid impact on the Earth helps to visually display the hazard range in a geographic information system, thus better providing support for planetary defense.

[0003] The hazard of an asteroid impact on the Earth is characterized by low probability, great harm, and randomness. Therefore, the hazard area needs to include any possible hazard points. In hazard assessment, engineering models are usually used to obtain any airburst / impact points and hazard radii within the range of input parameters such as the impact location, kinematic parameters, asteroid characteristics, entry into the atmosphere, and impact on the surface model. Each airburst / impact point and hazard radius determines the local hazard area, and the hazard assessment will give a huge number (tens of millions) of combinations of airburst / impact points and hazard radii. As Figure 10 shown ( Figure 10 in which A is the boundary of the hazard area, P i is the i th airburst / impact point, R i is the hazard radius of the i th airburst / impact point R, P bj is the boundary point on the j th hazard area), the boundary of the hazard area A is the boundary of the global area determined by the combination of the above-mentioned huge number of n airburst / impact points P i ( i = 1, 2... n) and the hazard radius R i ( i = 1, 2... n), but in practical applications, there are two difficulties in extracting the boundary of the hazard area: Firstly, the number of airburst / impact points is huge, and it is difficult to shorten the calculation time; Secondly, the shape of the boundary of the hazard area is complex, and the extracted boundary needs to have a high fidelity.

[0004] Currently, there is no work on extracting the boundary of the hazard area of an asteroid impact on the Earth in China. Summary of the Invention

[0005] An object of the present invention is to solve at least the above problems and / or deficiencies and provide at least the advantages described hereinafter.

[0006] To achieve these objects and other advantages of the present invention, a method for extracting the boundary of the hazardous area of an asteroid hitting the earth is provided, including: S1. Obtain the longitude, latitude and hazard radius data of all air bursts and / or impact points; S2. According to the preset resolution, evenly divide the longitude and latitude into multiple intervals to obtain the 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 the latitude interval in the map grid to obtain the corresponding local grid interval; In the local grid interval, determine whether the center of each grid falls within the hazard radius of the input data. If it falls within, set the corresponding grid state variable to 1; S4. Based on the judgment of whether the adjacent grid state variables are equal, determine the preliminary boundary value of whether the corresponding grid is a hazardous area, and traverse the longitude grid line to obtain the corresponding preliminary boundary data set; S5. Take the preliminary boundary data set as the input, and use the concave hull algorithm to obtain the final boundary data set of the hazardous area, and complete the boundary extraction of the hazardous area.

[0007] Preferably, in S2, the multiple intervals refer to evenly dividing the longitude in the full range of [-180°, 180°] into ni - 1 intervals, and evenly dividing the latitude in 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 represented by mgv i,j for characterization.

[0008] Preferably, in S3, the way to obtain the index value is as follows: S30. For the i th input data, based on the corresponding air burst 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 intervalB i - dB , B i + dB ; Among them, dL, dB are respectively represented by the following formulas: , In the above formulas, dL is half of 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; S31. Based on the longitude interval and latitude interval calculated from S30, calculate the index values corresponding to the longitude interval and latitude interval in the map grid.

[0009] Preferably, in S5, after taking the preliminary boundary data set as input, it further includes: S50. Divide these input data into m parts, perform the convex hull algorithm on each part of the data, and obtain the sub-hazard boundaries; S51. Combine all the sub-hazard boundaries and perform the convex hull algorithm again to obtain the final boundary data set; Among them, 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 counterclockwise or clockwise order on the map so as to be directly input into the geographic information system for boundary drawing in the later stage.

[0010] The present invention has at least the following beneficial effects: The present invention can be regarded as the post-processing process of the engineering calculation of the asteroid impact on the earth's hazard assessment, providing support for visualizing the hazard range in the geographic information system. In actual implementation, taking the huge number of air burst / impact points and hazard radii obtained from the engineering model of the asteroid impact on the earth as input, first use the map grid state variables to extract the preliminary boundary of the hazard area, and then use the convex hull algorithm to obtain the final boundary. In this process, parallel computing is used to handle the huge amount of calculation brought by the huge amount of input.

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

[0012] Figure 1 is the flow chart of the method for extracting the boundary of the asteroid impact on the earth's hazard area; Figure 2 is the flow chart of the method for calculating the map grid state variables based on openMP parallelism; Figure 3 It is a flow chart of a method for calculating the preliminary boundary of the hazard area; Figure 4 It is a flow chart of a method for calculating the final boundary of the hazard area based on region decomposition and convex hull algorithm; Figure 5 It is an example diagram of a map grid variable cloud map; Figure 6 It is an example diagram of the preliminary boundary of the hazard area; Figure 7 It is Figure 6 a partially enlarged view of; Figure 8 It is an example diagram of the final boundary of the hazard area; Figure 9 It is Figure 8 a partially enlarged view of; Figure 10 It is a schematic diagram of the boundary of the hazard area caused by an asteroid hitting the earth. Specific implementation manners

[0013] The present invention will be further described in detail below with reference to the accompanying drawings, so that those skilled in the art can implement it according to the description in the specification.

[0014] As Figure 1 shown, the present invention provides a method for extracting the boundary of the hazard area caused by an asteroid hitting the earth, including: S1. Read in the longitude, latitude and hazard radius data of all air bursts / impact points (a total of n pieces of data). The longitude, latitude and hazard radius of the air burst / impact point i in the P i th piece of data are respectively denoted as L i , B i , R i .

[0015] S2. According to the preset resolution, evenly divide the full longitude range [-180°, 180°] into ni - 1 intervals, and evenly divide the full latitude range [-90°, 90°] into nj - 1 intervals, and set the map grid status variable mgv i,j = 0 ( i = 1, 2... ni, j = 1, 2... nj).

[0016] In S2, the method of dividing the map grid is as follows: S20. Let the map grid length resolution be d , and according to the average radius of the earth r e obtain the angular resolution of the gridδ , and δ is characterized by the following formula: S21. According to the 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: In the above formula, int() is the rounding function; S22. Calculate the grid point data respectively, that is, calculate the longitude discrete value through the following formula Lc i and latitude Bc j : , , S23. Set the initial value of the state variable mgv i,j = 0 ( i = 1, 2... ni, j = 1, 2... nj).

[0017] S3. Calculate the state variables in all map grids mgv i,j ( i = 1, 2... ni, j = 1, 2... nj), if the center of the grid ( i, j ) is within the local hazard radius, then set mgv i,j to 1. On the one hand, this step only conducts calculations within the local grid domain of each input airburst / impact point, avoiding the ineffective calculation of searching the global grid; on the other hand, multi-core parallel computing is applied to shorten the calculation time; As shown in Figure 2 , the steps to obtain the map grid state variable mgv i, are as follows: S30. For the i th input data, that is, the longitude P i , latitude L i and hazard radius B i of the airburst and / or impact point R i, calculate the local longitude range L i - dL , L i + dL , and the local latitude range B i - dB , B i + dB , where dL 、 dB are represented by the following formula: , Then use the following formula to calculate the index of the longitude range in the map grid i1 , i2 , and the index of the latitude range in the map grid j1 , j2 .

[0018] , , S31. In the local map grid area formed by the longitude range and the latitude range, that is, for each grid ( iL , jL )( i ∈ i1 , i2 , j ∈ j1 , j2 ), calculate the distance iL , jL ) between the center of the grid and the air burst and / or impact point P i . If D ≤ D i , then the grid center is within the hazard range, and let R iL,jL mgv = 1. The formula for the distance D is: ​ Among them, ( x , y , z ) is the coordinate of the center of the grid ( iL , jL ) in the geocentric coordinate system. ( x i , y i , z i ) is the coordinate of the airburst and / or impact point P i in the geocentric coordinate system, Lc iL and Bc jL are the coordinates of the grid ( iL , jL ) calculated in step S22. r e is the average radius of the earth, δ is the grid angular resolution set in step S20.

[0019] S32. Traverse all input data and execute S30 - S31. During this process, use the shared memory parallel programming (openMP) technology to accelerate the loop. Specifically, add some parallel directive statements in the algorithm, such as $ompparallel, $omp do, etc.

[0020] S4. Calculate the preliminary boundary of the hazard area. Based on the adjacent map grid state variable mgv i,j value, obtain the preliminary boundary of the hazard area.

[0021] As Figure 3 shown, the steps to calculate the preliminary boundary of the hazard area are as follows: S40. For the iL th longitude grid line, find the boundary points. The characteristic of the boundary points is that the adjacent grid state variables are not equal, then consider the center of the map grid where ( iL , jL) ) is located as the preliminary boundary of the hazard, and record it.

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

[0023] S5. Calculate the final boundary of the hazard area. Use the preliminary boundary calculated as the input and obtain the final boundary using the concave hull algorithm. This step applies multi-core parallel computing to shorten the calculation time.

[0024] As Figure 4 shown, the steps for calculating the final boundary of the hazard area are as follows: S50. Using the preliminary boundary of the hazard area as input, divide this input data into m parts (corresponding to m CPU cores).

[0025] S51. Perform the convex hull algorithm on each part of the data to obtain sub-hazard boundaries.

[0026] S52. Combine all the sub-hazard boundaries and perform the convex hull algorithm again to obtain the final boundary.

[0027] Example: For a method for extracting the boundary of the hazard area caused by an asteroid hitting the Earth of the present invention, in order to more clearly illustrate the technical solution of the present invention, it is described through examples. At the same time, this example is based on the engineering calculation results of the hazard assessment of the 6th International Academy of Astronautics (IAA) Planetary Defense Conference desktop exercise, and the method of the present invention is used to extract the hazard boundary.

[0028] The specific implementation process is as Figure 1 shown and includes: S1. Read in the hazard assessment result data These data are given in the form of text or binary data, and the corresponding hazard assessment result data table (as shown in Table 1) is obtained: Table 1 Table 1 gives 9 rows of data of the hazard results, and each row represents a possible hazard result. Each column of data in Table 1 represents different variables, including the longitude L i , latitude B i and the hazard radius R i data. There are 10 million rows of data in this example.

[0029] S2. Divide the longitude and latitude grids Set the grid length resolution d to 2.5 km, then the grid angular resolution is 0.0225°, and ni and nj points are evenly distributed within the Earth's longitude and latitude ranges at this interval, where ni = 16013 and nj = 8007.

[0030] S3. As Figure 2 shown, calculate the state variables in all map grids mgv i,j ( i = 1, 2... ni, j = 1, 2... nj), and obtain as Figure 5The results shown, where Figure 5 in mgv i,j the magnitude values are represented by colors (see the legend on the right), and the range of the hazard area (green) can be seen from the figure. It can be seen from the calculation process that under a 10-core CPU, the parallel time using openMP is 5 times shorter than the serial time.

[0031] S4. Use the process as Figure 3 shown to calculate the preliminary boundary of the hazard area, and obtain the results as Figure 6 shown. From Figure 6 it can be clearly seen the boundary of the hazard area, but from Figure 7 the enlarged schematic diagram, it can be seen that there are redundant boundary points, and these data are not all arranged counterclockwise or clockwise, and cannot be used for GIS illustration.

[0032] S5. Use the process as Figure 4 shown to calculate the final boundary of the hazard area, and obtain the results as Figure 8 shown. All the calculation processes in this embodiment use 10 CPU cores of intel i9-9900K. For input data in the tens of millions, the high-fidelity hazard boundary can be obtained within 5 minutes. Comparing Figure 8 with Figure 6 , it can be seen 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, viewing Figure 9 , there is no redundancy in the boundary data. Figure 8 The data shown is ordered data, that is, all data points are arranged in a completely counterclockwise or completely clockwise order, while Figure 10 in P bj gives an example of an ordered arrangement, so that it can be directly input into a geographic information system (that is, various navigation maps on the market) for illustration.

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

[0034] Although the embodiments of the present invention have been disclosed as above, it is not limited to the applications listed in the specification and embodiments. It can be fully applied to various fields suitable for the present invention. For those familiar with the field, additional modifications can be easily made. Therefore, without departing from the general concept defined by the claims and the equivalent scope, the present invention is not limited to the specific details and the illustrated examples here.

Claims

1. A method for extracting the boundary of the hazardous area caused by an asteroid hitting the Earth, characterized in that, Including: S1. Obtain the longitude, latitude, and hazard radius data of all air bursts and / or impact points; S2. According to the preset resolution, evenly divide the longitude and latitude into multiple intervals to obtain corresponding map grids, and set the initial values of the map grid status 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, judge whether the center of each grid falls within the hazard radius of the input data. If it falls within, set the corresponding grid status variable to 1; S4. Based on the judgment of whether the adjacent grid status variables are equal, determine the preliminary boundary values of whether the corresponding grids are hazard areas, and traverse the longitude grid lines to obtain the corresponding preliminary boundary data set; S5. Use the preliminary boundary data set as input and use the concave hull algorithm to obtain the final boundary data set of the hazard area, completing the boundary extraction of the hazard area.

2. The method for extracting the boundary of the hazardous area of an asteroid impact on the Earth according to claim 1, wherein In S2, the multiple intervals refer to evenly dividing the longitude in the full range of [-180°, 180°] into ni - 1 intervals, and evenly dividing the latitude in 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 characterized by mgv i,j .

3. The method for extracting the boundary of the hazardous area of an asteroid hitting the earth according to claim 1, wherein In S3, the method for obtaining the index value is: S30. For the i th input data, based on the corresponding air burst 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 interval B i - dB , B i +dB ; Among them, dL, dB are respectively characterized by the following formulas: ， In the above formula, dL is half of 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; 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 hazardous area caused by an asteroid hitting the Earth according to claim 1, wherein In S5, after using the preliminary boundary data set as input, it further includes: S50. Divide these input data into m parts, perform the concave hull algorithm on each part of the data, and obtain the sub-hazard boundaries; S51. Combine all the sub-hazard boundaries and perform the concave hull algorithm again to obtain the final boundary data set; Among them, 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 counterclockwise or clockwise order on the map so that they can be directly input into the geographic information system for boundary drawing in the later stage.

Images