Astronomical cross identification method based on spherical trapezoidal division

Abstract

The invention discloses an astronomical cross identification method based on spherical trapezoidal division, and the method comprises the steps: S1, taking a unit spherical surface as a celestial sphere, dividing the celestial sphere into a plurality of strips parallel to an equatorial line, and determining the right ascension span of each strip for the transverse segmentation of the strip; s2, dividing each strip into a same spherical isosceles trapezoid with two waists perpendicular to the equator, and taking the same spherical isosceles trapezoid as a sky area; s3, determining an index of each sky area; s4, dividing a star catalogue A and a star catalogue B, and creating a zipper type hash table based on the row number of the star catalogue B; s5, mapping each celestial body in the star catalogue B into a linked list in a hash table in a multi-thread parallel mode, and inserting celestial body data into the linked list by means of a two-stage insertion strategy; s6, mapping each celestial body in the star catalogue A into a linked list in a group of hash tables in a multi-thread parallel mode, and traversing to obtain celestial body IDs of which the spherical distance is smaller than the error radius; the method not only can completely avoid the problem of result omission, but also can greatly improve the astronomical cross identification calculation performance.

Description

Technical Field

[0001] The present invention relates to the field of astronomical information technology in astronomy and computer science, and particularly relates to an astronomical cross-identification method based on spherical trapezoid division. Background Art

[0002] In the research of astronomical data fusion, cross-identification calculation refers to an important operation of matching the corresponding data of the same celestial body or region in multiple astronomical observation catalog datasets. When observing with different types of observation devices, due to different visible bands or imaging methods, the same celestial body may be observed by one or several of these methods. In addition, under different observation methods, observation projects or observation times, due to factors such as errors generated during observation and atmospheric interference, parameters such as the coordinates of the same celestial body in different astronomical observation datasets may have certain differences. Through the cross-identification process, different observation catalogs can be matched and fused to obtain multiple observation data of the same celestial body, and the observation data from different sources can be integrated to form a more comprehensive and accurate astronomical catalog.

[0003] With the development of astronomical observation means and the implementation of large-scale sky survey projects, astronomical catalog data has been rapidly accumulating and growing. The current scale of large catalogs has reached the billion level. To accelerate cross-identification calculation, existing work often uses sky region division and indexing methods to perform local matching on each divided region. However, existing sky region division and indexing methods all have at least one or more of the following problems: First, there will be a boundary leakage problem, resulting in missing identification results; Second, static or hierarchical division methods cannot flexibly adapt to unevenly distributed or different-scale catalog data, affecting the efficiency of cross-identification calculation; Third, after the sky region index is mapped to a hash table container, physically adjacent regions are not adjacent on the hash table, and the unfriendly memory access operations will reduce performance. Therefore, it is necessary to design an efficient astronomical cross-identification method based on a new sky region division method. Summary of the Invention

[0004] The purpose of the present invention is to provide an astronomical cross-identification method based on spherical trapezoid division to solve the above technical problems.

[0005] To this end, the technical solution of the present invention is as follows:

[0006] An astronomical cross-identification method based on spherical trapezoid division, the steps are as follows:

[0007] S1. Taking the unit sphere as the celestial sphere, dividing the unit sphere into a plurality of strips parallel to the equator line, and for each strip, determining the right ascension span for transverse cutting of the strip;

[0008] S2. For each strip, divide it into identical spherical isosceles trapezoids with their two waists perpendicular to the equator. Each spherical isosceles trapezoid is a celestial region. Among them, the base length of the spherical isosceles trapezoid is the maximum right ascension span calculated for this strip in step S1.

[0009] S3. For each celestial region divided in step S2, use the strip number where it is located and the celestial region number within the strip together as the index of this celestial region.

[0010] S4. For the two star catalogs to be cross-identified, take one star catalog as star catalog A and the other as star catalog B, and create a zippered hash table of appropriate size based on the number of rows in star catalog B.

[0011] S5. In a multi-threaded parallel manner, for each celestial body in star catalog B, find the index value of the celestial region where it is located, and use the rectangular aggregation hashing method to map the index value to a linked list in the hash table. Then, with the help of a two-stage insertion strategy, insert the celestial body data into the linked list.

