CN109558568A - A kind of target RCS calculation method based on CUDA - Google Patents
A kind of target RCS calculation method based on CUDA Download PDFInfo
- Publication number
- CN109558568A CN109558568A CN201811389168.0A CN201811389168A CN109558568A CN 109558568 A CN109558568 A CN 109558568A CN 201811389168 A CN201811389168 A CN 201811389168A CN 109558568 A CN109558568 A CN 109558568A
- Authority
- CN
- China
- Prior art keywords
- grid
- pointer
- triangular surface
- max
- radar
- 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.)
- Granted
Links
- 238000004364 calculation method Methods 0.000 title claims abstract description 31
- HPTJABJPZMULFH-UHFFFAOYSA-N 12-[(Cyclohexylcarbamoyl)amino]dodecanoic acid Chemical compound OC(=O)CCCCCCCCCCCNC(=O)NC1CCCCC1 HPTJABJPZMULFH-UHFFFAOYSA-N 0.000 title abstract 4
- 238000000034 method Methods 0.000 claims abstract description 33
- 239000013598 vector Substances 0.000 claims description 24
- 239000011159 matrix material Substances 0.000 claims description 12
- 230000006870 function Effects 0.000 claims description 5
- 230000010287 polarization Effects 0.000 claims description 4
- 230000009466 transformation Effects 0.000 claims description 2
- 238000010276 construction Methods 0.000 abstract 1
- 230000001186 cumulative effect Effects 0.000 abstract 1
- 238000005192 partition Methods 0.000 abstract 1
- 238000004088 simulation Methods 0.000 description 5
- 230000007547 defect Effects 0.000 description 2
- 230000005672 electromagnetic field Effects 0.000 description 2
- 238000006467 substitution reaction Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 239000003086 colorant Substances 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 238000009877 rendering Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/88—Radar or analogous systems specially adapted for specific applications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T1/00—General purpose image data processing
- G06T1/20—Processor architectures; Processor configuration, e.g. pipelining
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Data Mining & Analysis (AREA)
- Computational Mathematics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Databases & Information Systems (AREA)
- Electromagnetism (AREA)
- Computing Systems (AREA)
- Algebra (AREA)
- Computer Networks & Wireless Communication (AREA)
- Computer Graphics (AREA)
- Geometry (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
The present invention relates to electromagnetism computing technique fields, disclose a kind of target RCS calculation method based on CUDA.This method comprises the following steps: first by the Surface Partition of radar mockup at multiple Triangular patch, CUDA is recycled to project each Triangular patch to grid needed for the x'Oy' plane of radar fix system and determining all face elements of covering, by calculating the corresponding depth value of each Triangular patch and bin number, each Triangular patch is put into corresponding grid to obtain clear zone face element, the RCS value of whole clear zone face elements is carried out the cumulative RCS for obtaining radar target by the RCS value for finally utilizing the corresponding clear zone face element of each grid of CUDA parallel computation.Face element can be put into corresponding grid using the raster data structure of no lock when obtaining clear zone face element by the present invention, it can be avoided when being compared each face element with remaining all face element and face element being put into grid using there is lock construction, and then can reduce the complexity of calculating, it reduces and calculates the time.
Description
Technical Field
The invention relates to the technical field of electromagnetic computing, in particular to a target Radar scattering Cross-sectional area (RCS) computing method based on a general parallel computing Architecture (CUDA).
Background
When an object in space is irradiated by electromagnetic waves, energy of the object is scattered towards all directions, a scattered field and an incident field are superposed to form a space total electromagnetic field, the actual scene is that the radar transmits the electromagnetic waves to detect a target, the space total electromagnetic field can be expressed as the power density of the radar waves multiplied by an effective area, and the effective area is the target RCS. The Physical Optics (PO) method is a high-frequency approximation of RCS of a computational object, and represents a scattering field by the integral of induced current on the surface of a scattering body.
The PO method is characterized in that the surface of the three-dimensional model corresponding to the target is dispersed into a plurality of surface elements, the RCS of each surface element is calculated and superposed, and the RCS of the whole target is obtained. The PO method assumes that the portion of the target not illuminated by the radar waves, i.e., the shaded portion, does not contribute to the RCS of the entire target. Therefore, the PO method only needs to calculate the RCS of the bright area of the three-dimensional model corresponding to the target, so that the three-dimensional model needs to be divided into a plurality of surface elements, and then the bright area surface elements are obtained through the shielding judgment among the surface elements of the three-dimensional model. When the occlusion judgment is performed, the existing method needs to compare each surface element with all other surface elements, and the calculation complexity is the square of the number of the surface elements. When the RCS of an electrically large complex target is calculated, the bins need to be divided more finely, the number of the bins is increased greatly, which results in complicated calculation and greatly increased calculation time. In order to solve the above problems, some researchers propose to assign different colors to each bin of the three-dimensional model corresponding to the target and then use an Open Graphics library (OpenGL) rendering result to acquire the information of the bin in the bright area, but in this method, when the included angle between the sight line and the bin is small, the condition that two regions of bins are missed is caused.
In summary, the inventors found that in the prior art, in the process of calculating the RCS of the three-dimensional model corresponding to the high-frequency target by using the PO method, the occlusion determination between the bins is complex to calculate, the calculation time is long, and false determination of the bright-area bin may occur.
Disclosure of Invention
In view of this, embodiments of the present invention provide a method for calculating a target RCS of a CUDA, which can place bins into corresponding grids by using a grid data structure without lock when obtaining bright-area bins, and can avoid comparing each bin with all other bins and using a lock structure when placing bins into grids, thereby reducing the complexity of calculation and reducing the calculation time.
In order to achieve the above purpose, the embodiment of the invention adopts the following technical scheme:
the method for calculating the target RCS based on the CUDA comprises the following steps:
step 1, dividing the outer surface of a radar target three-dimensional model into a plurality of triangular surface elements, and converting vertex coordinates of the triangular surface elements into a radar irradiation coordinate system in parallel by using a CUDA (compute unified device architecture), so as to obtain coordinates of the vertexes of the triangular surface elements in the radar irradiation coordinate system.
And 2, according to the coordinates of the vertexes of the triangular surface elements in a radar coordinate system, projecting the triangular surface elements to an x 'Oy' plane of the radar coordinate system by using a CUDA (compute unified device architecture), determining the number of columns and the number of rows of grids required by the triangular surface elements projected in the x 'Oy' plane and the lengths of the single grids on an x 'axis and a y' axis, further setting a corresponding number of grids in the x 'Oy' plane, and numbering the grids respectively.
And 3, calculating the depth value of each triangular surface element and the number of the grid covered by each triangular surface element in parallel by using the CUDA, wherein the plurality of triangular surface elements with the consistent projection directions cover the same grid.
For each grid in the x 'Oy' plane, respectively determining a triangular surface element with the largest depth value in triangular surface elements covering the grid, and determining the triangular surface element with the largest depth value as a bright area surface element corresponding to the grid.
And 4, utilizing the CUDA to calculate the RCS values of the bright area elements corresponding to the grids in parallel, and accumulating the RCS values of all the bright area elements to obtain the RCS of the radar target. The method for calculating the RCS of the target based on the CUDA comprises the steps of firstly dividing the surface of a three-dimensional model of a radar target into a plurality of triangular surface elements, then calculating the coordinates of the triangular surface elements in a radar irradiation coordinate system in parallel by using the CUDA, then projecting the triangular surface elements to a rectangular area in an x 'Oy' plane in the radar coordinate system, and setting a grid matrix in the rectangular area; then, the depth values of the triangular surface elements and the numbers of the grids covered by the triangular surface elements are calculated in parallel by using a CUDA (compute unified device architecture), for each grid, the triangular surface element with the largest depth value in the triangular surface elements covering the grid is respectively determined, and the triangular surface element with the largest depth value is determined as a bright area surface element corresponding to the grid; and finally, parallel computing the RCS values of the bright area elements corresponding to the grids by using the CUDA, and accumulating the RCS values of all the bright area elements to obtain the RCS of the radar target.
According to the method, when the bright area surface element is obtained, the triangular surface element is firstly placed into the corresponding grid, so that when the surface element shielding judgment is carried out, only different triangular surface elements corresponding to the same grid need to be compared, each surface element does not need to be compared with all other surface elements, the calculation complexity is reduced, and the calculation time is shortened. In addition, when the triangular surface element is placed into the corresponding grid, the method adopts a grid data structure without lock, thereby overcoming the defect of longer waiting time caused by using the grid data structure with lock. Therefore, the method has the advantages of low calculation complexity and short calculation time.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly described below, it is obvious that the drawings in the following description are only some embodiments of the present invention, and for those skilled in the art, other drawings can be obtained according to the drawings without creative efforts.
Fig. 1 is a schematic flowchart of a target RCS calculation method based on CUDA according to an embodiment of the present invention;
FIG. 2 is a comparison graph of the results of calculating RCS of a smooth vertebral body model by using FEKO software and the method of the present invention.
Detailed Description
The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
Fig. 1 is a schematic flow chart of a target RCS calculation method based on CUDA according to an embodiment of the present invention.
Referring to fig. 1, the method for calculating the target RCS based on the CUDA according to the present invention includes the following steps:
step 1, dividing the outer surface of a radar target three-dimensional model into a plurality of triangular surface elements, and converting vertex coordinates of the triangular surface elements into a radar irradiation coordinate system in parallel by using a CUDA (compute unified device architecture), so as to obtain coordinates of the vertexes of the triangular surface elements in the radar irradiation coordinate system.
It should be noted that the CUDA is used as a general computing architecture of a Graphics Processor (GPU), and includes a CUDA instruction set architecture and a parallel computing engine, and has a stronger computing capability in Processing a large amount of homogeneous data than a Central Processing Unit (CPU).
Further, step 1 comprises the following steps:
step 1.1, the surface of the three-dimensional model of the target is divided into K triangular surface elements which are numbered as 1, 2 and … … K in sequence, and a corresponding thread is allocated to each surface element.
All K threads execute the following steps 1.2 to 1.4 in parallel:
step 1.2, the kth thread in all the K threads obtains the serial numberCoordinate (x) of three vertexes of the k triangular surface element in a three-dimensional world coordinate system1 k，y1 k，z1 k)，(x2 k，y2 k，z2 k)，(x3 k，y3 k，z3 k) (ii) a Wherein K is in the form of {1, 2.
Step 1.3, calculating the rotation axis vector
Wherein,theta is the pitch angle of the radar target relative to the radar,is the azimuth angle of the radar target relative to the radar;representing a vectorSum vectorThe outer product of (d); omega1、Ω2、Ω3Respectively represent the rotation axis vectorsThe first element, the second element, and the third element.
Step 1.4, calculating a transformation matrix M according to a Rodrigue rotation formula by using the rotation axis vector:
wherein I isAn identity matrix of order 3X 3, β beingTo knowAngle of (Q)JIs a vectorThe cross-multiplication matrix of (a) is,
it should be noted that the axial vector is easily understood by those skilled in the artThe cross-multiplication matrix of (a) means: order toWill be provided withWritten in vector formTransformed form to obtainThen omegaJIs called asCross-product matrix of (a).
Step 1.5, multiplying the coordinates of the three vertexes of the triangular surface element with the number of K in the three-dimensional world coordinate system by the projection matrix M through the kth thread in all the K threads to obtain the coordinates of the three vertexes of the triangular surface element with the number of K in the radar irradiation coordinate system Wherein,
and 2, according to the coordinates of the vertexes of the triangular surface elements in the radar coordinate system, projecting the triangular surface elements to an x 'Oy' plane in the radar coordinate system by using a CUDA (compute unified device architecture), determining a rectangular area covering the projection of all the triangular surface elements in the x 'Oy' plane, and setting a grid matrix in the rectangular area.
Further, step 2 comprises the following steps:
all K threads execute the following steps 2.1 to 2.4 in parallel:
step 2.1, the kth thread in all the K threads projects the triangular surface element with the number of K to an x ' Oy ' plane, and the surface element length D of the triangular surface element with the number of K on an x ' axis is calculatedx kAnd bin length D in the y' axisy k。
Where max (. cndot.) represents taking the maximum value,
step 2.2, the K-th thread in all the K threads calculates the maximum value max of the triangular surface element with the number of K on the x' axisx kAnd minimum min on the x' axisx kMaximum value max on the y' axisy kAnd minimum min on the y' axisy k。
Wherein，
Step 2.3, calculate the grid length in the x' axisLength in y' axisWherein factor represents a grid precision parameter, Dx＝max(Dx 1，Dx 2，...，Dx K)，Dy＝max(Dy 1，Dy 2，...，Dy K)。
Step 2.4, according to the length R of the grid on the x' axisxAnd the length R of the grid in the y' axisyDetermining the number of columns of the grid required to cover the projected triangular bin in the x 'Oy' planeAnd number of linesTherein, maxx＝max(maxx 1，maxx 2，...，maxx K)，minx＝min(minx 1，minx 2，...，minx K)，maxy＝max(maxy 1，maxy 2，...，maxy K)，miny＝min(miny 1，miny 2，...，miny K)，Indicating a rounding down.
Step 2.5, according to the number of columns and rows of the grid and the length R of the single grid in the x 'axis and the y' axisyFour coordinate points (min) in the x 'Oy' planex，miny)、(minx，maxy)、(maxx，miny)、(maxx，maxy) Setting M × N grids in the rectangular region, and calculating from the coordinate point (min)x，miny) The grids are numbered 1, 2 … … mxn in order in the positive x' direction.
And 3.1, calculating the depth value of each triangular surface element and the number of the grid covered by each triangular surface element in parallel by using the CUDA, wherein the plurality of triangular surface elements with consistent projection directions cover the same grid. For each grid in the x 'Oy' plane, respectively determining a triangular surface element with the largest depth value in triangular surface elements covering the grid, and determining the triangular surface element with the largest depth value as a bright area surface element corresponding to the grid.
Further, step 3 comprises the following steps:
all K threads execute the following steps 3.1 to 3.5 in parallel:
step 3.1, the kth thread in all K threads calculates the grid number ID covered by the triangular bin with the number of Kras kAnd Depth value Depthtr ik：
Wherein, IDras k∈{1，2，...，M×N}。
the term "projection directions are the same" means that there is a case where a bright area bin and a dark area bin are projected on the same grid in the radar irradiation coordinate system, and specifically, depth values of all bins projected on the same grid are calculated, only the bin with the largest depth value among all bins is irradiated by a radar wave, and the remaining bins are blocked by the bin, and when a target RCS is calculated, only the bin irradiated by the radar wave, that is, the bright area bin contributes to the target RCS, and the blocked bin does not participate in the calculation.
Step 3.2, apply for a pointer set P in the video memoryBAnd for storing a set of pointers PBOf a contiguous memory space B, set of pointers PBContains M × N + K pointers, the pointer values are p +1, p +2 … … p + M × N + K.
Wherein, the pointer set PBPointers with medium pointer values P +1 to P + MxN point to the raster data structure pointer P1，P2，...，PM×NEach grid data structure pointer corresponds to a grid and points to the grid data structure of the corresponding grid, the grid data structure comprises the surface element number and the depth value of a triangular surface element, and the initial values are 0 and 1.7 multiplied by 10 respectively-308. Set of pointers PBThe K pointers with middle pointer values from p + mxn +1 to p + mxn + K are respectively in one-to-one correspondence with the K triangular bins, each pointer points to a corresponding grid data structure, and the data in the grid data structure is the data of the triangular bin corresponding to the pointer. The pointer with the pointer value p + M multiplied by N + k corresponds to the triangular surface element with the number of k, the triangular surface element with the number of k is contained in the pointed raster data structure, and the Depth value Depth istri kDepth values for the triangular bin numbered k.
To be explainedThe grid data structure comprises a bin number of a triangular bin and a corresponding depth value, the data type of the bin number of the triangular bin is 32-bit unsigned integer, the data type of the corresponding depth value is double-precision floating point, wherein the grid data structure corresponding to the grid is assigned an initial value, the triangular bin number is the minimum value 0 of the 32-bit unsigned integer numerical value, and the corresponding depth value is the minimum value 1.7 multiplied by 10 of the double-precision floating point numerical value-308。
Step 3.3, judging the pointer value to beWhether the raster data structure corresponding to the pointer of (1) is being updated: if not, executing the step (3.4); if yes, after the raster data structure is updated, the step (3.4) is executed.
Step 3.4, the kth thread judges the Depth value Depth in the raster data structure corresponding to the pointer with the pointer value of p + M multiplied by N + ktri kWhether greater than the pointer value ofIf so, the pointer value is set to be the depth value in the grid data structure corresponding to the pointer in the step (b), and if so, the atomic CAS function of the CUDA is used to set the pointer value to beThe depth value in the raster data structure corresponding to the pointer is updated to the depth value in the raster data structure corresponding to the pointer with the pointer value of p + M multiplied by N + kThe bin number in the grid data structure corresponding to the pointer of (1) is updated to the bin number k in the grid data structure corresponding to the pointer with the pointer value of p + M × N + k.
And 3.5, recording all surface element numbers in the grid data corresponding to the pointers with the pointer values from p +1 to p + M multiplied by N, and determining the bright area elements corresponding to the grids.
It should be noted that, as will be readily understood by those skilled in the art, if a grid data structure is directly stored in the grid, a triangular bin number and corresponding depth value are included in the grid data structure. Because every grid corresponds a thread in the CUDA, the thread will be visited in an exclusivity form when visiting the grid data structure in the video memory, guarantee not disturbed by other threads, so will lock every grid data structure, but because the locking brings the time overhead such as locking, unblock, wait, so it is not efficient to have a grid surface element data structure of lock. The invention innovatively uses a lock-free grid data structure, the pointer in the grid data structure is stored in the grid, and the length of the pointer is smaller than that of the whole grid data structure, so that the operation can be completed by using an atomic CAS function, thereby saving the operation time and improving the overall operation efficiency.
And 4, utilizing the CUDA to calculate the RCS values of the bright area elements corresponding to the grids in parallel, and accumulating the RCS values of all the bright area elements to obtain the RCS of the radar target.
Further, step 4 comprises the following steps:
all threads corresponding to all bright area elements execute the following steps 4.1 to 4.2 in parallel:
step 4.1, calculating the RCS value of the bright area element corresponding to each grid:
wherein,is the square root of the RCS value of the bright area bin corresponding to the grid numbered j,is the normal vector of the bright area bin corresponding to the grid with grid number j,is a unit vector of the electric polarization direction of the far-field receiving device,is the incident field of the radar onto the phantom,is the direction vector of the reflected wave,is an orientation vector describing the M-th edge on the bright area element corresponding to the grid with grid number j, M is the number of edges of the bright area element corresponding to the grid with grid number j, M is 3, TjIs thatThe projection length on the plane where the bright area element corresponding to the grid with the grid number j is located, k is the free space wave number, k is 2 pi f/c, c is the light speed, f is the radar irradiation frequency,is the position vector of the centroid of the bright area bin corresponding to the grid with grid number j,
step 4.2, acquiring the thread number corresponding to the bright area element corresponding to each grid, and taking a 32-modulus mode for the thread number in parallel to acquire a thread bundle number I; adding the RCS values calculated by the threads with the consistent thread bundle numbers by using a atomic Add function to obtain added data RCSISaid RCSIFrom the GPU to the CPU, where I e {1, 2, …, 32 }.
Step 4.3, in the CPU, the added data RCSIAnd adding to obtain the RCS value of the calculated radar target:
preferably, when the radar irradiation angle and the radar irradiation frequency of the irradiated target change simultaneously, the formula in step 4.1 may be used to calculate the RCS values corresponding to all frequencies to be calculated at the same angle, that is, the RCS values of all frequencies corresponding to all surface elements, then the RCS value corresponding to each frequency to be calculated at the current radar irradiation angle is calculated in sequence, then the radar irradiation angle is changed, and the above procedure is repeatedly executed, so that all RCS values of the target at different radar irradiation angles and different radar irradiation frequencies can be obtained. Because the shielding relation is consistent under the same angle and different frequencies, namely the illuminated surface elements under the same angle are the same, the repeated acquisition of bright area surface elements can be avoided, and the calculation time is reduced.
According to the method, when the bright area surface element is obtained, the triangular surface element is firstly placed into the corresponding grid, so that when the surface element shielding judgment is carried out, only different triangular surface elements corresponding to the same grid need to be compared, each surface element does not need to be compared with all other surface elements, the calculation complexity is reduced, and the calculation time is shortened. In addition, when the triangular surface element is placed into the corresponding grid, the method adopts a grid data structure without lock, thereby overcoming the defect of longer waiting time caused by using the grid data structure with lock. Therefore, the method has the advantages of low calculation complexity and short calculation time.
Further, the beneficial effects of the method of the embodiment of the invention are verified by simulation experiments as follows:
the servers with CPUs configured as Zhiqiang E5-2863V 3 and GPUs as GTX1080 respectively use the software and FEKO software models of the method to perform calculation.
Experiment I,
Simulation conditions are as follows: and the irradiation frequency f of the S-band radar is equal to GHz frequency, the azimuth angle is 0 degree, the pitch angle is 0 degree to 180 degrees at an interval of 1 degree, and horizontal polarization is realized.
Simulation content and result analysis: the software and the professional electromagnetic calculation software FEKO which use the method of the invention calculate the smooth cone model, the calculation time of the software which uses the method of the invention is 300 milliseconds, the calculation time of the professional electromagnetic calculation software FEKO is 3600 milliseconds, and the calculation result is shown in figure 2.
Experiment two,
Simulation conditions are as follows: the irradiation frequency of the L-band radar wave is 1.1GHz to 1.2GHz, 64 frequencies are taken at equal intervals, the azimuth angle is 0 degree, the pitch angle is 0 degree to 180 degrees, 1 degree is taken at intervals, horizontal polarization is carried out, and the number of surface elements is 12800.
Simulation content and result analysis: the software adopting the method and the professional electromagnetic calculation software FEKO calculate the smooth cone model, the calculation time of the software adopting the method is 399 milliseconds, and the calculation time of the professional electromagnetic calculation software FEKO is 57 minutes, so that the calculation time can be obviously reduced.
Those of ordinary skill in the art will understand that: all or part of the steps for implementing the method embodiments may be implemented by hardware related to program instructions, and the program may be stored in a computer readable storage medium, and when executed, the program performs the steps including the method embodiments; and the aforementioned storage medium includes: various media that can store program codes, such as ROM, RAM, magnetic or optical disks.
The above description is only for the specific embodiments of the present invention, but the scope of the present invention is not limited thereto, and any person skilled in the art can easily conceive of the changes or substitutions within the technical scope of the present invention, and all the changes or substitutions should be covered within the scope of the present invention. Therefore, the protection scope of the present invention shall be subject to the protection scope of the appended claims.
Claims (5)
1. A RCS (radar cross section) calculation method for target radar reflection cross section based on a CUDA (compute unified device architecture) is characterized by comprising the following steps of:
step 1, dividing the outer surface of a radar target three-dimensional model into a plurality of triangular surface elements, and converting vertex coordinates of each triangular surface element into a radar irradiation coordinate system in parallel by using a CUDA (compute unified device architecture), so as to obtain coordinates of the vertexes of each triangular surface element in the radar irradiation coordinate system;
step 2, according to the coordinates of the vertexes of the triangular surface elements in a radar coordinate system, projecting the triangular surface elements to an x 'Oy' plane of the radar coordinate system by using a CUDA (compute unified device architecture), determining the number of columns and the number of rows of grids required by the triangular surface elements projected in the x 'Oy' plane and the lengths of the single grids on an x 'axis and a y' axis, further setting a corresponding number of grids in the x 'Oy' plane, and numbering the grids respectively;
step 3, utilizing the CUDA to calculate the depth value of each triangular surface element and the number of the grid covered by each triangular surface element in parallel, wherein a plurality of triangular surface elements with consistent projection directions cover the same grid;
for each grid in an x 'Oy' plane, respectively determining a triangular surface element with the largest depth value in triangular surface elements covering the grid, and determining the triangular surface element with the largest depth value as a bright area surface element corresponding to the grid;
and 4, utilizing the CUDA to calculate the RCS values of the bright area elements corresponding to the grids in parallel, and accumulating the RCS values of all the bright area elements to obtain the RCS of the radar target.
2. The method according to claim 1, wherein the step 1 comprises the steps of:
step 1.1, dividing the surface of a three-dimensional model of a target into K triangular surface elements, sequentially numbering the surface elements as 1, 2 and … … K, and simultaneously allocating a corresponding thread to each surface element;
all K threads execute the following steps 1.2 to 1.4 in parallel:
step 1.2, acquiring coordinates (x) of three vertexes of the triangular surface element with the number of K in a three-dimensional world coordinate system by the kth thread in all the K threads1 k，y1 k，z1 k)，(x2 k，y2 k，z2 k)，(x3 k，y3 k，z3 k) (ii) a Wherein K belongs to {1, 2.., K };
step 1.3, calculating the rotation axis vector
Wherein,theta is the pitch angle of the radar target relative to the radar,is the azimuth angle of the radar target relative to the radar;representing a vectorSum vectorThe outer product of (d); omega1、Ω2、Ω3Respectively represent the rotation axis vectorsThe first element, the second element, and the third element;
step 1.4, calculating a transformation matrix M according to a Rodrigue rotation formula by using the rotation axis vector:
where I is an identity matrix of order 3X 3, β isAndangle of (Q)JIs a vectorThe cross-multiplication matrix of (a) is,
step 1.5, multiplying the coordinates of the three vertexes of the triangular surface element with the number of K in the three-dimensional world coordinate system by the projection matrix M through the kth thread in all the K threads to obtain the coordinates of the three vertexes of the triangular surface element with the number of K in the radar irradiation coordinate system Wherein,
3. the method according to claim 2, wherein the step 2 comprises the steps of:
all K threads execute the following steps 2.1 to 2.4 in parallel:
step 2.1, the kth thread in all the K threads projects the triangular surface element with the number of K to an x ' Oy ' plane, and the surface element length D of the triangular surface element with the number of K on an x ' axis is calculatedx kAnd bin length D in the y' axisy k；
Where max (. cndot.) represents taking the maximum value,
step 2.2, calculating the maximum value max of the triangular surface element with the number of K on the x' axis by the kth thread in all the K threadsx kAnd minimum min on the x' axisx kMaximum value max on the y' axisy kAnd minimum min on the y' axisy k；
Wherein,
step 2.3, determining the length of the individual grids in the x' axisAnd length in the y' axis
Wherein factor represents a grid precision parameter, Dx＝max(Dx 1，Dx 2，...，Dx K)，Dy＝max(Dy 1，Dy 2，...，Dy K)；
Step 2.4, according to the length R of the grid on the x' axisxAnd the length R of the grid in the y' axisyDetermining the number of columns of the grid required to cover the projected triangular bin in the x 'Oy' planeAnd number of lines
Therein, maxx＝max(maxx 1，maxx 2，...，maxx K)，minx＝min(minx 1，minx 2，...，minx K)，maxy＝max(maxy 1，maxy 2，...，maxy K)，miny＝min(miny 1，miny 2，...，miny K)；Represents rounding down;
step 2.5, according to the number of columns and rows of the grid and the length R of the single grid in the x 'axis and the y' axisyFour coordinate points (min) in the x 'Oy' planex，miny)、(minx，maxy)、(maxx，miny)、(maxx，maxy) Setting M × N grids in the rectangular region, and calculating from the coordinate point (min)x，miny) The grids are numbered 1, 2 … … mxn in order in the positive x' direction.
4. The method according to claim 3, wherein the step 3 comprises the steps of:
all K threads execute the following steps 3.1 to 3.5 in parallel:
step 3.1, the kth thread in all K threads calculates the grid number ID covered by the triangular bin with the number of Kras kAnd Depth value Depthtri k：
Wherein,
IDras k∈{1，2，...，M×N}；
step 3.2, apply for a pointer set P in the video memoryBAnd for storing said set of pointers PBThe set of pointers P, the contiguous memory space B ofBThe system comprises M multiplied by N + K pointers, wherein the pointer values are p +1, p +2 … … p + M multiplied by N + K;
wherein the set of pointers PBPointers with medium pointer values P +1 to P + MxN point to the raster data structure pointer P1，P2，...，PM×NEach grid data structure pointer corresponds to a grid and points to the grid data structure of the corresponding grid, the grid data structure comprises the surface element number and the depth value of a triangular surface element, and the initial values are 0 and 1.7 multiplied by 10 respectively-308(ii) a The set of pointers PBK pointers with middle pointer values from p + MxN +1 to p + MxN + K are respectively in one-to-one correspondence with K triangular surface elements, each pointer points to a corresponding grid data structure, and data in the grid data structure is data of the triangular surface element corresponding to the pointer; the pointer with the pointer value p + M multiplied by N + k corresponds to the triangular surface element with the number of k, the triangular surface element with the number of k is contained in the pointed raster data structure, and the Depth value Depth istri kThe depth value of the triangular surface element with the number of k is obtained;
step 3.3, judging the pointer value to beWhether the raster data structure corresponding to the pointer of (1) is being updated: if not, executing the step (3.4); if yes, after the raster data structure is completely updated, executing the step (3.4);
step 3.4, the kth thread judges that the Depth value Depth in the raster data structure corresponding to the pointer with the pointer value of p + M multiplied by N + ktri kWhether greater than the pointer value ofIf so, the pointer value is set to be the depth value in the grid data structure corresponding to the pointer in the set, and if so, the atomic CAS function of the CUDA is used to set the pointer value to be the depth valueThe depth value in the raster data structure corresponding to the pointer is updated to the depth value in the raster data structure corresponding to the pointer with the pointer value of p + M multiplied by N + kUpdating the bin number in the raster data structure corresponding to the pointer to be the bin number k in the raster data structure corresponding to the pointer with the pointer value of p + M multiplied by N + k;
and 3.5, recording all surface element numbers in the grid data nodes corresponding to the pointers with the pointer values from p +1 to p + MXN, and determining the bright area surface elements corresponding to the grids.
5. The method according to claim 4, wherein the step 4 specifically comprises:
all threads corresponding to all bright area elements execute the following steps 4.1 to 4.2 in parallel:
step 4.1, calculating the RCS value of the bright area element corresponding to each grid:
wherein,of a bright area bin corresponding to the grid numbered jThe square root of the value of the RCS,is the normal vector of the bright area bin corresponding to the grid with grid number j,is a unit vector of the electric polarization direction of the far-field receiving device,is the incident field of the radar onto the phantom,is the direction vector of the reflected wave,is an orientation vector describing the M-th edge on the bright area element corresponding to the grid with grid number j, M is the number of edges of the bright area element corresponding to the grid with grid number j, M is 3, TjIs thatThe projection length on the plane where the bright area element corresponding to the grid with the grid number j is located, k is the free space wave number, k is 2 pi f/c, c is the light speed, f is the radar irradiation frequency,is the position vector of the centroid of the bright area bin corresponding to the grid with grid number j,
step 4.2, acquiring the thread number corresponding to the bright area element corresponding to each grid, and taking a 32-modulus mode for the thread number in parallel to acquire a thread bundle number I; RCS value calculated by utilizing atomic Add function to lead threads with consistent thread bundle numbersAdding to obtain added data RCSISaid RCSICopying from GPU to CPU; wherein I ∈ {1, 2, …, 32 };
step 4.3, in the CPU, the added data RCSIAnd adding to obtain the RCS value of the calculated radar target:
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811389168.0A CN109558568B (en) | 2018-11-21 | 2018-11-21 | Target RCS calculation method based on CUDA |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811389168.0A CN109558568B (en) | 2018-11-21 | 2018-11-21 | Target RCS calculation method based on CUDA |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109558568A true CN109558568A (en) | 2019-04-02 |
| CN109558568B CN109558568B (en) | 2022-09-16 |
Family
ID=65866871
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811389168.0A Active CN109558568B (en) | 2018-11-21 | 2018-11-21 | Target RCS calculation method based on CUDA |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109558568B (en) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112083415A (en) * | 2020-10-12 | 2020-12-15 | 吉林大学 | Millimeter wave radar model target visibility judgment method based on 3D information |
| CN114677494A (en) * | 2022-05-26 | 2022-06-28 | 中国人民解放军国防科技大学 | Method, device and equipment for calculating radar detection capability based on subdivision grids |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103576137A (en) * | 2013-09-27 | 2014-02-12 | 电子科技大学 | Multi-sensor multi-target location method based on imaging strategies |
| CN105243280A (en) * | 2015-10-30 | 2016-01-13 | 西安电子科技大学 | Time domain physical optics algorithm based on CPU (Central Processing Unit) and GPU (Graphics Processing Unit) hybrid asynchronous parallel way |
| JP2016095764A (en) * | 2014-11-17 | 2016-05-26 | 三菱電機株式会社 | Parallel processing device and parallel processing method |
| CN108446430A (en) * | 2018-02-05 | 2018-08-24 | 西安电子科技大学 | High-frequency electromagnetic shadowing method based on sciagraphy |
| CN108646228A (en) * | 2018-05-11 | 2018-10-12 | 电子科技大学 | A kind of Radar Cross Section Calculating based on PO methods |
| CN108711335A (en) * | 2018-06-15 | 2018-10-26 | 青岛擎鹰信息科技有限责任公司 | A kind of distributed large scene radar imagery emulation mode of realization and its system |
-
2018
- 2018-11-21 CN CN201811389168.0A patent/CN109558568B/en active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103576137A (en) * | 2013-09-27 | 2014-02-12 | 电子科技大学 | Multi-sensor multi-target location method based on imaging strategies |
| JP2016095764A (en) * | 2014-11-17 | 2016-05-26 | 三菱電機株式会社 | Parallel processing device and parallel processing method |
| CN105243280A (en) * | 2015-10-30 | 2016-01-13 | 西安电子科技大学 | Time domain physical optics algorithm based on CPU (Central Processing Unit) and GPU (Graphics Processing Unit) hybrid asynchronous parallel way |
| CN108446430A (en) * | 2018-02-05 | 2018-08-24 | 西安电子科技大学 | High-frequency electromagnetic shadowing method based on sciagraphy |
| CN108646228A (en) * | 2018-05-11 | 2018-10-12 | 电子科技大学 | A kind of Radar Cross Section Calculating based on PO methods |
| CN108711335A (en) * | 2018-06-15 | 2018-10-26 | 青岛擎鹰信息科技有限责任公司 | A kind of distributed large scene radar imagery emulation mode of realization and its system |
Non-Patent Citations (5)
| Title |
|---|
| C.PANG 等: "PARALLEL IMPLEMENTATION OF PHYSICAL OPTICS BY USING GPU WITH CUDA", 《IET INTERNATIONAL RADAR CONFERENCE 2013》 * |
| LONGXIANG LINGHU 等: "Parallel Computation of EM Backscattering from Large Three-Dimensional Sea Surface with CUDA", 《SENSORS》 * |
| 孟肖等: "基于并行双尺度射线追踪的海面电磁散射计算", 《电波科学学报》 * |
| 李蔚清等: "一种基于GPU的复杂目标电磁散射快速算法", 《系统仿真学报》 * |
| 翁寅侃等: "SAR辐射定标中角反射器RCS的快速求解", 《武汉大学学报(信息科学版)》 * |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112083415A (en) * | 2020-10-12 | 2020-12-15 | 吉林大学 | Millimeter wave radar model target visibility judgment method based on 3D information |
| CN114677494A (en) * | 2022-05-26 | 2022-06-28 | 中国人民解放军国防科技大学 | Method, device and equipment for calculating radar detection capability based on subdivision grids |
| CN114677494B (en) * | 2022-05-26 | 2022-09-16 | 中国人民解放军国防科技大学 | Method, device and equipment for calculating radar detection capability based on subdivision grids |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109558568B (en) | 2022-09-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US9569559B2 (en) | Beam tracing | |
| Tao et al. | GPU-based shooting and bouncing ray method for fast RCS prediction | |
| CN102542035B (en) | Polygonal rasterisation parallel conversion method based on scanning line method | |
| US20180182158A1 (en) | Beam tracing | |
| US10529119B2 (en) | Fast rendering of quadrics and marking of silhouettes thereof | |
| CN109558568B (en) | Target RCS calculation method based on CUDA | |
| CN111475979A (en) | Acoustic target intensity simulation method based on multi-GPU multi-resolution bouncing rays | |
| Gao et al. | Mapping the SBR and TW-ILDCs to heterogeneous CPU-GPU architecture for fast computation of electromagnetic scattering | |
| WO2021063169A1 (en) | Method for rendering on basis of hemispherical orthogonal function | |
| CN105701284B (en) | The parallel calculating method of time-varying multiple dimensioned TV university region sea electromagnetic scattering vector field | |
| Johnson et al. | The irregular z-buffer and its application to shadow mapping | |
| CN110083904B (en) | Quantum radar scattering cross section calculation method based on GPU acceleration | |
| Stojanović et al. | Performance improvement of viewshed analysis using GPU | |
| CN117350083B (en) | Method and device for calculating electromagnetic scattering characteristics of super-electric large-size structure | |
| CN101430376A (en) | Target radar scattering cross-section pre-estimation system with graphics electromagnetic computation accelerated by index information | |
| CN116842695A (en) | SAR echo rapid simulation method and system based on heterogeneous computing platform CUDA | |
| JP6482193B2 (en) | Video processing apparatus and method | |
| CN107301398B (en) | A kind of identification method of image target of synthetic aperture radar realized based on GPU | |
| Zhou et al. | A Parallel Scheme for Large‐scale Polygon Rasterization on CUDA‐enabled GPUs | |
| CN113075659B (en) | Self-adaptive partitioning preprocessing method and system for grid model of flat scene | |
| CN117218259A (en) | Sample generation method based on differentiable SAR image renderer | |
| Meng et al. | An Improved SBR-PTD Method for EM Scattering from Moving Target | |
| Zhang et al. | Bistatic radar cross section prediction of 3-D target based on GPU-FDTD method | |
| Xiao et al. | Shell: A spatial decomposition data structure for ray traversal on GPU | |
| CN111598991A (en) | Computer-based method for drawing multi-thread parallel unstructured grid volume |
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 |