[0012] S6. In a multi-threaded parallel manner, for each celestial body in star catalog A, according to the celestial region where it is located, determine a set of index values of the celestial regions to be queried by the four-point method, and map this set of index values to a set of linked lists in the hash table using the rectangular aggregation hashing method and traverse them. Output the pairs of celestial body IDs with the obtained spherical distance less than the error radius.

[0013] Further, the declination span of the strip is 2 times the error radius.

[0014] Further, in step S1, the calculation formula for the maximum right ascension span Δra is:

[0015]

[0016] In the formula, r is the error radius given for the cross-identification task, and dec zone is the declination of the strip boundary farthest from the equator line in the strip.

[0017] Further, in step S3, the specific setting method for the index of the region is:

[0018] S301. From the S pole to the N pole of the unit sphere, starting from 0, sequentially and continuously number the strip numbers divided in step S1.

[0019] S302. For each strip, starting from the 0-degree meridian, sequentially and continuously number the spherical trapezoid celestial region numbers divided in step S2.

[0020] S303. For each celestial region, the strip where it is located has the number id obtained from step S301 row , and the number id obtained from step S302col , then the number of this celestial region is the binary tuple (id row , id col ).

[0021] Furthermore, in step S4, the specific method for creating the hash table is described in detail as follows:

[0022] S401. Find the smallest prime number not less than the number of rows of star catalog B as the size size of the hash table table ;

[0023] S402. Create and arrange size table linked lists in sequence to form a chaining hash table.

[0024] Furthermore, the specific processing steps of step S5 are as follows:

[0025] S501. For the given number of worker threads count thread , divide star catalog B into count thread continuous data segments on average, and each thread is assigned a data segment as its insertion task;

[0026] S502. Create a container matrix C composed of linear lists, and its size is count thread ×count thread ;

[0027] S503. For each single worker thread, generate an insertion request for each celestial body in its insertion task in sequence, and insert the insertion request into the container matrix;

[0028] S504. When a certain number of insertion requests are accumulated in the container matrix, all worker threads pause generating insertion requests and instead traverse the container matrix and insert the insertion requests into the linked list at the corresponding position in the hash table;

[0029] S505. Repeat step S503 and step S504 until all insertion tasks of a single worker thread are completed.

[0030] Furthermore, in step S503, the specific operation steps for each single worker thread to insert the insertion request into the corresponding position of the container matrix are as follows:

[0031] S5031. Traverse each celestial body in the insertion task of a single worker thread in sequence;

[0032] S5032. For each celestial body, find the index binary tuple (id row , id col ) of the celestial region where it is located according to its spherical coordinates, and its specific calculation formula is:

[0033]

[0034] Among them, ra is the right ascension of the celestial body, dec is the declination of the celestial body, and r is the error radius. is the maximum right ascension span calculated in step S1 for the strip where the celestial body is located with the serial number id row ;

[0035] S5033. Adopt the rectangular aggregation hashing method to map the index binary group of the celestial sphere determined in step S5032 to an integer id table within [0, size hash , which is used as its insertion position in the hash table. The specific calculation formula of the rectangular aggregation hashing method is:

[0036]

[0037] Among them, represents converting the integer parts of and into binary and performing bit splicing;

[0038] S5034. Calculate the length L thread for evenly dividing the hash table into count chunk consecutive segments. Its specific calculation formula is:

[0039]

[0040] S5035. For the insertion request generated by thread id thread to be inserted into the hash table at position id hash , calculate the segment id chunk to which this insertion position belongs in the hash table; furthermore, insert the celestial body data of the star catalog B and id hash together into the linear list at the position (id thread , id chunk ) of the container matrix.

[0041] Furthermore, the specific implementation method of step S504 is: for thread id thread , it sequentially traverses the linear lists in the (0, id thread ) to (count thread - 1, id thread ) column of the container matrix, takes out the celestial body data and id hash , and inserts the celestial body data in the star catalog B into the linked list at the position id hash of the hash table.

[0042] Furthermore, the specific implementation steps of step S6 are:

[0043] S601. For a given number of worker threads count thread , evenly divide the star catalog A into count thread consecutive data segments, and allocate one data segment to each thread as its query task;

[0044] S602. For a single worker thread, sequentially traverse the celestial bodies in its query task. For each celestial body, use the four - point method to determine a set of sky area indices that need to be judged one by one with this celestial body;

[0045] S603. Use the rectangular aggregation hashing method to calculate the hash table positions corresponding to this set of indices, and obtain a hash table position set P after removing duplicate positions;

[0046] S604. Sequentially traverse the linked list of each position in each hash table in the set P. For each celestial body data in the star catalog B during the traversal process, calculate its spherical distance from the current celestial body in the star catalog A. If the spherical distance is less than the given error radius r, then output the ID pair of the two celestial bodies.

[0047] Furthermore, the specific implementation steps of the four - point method in step S602 are as follows:

[0048] S6021. Read a celestial body in the star catalog A in the query task, whose spherical coordinates are (α, δ), and obtain its maximum right ascension span Δra calculated according to step S1 by looking up a table;

[0049] S6022. Calculate the coordinates of the four corner points of the spherical square with (α, δ) as the center and Δra as the side length: and

[0050] S6023. According to the four spherical coordinates obtained in step S6022, calculate the indices of the sky areas where they are located respectively, thereby determining a set of sky area indices that need to be judged one by one with this celestial body.

[0051] Compared with the prior art, the astronomical cross - identification method based on spherical trapezoid division can flexibly adapt to different celestial body distribution situations and cross - identification tasks with different error radii. While ensuring calculation accuracy and integrity, it makes full use of the parallel computing power of the computer, improves the calculation performance of the cross - identification algorithm, and has good multi - thread parallel acceleration ability, realizing the improvement of the overall calculation performance of astronomical cross - identification. BRIEF DESCRIPTION OF THE DRAWINGS

[0052] Figure 1 is a flowchart of the astronomical cross - identification method based on spherical trapezoid division of the present invention;

[0053] Figure 2This is the flowchart of the specific implementation method of step S5 in the astronomical cross-identification method based on spherical trapezoid division of the present invention;

[0054] Figure 3 This is the flowchart of the specific implementation method of step S6 in the astronomical cross-identification method based on spherical trapezoid division of the present invention. Specific Embodiment

[0055] The following further describes the present invention in conjunction with the accompanying drawings and specific embodiments, but the following embodiments are by no means any limitation to the present invention.

[0056] See Figure 1 , the specific implementation steps of the astronomical cross-identification method based on spherical trapezoid division are described as follows.

[0057] S1. Taking the unit sphere as the celestial sphere, divide the unit sphere into a number of strips parallel to the equator line, and the declination span of the strip is 2 times the error radius; for each strip, calculate the right ascension span for horizontally dividing the strip according to the formula;

[0058] Among them, the calculation formula for the maximum right ascension span Δra is:

[0059]

[0060] In the formula, r is the error radius given by the cross-identification task, and dec zone is the declination of the strip boundary farthest from the equator line in the strip; r is set to This value is also equal to the radian representation of 1 arcsecond; at this time, the celestial sphere is divided into 648,000 strips, so a real number array needs to be created to store the Δra values of each strip;

[0061] The advantage of this division method is that for celestial bodies located in a certain spherical trapezoid, celestial bodies with a spherical distance less than the error radius r from it can only be distributed within the same spherical trapezoid and no more than 3 adjacent spherical trapezoids, which can relatively simply solve the problem of omission of identification results.

[0062] It should be noted that this method is based on the astronomical standard celestial coordinate system, with the geocenter as the origin, and uses right ascension and declination to locate any point on the unit sphere (the abstract model of the celestial sphere with a radius of 1). Specifically, it is defined as: right ascension is the eastward azimuth angle measured from the vernal equinox direction (the northward intersection of the celestial ecliptic and the equator) in the equatorial plane, and the value range is [0°, 360°); declination is the polar angle with the celestial equatorial plane as the reference plane and the north celestial pole direction as positive, and the value range is [-90°, 90°], where the positive and negative values respectively represent the northern and southern hemisphere celestial regions.

[0063] In addition, the key reference objects are defined as follows: North Celestial Pole (N Pole): The only position point on the unit sphere with a declination value of +90°. Its spatial position is defined by the intersection of the northward extension of the Earth's axis of rotation and the celestial sphere. Right ascension has no azimuth significance at this point; South Celestial Pole (S Pole): The only position point on the unit sphere with a declination value of -90°, corresponding to the intersection of the southward extension of the Earth's axis of rotation and the celestial sphere. Right ascension has no azimuth significance at this point; Equator line: The set of all position points on the unit sphere with a declination value of 0°, forming the intersection line of the equatorial plane and the celestial sphere. Right ascension is continuously distributed in the range of [0°, 360°) along this line, and its geometric shape is a perfect great circle passing through the center of the sphere.

[0064] S2. For each strip, divide it into identical spherical isosceles trapezoids with two waists perpendicular to the equator. Each spherical isosceles trapezoid is a celestial region; where the base length of the spherical isosceles trapezoid is the maximum right ascension span calculated for the strip in step S1.

[0065] S3. For each celestial region divided in step S2, use the strip number where it is located and the celestial region number within the strip together as the index of this celestial region.

[0066] In step S3, the specific setting method of the index of the celestial region is as follows:

[0067] S301. From the S Pole to the N Pole of the unit sphere, starting from 0, continuously number the strips divided in step S1 in sequence; where if there is an incomplete strip near the N Pole, it is also regarded as an independent strip and numbered;

[0068] S302. For each strip, starting from the 0-degree meridian, continuously number the spherical trapezoid celestial regions divided in step S2 in sequence; where if there is an incomplete spherical trapezoid at the 360-degree meridian, it is also regarded as an independent celestial region and numbered;

[0069] S303. For each celestial region, the strip where it is located has the number id obtained from step S301 row , and the number id obtained from step S302 col , then the number of this celestial region is the binary tuple (id row , id col ).

[0070] S4. For two star catalogs to be cross-identified, take one of the star catalogs as star catalog A and the other as star catalog B, and create a zippered hash table of appropriate size based on the number of rows in star catalog B.

[0071] In step S4, the specific creation method of the hash table is described as follows:

[0072] S401. Find the smallest prime number not less than the number of rows of star catalog B as the size of the hash table, denoted as size. table ;

[0073] S402. Create and arrange size table linked lists in sequence to form a chained hash table.

[0074] In specific applications, to reduce the size of the hash table and improve efficiency, it is usually preferred to select the star catalog with more rows as star catalog A and the star catalog with fewer rows as star catalog B according to experience; each position of the hash table uses a linked list without a head node, and each node in the linked list stores a pointer to the next node and the information necessary for performing cross-identification calculations to form a tuple (ra, dec, id), which are the right ascension, declination, and integer ID of the celestial body respectively.

[0075] S5. In a multi-threaded parallel manner, for each celestial body in star catalog B, find the index value of its celestial region, and use the rectangular aggregation hash method to map the index value to a linked list in the hash table, and then insert the celestial body data into the linked list by means of a two-stage insertion strategy.

[0076] See Figure 2 , the specific implementation steps of this step S5 are as follows:

[0077] S501. For a given number of worker threads count thread , divide star catalog B into count thread continuous data segments on average, and allocate one data segment to each thread as its insertion task;

[0078] S502. Create a container matrix C composed of linear lists, with a size of count thread ×count thread to be used for temporarily storing insertion requests;

[0079] S503. For each single worker thread, generate an insertion request for each celestial body in its insertion task in sequence, and insert the insertion request into the corresponding position of the container matrix.

[0080] Among them, in this step S503, the specific method for each single worker thread to insert the insertion request into the corresponding position of the container matrix is as follows:

[0081] S5031. Traverse each celestial body in the insertion task of a single worker thread in sequence;

[0082] S5032. For each celestial body, find the index tuple (id row , id col ) of its celestial region according to its spherical coordinates, and its specific calculation formula is:

[0083]

[0084] Wherein, ra is the right ascension of the celestial body, dec is the declination of the celestial body, and r is the error radius, which is in this embodiment, for the strip numbered id where the celestial body is located row the maximum right ascension span that has been calculated and saved in step S1;

[0085] S5033. Adopt the rectangular aggregation hashing method to map the index binary tuple of the sky area determined in step S5032 to an integer id within [0, size table -1] as its insertion position in the hash table. The specific calculation formula of the rectangular aggregation hashing method is: hash

[0086]

[0087] Wherein, means converting the integer parts of and into binary and performing bit concatenation;

[0088] The advantage of this rectangular aggregation hashing method is that the spherical trapezoidal regions adjacent in position on the sphere are combined into larger approximately rectangular-shaped block areas. The spherical trapezoidal regions in the same block area will be mapped to adjacent positions on the hash table after rectangular aggregation hashing, which can improve the locality of memory access and enhance the computing performance.

[0089] S5034. Calculate the length L of dividing the hash table into count thread consecutive segments. Its specific calculation formula is: chunk

[0090]

[0091] S5035. For the insertion request generated by thread id thread to be inserted into the hash table at position id hash , calculate the segment to which this insertion position belongs in the hash table Furthermore, insert the celestial body data of the star catalog B and id hash together into the linear list at the position (id thread , id chunk ) of the container matrix.

[0092] S504. After a certain number of insertion requests have accumulated in the container matrix, all worker threads pause generating insertion requests. Then all threads start traversing the container matrix and insert the insertion requests into the linked list at the corresponding position in the hash table. Among them, according to experience, the certain number is generally set as follows: when each single thread has generated 65,536 insertion requests, the generation of insertion requests is paused at this time.

[0093] In this step S504, the specific operation steps for inserting the insertion requests into the hash table are as follows:

[0094] For thread id thread , it sequentially traverses the linear list in the (0, id thread ) to (count thread - 1, id thread ) column of the container matrix, extracts the celestial body data and id hash , and inserts the celestial body data in star catalog B into the linked list at the id hash position in the hash table. In this embodiment, head insertion is adopted to ensure that each insertion operation is a constant-time operation.

[0095] S505. Repeat the request generation process in step S503 and the insertion operation in step S504 until all insertion tasks of a single worker thread are completed.

[0096] In this step S5, a two-stage insertion strategy of first inserting into the container matrix and then inserting into the linked list is adopted. The insertion process is decomposed into two stages and is repeatedly executed alternately until all celestial bodies in star catalog B are inserted into the hash table. The advantage of this strategy is that when the scale of the star catalog is large, the size table of the hash table is also relatively large. At this time, for a single insertion operation, the probability of collision is very low. Therefore, it is inefficient to use atomic operations for each insertion operation. However, according to the law of large numbers, for a large number of insertion operations, collisions are almost certain to occur. Therefore, measures must be taken to avoid parallel collisions. The above two-stage insertion strategy can avoid write conflicts with relatively low additional performance overhead by introducing a container matrix.

[0097] S6. In a multi-threaded parallel manner, for each celestial body in star catalog A, according to its sky area, a set of index values of the sky area to be queried is determined by the four-point method, and this set of index values is mapped into a set of linked lists in the hash table using the rectangular aggregation hashing method and traversed, and the pairs of celestial body IDs with the obtained spherical distance less than the error radius are output.

[0098] The specific implementation steps of step S6 are as follows:

[0099] S601. For the given number of worker threads count thread , divide star catalog A into count equal partsthread a continuous data segment, and each thread is assigned a data segment as its query task;

[0100] S602. For a single worker thread, traverse the celestial bodies in its query task in sequence. For each celestial body, use the four-point method to determine a set of sky area indexes that need to be judged one by one with this celestial body;

[0101] The specific implementation steps of the four-point method described in step S602 are as follows:

[0102] S6021. Read a celestial body in a star catalog A in the query task, whose spherical coordinates are (α, δ), and obtain its maximum right ascension span Δra calculated according to step S1 by looking up the table;

[0103] S6022. Calculate the coordinates of the 4 corner points of the spherical square with (α, δ) as the center and Δra as the side length. Specifically, the coordinates of these 4 points are respectively and

[0104] S6023. Calculate the indexes of the sky areas where the 4 spherical coordinates obtained in step S6022 are located respectively according to step S5032, thereby determining a set of sky area indexes that need to be judged one by one with this celestial body.

[0105] S603. Use the rectangular aggregation hashing method described in step S5033 to calculate the hash table positions corresponding to this set of indexes, and obtain a hash table position set P after removing duplicate positions;

[0106] S604. Traverse the linked list of each position in each hash table in set P in sequence. For each celestial body data in star catalog B during the traversal process, calculate its spherical distance from the current celestial body in star catalog A. If the spherical distance is less than the given error radius r, then output the ID pair of the two celestial bodies.

[0107] In the actual identification task, after this embodiment is executed, the ID pairs of celestial bodies with a spherical distance less than the given error radius will be continuously output.

[0108] In practical applications, the identification method of this application has the following advantages compared with existing identification methods: 1) By dividing the spherical surface into trapezoidal regions that are easy to process and obtaining the correct sky area to be queried with the four-point method, the problem of missing sources in other existing methods is effectively solved; 2) The division process is non-hierarchical, and the sky area index is easier to calculate, which can improve the sky area indexing process and flexibly adapt to different celestial body distributions and error radii; 3) A rectangular aggregation hashing strategy is proposed, which can further optimize the access performance of the hash table by utilizing the local aggregation characteristics of astronomical observation catalogs, reduce the number of random memory accesses, and improve the overall performance of the algorithm; 4) In terms of multi-threaded parallel acceleration, this chapter proposes a two-stage insertion method to solve the small-probability hash collision problem in the insertion stage in a multi-threaded environment with relatively low additional performance overhead.

[0109] In terms of the shape of the sky area division, the sky area division described in steps S1 and S2 of this application can divide the spherical surface into spherical trapezoids that are easy to judge the positional relationship with the spherical small circle, and the divided trapezoids are all of equal height. When the latitude changes, except for the change in width, the trapezoids will not produce shape distortions that affect the cross-identification accuracy. In existing division methods, HEALPix divides the spherical surface into rhombuses, and HTM divides the spherical surface into triangles. On the one hand, it is difficult to judge the positional relationship between the small pieces and the spherical small circle due to the irregularity of the small piece shapes themselves. On the other hand, the shapes of the small pieces at different latitudes will be different, and this shape distortion will be more serious near the poles. These two aspects together lead to the problem of difficult-to-handle missing sources, so it is difficult to be applied to the precise matching process based on spherical positions.

[0110] In terms of the division process and sky area index, step S1 of this application pre-stores the processes that require complex trigonometric function and square root function calculations in a table, and only simple operations such as floating-point division are required in the subsequent process of positioning the sky area according to spherical coordinates; the sky area division and numbering processes described in steps S1, S2, and S3 are only related to the error radius, and can be flexibly divided into different numbers and sizes of small pieces to adapt to different data scales or error requirements. For sky area division and indexing methods such as HEALPix and HTM, relatively complex trigonometric functions need to be calculated during operation, and the spherical surface can only be divided according to a certain hierarchy, and the number of sky areas can only be adjusted in integer powers of 4, making it difficult to flexibly adapt to different celestial body distribution densities and different identification error radii.

[0111] In terms of the hashing process, since using the four-point method for each celestial body in catalog A will generate 4 sky regions to be accessed, when using the traditional locality-sensitive hashing method, this process will always result in 4 memory accesses at random positions. At this time, the random access of memory will become a performance bottleneck, restricting the overall performance of cross-identification. This application proposes a rectangular aggregation hashing strategy, which can better utilize the local aggregation property of astronomical observation catalog data. After actual test verification, compared with locality-sensitive hashing, the performance of the query process of cross-identification using the rectangular aggregation hashing strategy has been significantly improved. The main idea of the rectangular aggregation hashing method is to use the characteristics of the sky region division method proposed in this application, take a larger approximate rectangle composed of multiple adjacent spherical trapezoids as a larger sky region, and map them to adjacent positions on the hash table when performing hash mapping.

[0112] In terms of multi-threaded parallel acceleration, step S5 involves the insertion operation of the hash table. Therefore, it is necessary to avoid multiple threads simultaneously performing insertion operations on the linked list at the same position of the hash table. When multiple threads simultaneously perform linked list insertion operations, the most common method to resolve conflicts is to use mutex locks to protect the critical section. However, in large-scale astronomical cross-identification calculations, due to the excessively large amount of data to be inserted, if a mutex lock is used for each catalog record and a lock-unlock operation is performed once, the algorithm performance will be extremely poor. An optimization for mutex locks is to utilize the atomic instructions provided in modern computer architectures to implement linked list insertion in a lock-free manner through atomic operations. The underlying principle of atomic operations is usually to lock the bus at the hardware level. Therefore, although atomic operations are lighter than mutex locks, they are still much slower than regular instructions. In addition, when the catalog scale is large, the size of the hash table is also relatively large. At this time, for a single insertion operation, the probability of conflict is very low. Therefore, using atomic operations for each insertion operation is inefficient. However, for a large number of insertion operations, conflicts are almost certain to occur, and measures must be taken to avoid parallel conflicts. Therefore, this application proposes a two-stage insertion method, as Figure 1 shown, that is, before actually inserting into the hash table, first summarize the elements to be inserted, temporarily store the insertion requests in the container matrix, and then insert them into the hash table in parallel without interference. This method avoids the performance overhead of being forced to use global atomic operations due to low-probability conflict events and achieves higher computational efficiency.

Claims

1. An astronomical cross-verification method based on spherical trapezoidal division, characterized in that: Here are the steps: S1. Taking the unit sphere as the celestial sphere, the unit sphere is divided into several strips parallel to the equator, and for each strip, the right ascension span used for horizontally dividing the strip is determined; S2. For each strip, divide it into two identical spherical isosceles trapezoids whose sides are perpendicular to the equator, each spherical isosceles trapezoid being a sky region; wherein the base length of the spherical isosceles trapezoid is the maximum right ascension span of the strip calculated in step S1; S3, for each sky area divided in step S2, the strip number where the sky area is located and the sky area number in the strip are used as the index of the sky area; S4. For the two star tables to be cross-certified, one of the star tables is taken as star table A, and the other star table is taken as star table B, and a zipper hash table of equivalent size is created based on the number of rows in star table B; S5. For each celestial body in the star table B, the index value of the celestial region where the celestial body is located is calculated in a multi-threaded parallel manner, and the index value is mapped into a linked list in the hash table using the rectangular cluster hash method, and then the celestial body data is inserted into the linked list using the two-stage insertion strategy; S6. For each celestial body in the star catalog A, a set of index values ​​of the celestial region to be queried is determined by the four-point method according to the celestial region where the celestial body is located in a multi-threaded parallel manner, and the set of index values ​​is mapped into a linked list in a hash table using the rectangular cluster hash method and traversed, and the celestial body ID pairs whose spherical distance is less than the error radius are output.

2. The astronomical cross-verification method based on spherical trapezoidal division according to claim 1 is characterized in that: In step S1, the declination span of the strip is twice the error radius.

3. The astronomical cross-verification method based on spherical trapezoidal division according to claim 1, characterized in that: In step S1, the calculation formula of the maximum right ascension span Δra is: Where r is the error radius given by the cross-verification task, dec zone is the declination of the strip boundary farthest from the equator.

4. The astronomical cross-verification method based on spherical trapezoidal division according to claim 1, characterized in that: In step S3, the specific method for setting the zone index is: S301, numbering the strips divided in step S1 in sequence from the S pole to the N pole of the unit sphere starting from 0; S302, for each strip, numbering the spherical trapezoidal sky areas divided in step S2 in sequence starting from the 0 degree longitude; S303: for each sky region, the strip in which it is located is numbered as id according to the number obtained in step S301. row , the number obtained in step S302 is id col , then the number of the area is a two-tuple (id row ,id col ).

5. The astronomical cross-verification method based on spherical trapezoidal division according to claim 1, characterized in that: In step S4, the specific method of creating the hash table is described as follows: S401. Find the smallest prime number that is not less than the number of rows in the star table B, and use it as the size of the hash table. table ; S402, create and arrange sizes in sequence table A linked list forms a zipper hash table.

6. The astronomical cross-verification method based on spherical trapezoidal division according to claim 1, characterized in that: The specific processing steps of step S5 are: S501. For a given number of worker threads count thread , divide the star list B into count thread consecutive data segments, each thread is assigned a data segment as its insertion task; S502, create a container matrix C consisting of linear lists, whose size is count thread ×count thread ; S503, for each single working thread, generating an insertion request for each celestial body in its insertion task in turn, and inserting the insertion request into the container matrix; S504, when a certain number of insert requests are accumulated in the container matrix, all worker threads stop generating insert requests, switch to traversing the container matrix and actually inserting the insert requests into the linked list at the corresponding position of the hash table; S505. Repeat step S503 and step S504 until all insertion tasks of a single working thread are completed.

7. The astronomical cross-verification method based on spherical trapezoidal division according to claim 6, characterized in that: In step S503, the specific operation steps for each single working thread to insert the insertion request into the corresponding position of the container matrix are: S5031, traverse each celestial body inserted into the task by a single working thread in turn; S5032. For each celestial body, find the index tuple (id row ,id col ), the specific calculation formula is: Among them, ra is the right ascension of the celestial body, dec is the declination of the celestial body, r is the error radius, Δra idrow The celestial body is numbered as id row The maximum right ascension span of the strip calculated in step S1; S5033: Use the rectangular cluster hashing method to map the index tuple of the sky area determined in step S5032 to [0, size table -1] in integer id hash , as its insertion position in the hash table, the specific calculation formula of the rectangular aggregation hash method is: in, Indicates that and The integer part of is converted to binary and bit-joined; S5034, calculate and divide the hash table into count thread The length of the consecutive segments is L chunk , the specific calculation formula is: S5035. For thread id thread The generated id should be inserted into the hash table hash The insertion request for the position is calculated to calculate the segment ID of the insertion position in the hash table. chunk ; Then, the celestial body data and id of the star catalog B hash Inserted into the container matrix together (id thread ,id chunk ) positions in a linear table.

8. The astronomical cross-verification method based on spherical trapezoidal division according to claim 6, characterized in that: The specific implementation method of step S504 is: for thread id thread , which traverses the container matrix (0, id thread )~(count thread -1,id thread ) Linear table of this column, extract celestial body data and ID hash , and insert the celestial body data in star table B into the hash table id hash Position in the linked list.

9. The astronomical cross-verification method based on spherical trapezoidal division according to claim 1, characterized in that: The specific implementation steps of step S6 are: S601. For a given number of worker threads, count thread , divide the star list A into count thread continuous data segments, and each thread is assigned a data segment as its query task; S602: For a single working thread, traverse the celestial bodies in its query task in turn, and for each celestial body, use the four-point method to determine a group of sky region indexes that need to be judged one by one with the celestial body; S603, using a rectangular cluster hashing method, calculating the hash table position corresponding to the group index, and obtaining a hash table position set P after removing duplicate positions; S604, traverse the linked list of each position of each hash table in the set P in turn, and for each celestial body data in the star table B during the traversal process, calculate its spherical distance with the celestial body in the current star table A, and if the spherical distance is less than the given error radius r, output the ID pair of the two celestial bodies.

10. The astronomical cross-verification method based on spherical trapezoidal division according to claim 9, characterized in that: The specific implementation steps of the four-point method in step S602 are: S6021, read a celestial body in the star table A in the query task, whose spherical coordinates are (α, δ), and obtain the maximum right ascension span calculated according to step S1 as Δra by looking up the table; S6022. Calculate the coordinates of the four corner points of a spherical square with (α, δ) as the center and Δra as the side length: and S6023. According to the four spherical coordinates obtained in step S6022, the indexes of the celestial regions are calculated respectively, thereby determining a group of celestial region indexes that need to be judged one by one with the celestial body.

Images