Summary of the invention
Based on this, be necessary the defect for above-mentioned pulsed radar moving target detection technique existence, provide a kind of degree of accuracy the high little a kind of pulsed radar moving target detection method of operand simultaneously.
A pulsed radar moving target detection method, comprising: radar system parameter is set to obtain clutter map size; According to described radar system calculation of parameter threshold value; Described clutter map is carried out to two-dimensional dct transform; According to described threshold value, the matrix data through two-dimensional dct transform is carried out to filtering; To carrying out idct transform through filtered matrix data, reduction clutter map; And the clutter map after described reduction is detected.
In the present embodiment, wherein, radar system parameter is set to obtain clutter map, comprises the steps: according to radar performance and sensing range selected frame size; According to the motion feature design sampling system of object to be detected; And the sampling precision needing according to described sampling system, obtain clutter map size.
In the present embodiment, described threshold value comprises: longitudinal low frequency threshold value Tl, longitudinal high frequency threshold value Th, horizontal low frequency threshold value Tsl and horizontal high frequency threshold value Tsh; Wherein, the computing method of described longitudinal low frequency threshold value Tl are: Tl=f
low* M/ (f
v/ 2)
In formula, f
lowfor object motion frequency to be detected, the line number that M is described clutter map, f
vfor slow time frequency acquisition; Wherein, the computing method of described longitudinal high frequency threshold value Th are:
Th=f
hi*M/(f
V/2)
In formula, f
hifor object motion frequency to be detected, the line number that M is described clutter map, f
vfor slow time frequency acquisition; Wherein, the computing method of described horizontal low frequency threshold value Tsl are:
Tsl=f
sl*N/(f
H/2)
In formula, f
slfor real pulse shape, the columns that N is described clutter map, f
hfor fast time frequency acquisition;
Wherein, the computing method of described horizontal high frequency threshold value Tsh are:
Tsh=f
sh*N/(f
H/2)
In formula, f
shfor real pulse shape, the columns that N is described clutter map, f
hfor fast time frequency acquisition.
In the present embodiment, wherein, described clutter map is carried out to two-dimensional dct transform and according to following constructive formula, realizes:
Wherein, in formula, p represents the ordinate of pixel, and q represents the horizontal ordinate of pixel
0≤p≤M-1,0≤q≤N-1, a
p, a
qbe one group of coefficient of dct transform, m represents the horizontal ordinate of original image pixel, and n represents the ordinate of original image pixel, M represents the quantity of original image pixels across point, N represents the quantity of the longitudinal pixel of original image, and Amn represents the value of the pixel that in original image, coordinate is (m, n).
In the present embodiment, wherein, according to described threshold value, the matrix data through two-dimensional dct transform is carried out to filtering and realizes according to following constructive formula:
In formula, B (i, j) represents the value of the pixel that in the transform domain image after two-dimensional dct transform, coordinate is (i, j).
In the present embodiment, wherein, to carrying out idct transform through filtered matrix data, reduction clutter map is realized according to following constructive formula:
Wherein, in formula, p represents the ordinate of transform domain pixel, and q represents the horizontal ordinate of transform domain pixel
0≤p≤M-1,0≤q≤N-1, a
pa
qit is one group of coefficient of dct transform.M represents the horizontal ordinate of original image pixel, and n represents the ordinate of original image pixel, and M represents the quantity of original image pixels across point, and N represents the quantity of the longitudinal pixel of original image, and Amn represents the value of the pixel that in original image, coordinate is (m, n).
In the present embodiment, wherein, the clutter map after described reduction is detected to the target function value that comprises the steps: to calculate each point listing on the clutter map after reduction; Record surpasses the pixel position corresponding to target function value of amplitude threshold; Add up above-mentioned pixel position, and calculate its probability density; And select the point that all density is large, be target location.
In the present embodiment, described target function value is variance yields, and the computing method of described variance yields are:
In formula, x
i(t) represent the value of the file pixel on the rear image of i.e. reduction on respective distances picture, mean (x) represents the average of all samples of this file pixel.
Described target function value is the kurtosis value of each row, and the computing method of described kurtosis value are:
Wherein: T is the sampling period, and μ gets 60*T, x conventionally
i(t) represent the value of the file pixel on the rear image of i.e. reduction on respective distances picture,
represent kurtosis.
Above-mentioned a kind of pulsed radar moving target detection method is by carrying out two-dimensional dct transform to clutter map, and dct transform calculated amount is less, and has clear and definite physical significance, and degree of accuracy is high; The threshold filter simultaneously adopting, not for characteristic frequency, has better adaptability.By adopting the pulsed radar moving target detection method of two-dimensional dct filtering can simplify calculating, there is better adaptability simultaneously.
Embodiment
Refer to Fig. 1, the schematic flow sheet of a kind of pulsed radar moving target detection method providing for the embodiment of the present invention, comprises the steps:
Step S10: radar system parameter is set to obtain clutter map size.
Refer to Fig. 2, for what the embodiment of the present invention provided, radar system parameter is set to obtain the schematic flow sheet of clutter map, step S10 is specially:
Step S11: according to radar performance and sensing range selected frame size.
Radar performance refers to that radar itself is to the sampling rate of signal and the storage capacity to information.The minor increment resolution of supposing radar is R, and storage depth is M, and the distance of a frame data representative is R*M; If R and M are adjustable, adjust as required the value of R and M, make it meet real needs.
Step S12: according to the motion feature design sampling system of object to be detected.
Sampling system should design according to practical application.For the application of appointment, if the slow time maximum frequency deviation that the motion feature of known object to be detected causes is f, choose sample frequency and be f 3~5 times, and without pursuing accurately sampling radar intermediate frequency signal; And the performance by radar itself determines for fast time-sampling frequency, corresponding with the fastest time that radar hardware can be made a response.
Step S13: the sampling precision needing according to sampling system, obtains clutter map size.
According to step S12, obtain slow time-sampling frequency, and fast time-sampling precision is determined by practical application, considers range resolution R, the frame number num-M representing on fast time shaft, the parameters such as frequency resolution Rm on slow time shaft.If slow temporal frame per second is nr, in order to reach the frequency resolution Rm on slow time shaft, the one-tenth frame number N=nr/Rm needing; Clutter map is of a size of: M * num-M * N, and wherein M is counting of fast time one frame, N is image longitudinal size.
Step S20: according to radar system calculation of parameter threshold value.
Radar system parameter can comprise: fast time-sampling frequency is f
h, slow time-sampling frequency is f
vand the clutter map size of intercepting, be designated as the capable N row of M; Threshold value comprises: longitudinal low frequency threshold value Tl, longitudinal high frequency threshold value Th, horizontal low frequency threshold value Tsl and horizontal high frequency threshold value Tsh; Here getting object motion frequency to be detected is that 0.2Hz~5Hz is also designated as respectively f
low~f
hi, real pulse shape f
sl~f
sh:
Tl=f
low*M/(f
V/2)
Th=f
hi*M/(f
V/2)
Tsl=f
sl*N/(f
H/2)
Tsh=f
sh*N/(f
H/2)。
Step S30: clutter map is carried out to two-dimensional dct transform.
DCT is a kind of real number field conversion, and its transformation kernel is real number cosine function; Piece image is carried out after discrete cosine transform, and the important information of many relevant images all concentrates in the sub-fraction coefficient of dct transform.In embodiment provided by the invention, clutter map is carried out to two-dimensional dct transform and according to following constructive formula, realizes:
Wherein, in formula, p represents the ordinate of transform domain pixel, and q represents the horizontal ordinate of transform domain pixel
0≤p≤M-1,0≤q≤N-1,
A
p, a
qbe one group of coefficient of dct transform, m represents the horizontal ordinate of original image pixel, and n represents the ordinate of original image pixel, M represents the quantity of original image pixels across point, N represents the quantity of the longitudinal pixel of original image, and Amn represents the value of the pixel that in original image, coordinate is (m, n).
Step S40: the matrix data through two-dimensional dct transform is carried out to filtering according to threshold value.
DCT is a kind of real number field conversion, and its transformation kernel is real number cosine function.Piece image is carried out after discrete cosine transform, and the important visual information of many relevant images all concentrates in the sub-fraction coefficient of dct transform.Known according to the meaning of dct transform, the energy of whole two field picture is assembled to the upper left corner, represents that frequency is lower the closer to the upper left corner.Therefore according to the parameter of radar system, can calculate the wave filter that will design.Remember that fast time-sampling frequency is f
h, slow time-sampling frequency is f
v, the clutter map of intercepting is of a size of the capable N row of M; Here the motion frequency 0.2Hz~5Hz that gets object to be detected is designated as f
low~f
hi, real pulse shape f
sl~f
sh:
Tl=f
low*M/(f
V/2)
Th=f
hi*M/(f
V/2)
Tsl=f
sl*N/(f
H/2)
Tsh-f
sh*N/(f
H/2)
In formula, Tl is longitudinal low frequency threshold value, and Th is longitudinal high frequency threshold value, and Tsl is horizontal low frequency threshold value, and Tsh is horizontal high frequency threshold value, and B (i, j) represents the value of the pixel that in the transform domain image after two-dimensional dct transform, coordinate is (i, j).The new figure DCT coefficient of clutter map being formed according to above formula carries out filtering, specifies the DCT coefficient in area to be cleared, and has retained the data in the frequency range of our needs.Certainly if carry out more complicated filtering, only need calculate corresponding size and the position of frequency range of wanting filtering according to above formula, then corresponding DCT coefficient be carried out to zero clearing or specific modification can form new algorithm, the method adaptability is very wide.
Step S50: to carrying out idct transform through filtered matrix data, reduction clutter map.
In embodiment provided by the invention, to carrying out idct transform through filtered matrix data, reduction clutter map is realized according to following constructive formula:
Wherein, in formula, p represents the ordinate of transform domain pixel, and q represents the horizontal ordinate of transform domain pixel
0≤p≤M-1,0≤q≤N-1; a
p, a
qit is one group of coefficient of dct transform.M represents the horizontal ordinate of original image pixel, and n represents the ordinate of original image pixel, and M represents the quantity of original image pixels across point, and N represents the quantity of the longitudinal pixel of original image, and Amn represents the value of the pixel that in original image, coordinate is (m, n).
Step S50: the clutter map after reduction is detected.
Refer to Fig. 3, the schematic flow sheet that the clutter map after reduction is detected for the embodiment of the present invention provides, comprises the steps:
Step S51: the target function value that calculates each point listing on the clutter map after reduction.
For the clutter map through DCT filtering inverse transformation, its size is still N * M, the capable M row of N.The N that each is listed point carries out the value of statistical computation objective function, and the covariance of take here describes as example: establishing the pixel that i lists is x
i(t) t ∈ [1, N], the variance on this Range Profile
In formula, x
i(t) represent the value of the file pixel on the rear image of i.e. reduction on respective distances picture, mean (x) represents the average of all samples of this file pixel.
Be appreciated that except covariance, can also each kurtosis being listed as calculate as objective function.Kurtosis representative be the intensity that one group of data distributes, when data sample is very large, or when unknown, should meet normal distribution, its kurtosis is 3; If change because some reason causes the distribution of data, kurtosis must depart from 3, so we also can judge whether target by calculating its kurtosis.The calculating formula of kurtosis can adopt the iterative algorithm of simplification:
Wherein: T is the sampling period; μ gets 60*T, x conventionally
i(t) represent the value of the file pixel on the rear image of i.e. reduction on respective distances picture,
represent kurtosis.This be one iterative, by this formula traversal sample, can try to achieve numerical value.。
Step S52: record surpasses the pixel position corresponding to target function value of amplitude threshold.
For adopting variance, as the situation of objective function, choose peaked 80% left and right threshold value, calculate the variance yields on all M Range Profile, get it peaked 80% for threshold value, variance is greater than the distance picture point i of threshold value, i ∈ [1, M] is recorded in a new array Index.
For the situation that adopts kurtosis as objective function, kurtosis value and 3 is done to subtraction, then get its absolute value, select the Range Profile that is greater than maximal value 80% on all Range Profiles, indicate target appearance, and be less than 80% data, abandon.
Step S53: add up above-mentioned pixel position, and calculate its probability density.
In embodiment provided by the invention, for guaranteeing robustness, repeatedly repeat the method for above-mentioned pixel position; Multiplicity is advisable between 4~10 times.The number of times that in statistics Index, each i occurs, calculates its probability (occurrence number is divided by multiplicity).
Step S54: select the point that all density is large, be target location.
Select N the point that density is large, be target location: the probability density in Index is got to peak value, have how many peaks just to have how many moving targets, and its corresponding Range Profile i, be the distance and position of target.
The size that is appreciated that density depends on artificial selection, and in the present embodiment, being preferably more than 80% of maximum density values is large density.
By specific embodiment, to of the present invention, be further described below.
For convenience of description, suppose that upper every frame storage depth of fast time is 512 points, corresponding distance is between points about 8 millimeters; Slow time shaft adopts the sample frequency of 20Hz.We choose 200 frames is a decision window, every 10 seconds, can obtain 200 row, the clutter map of 512 row, be designated as Frame_org, wherein each pixel is designated as f (i, n), wherein i represents the position on fast time shaft, and the span of i is [1,512]; N represents the position on slow time shaft, and the span of n is [1,200].
Refer to Fig. 4, for what the present embodiment provided, to each frame, go average value processing to obtain new clutter map, be designated as Frame.This figure is the original image that radar is obtained, and transverse axis represents fast time shaft, has 512 points here on transverse axis, and on ordinate, we adopt be a window at 256, has shared 4 windows and has spliced, altogether 1024 points.
Refer to Fig. 5, for provided by the invention, clutter map is carried out to the clutter map B that two-dimensional dct transform obtains.Discrete cosine transform (DCT) is a kind of real number field conversion, and its transformation kernel is real number cosine function; Piece image is carried out after discrete cosine transform, and the important information of many relevant images all concentrates in the sub-fraction coefficient of dct transform.The two-dimensional dct of a matrix is defined as follows:
Wherein, in formula, p represents the ordinate of pixel, and q represents the horizontal ordinate of pixel
0≤p≤M-1,0≤q≤N-1, a
p, a
qbe one group of coefficient of dct transform, m represents the horizontal ordinate of original image pixel, and n represents the ordinate of original image pixel, M represents the quantity of original image pixels across point, N represents the quantity of the longitudinal pixel of original image, and Amn represents the value of the pixel that in original image, coordinate is (m, n).。
In clutter map B, horizontal ordinate and ordinate are the coefficient coordinate of original image after dct transform, are expressed as the frequency of cosine function.From clutter map, B finds out, adopts above-mentioned constructive formula (1) to carry out two-dimensional dct to Frame, obtain having onesize DCT coefficient with former clutter map, and the important information of image all concentrates on the upper left quarter of image.
After obtaining above-mentioned pattern, according to threshold value, the matrix data through two-dimensional dct transform is carried out to filtering.Specific as follows: to remember that fast time-sampling frequency is f
h, slow time-sampling frequency is f
v, the clutter map of intercepting is of a size of the capable N row of M; Here get human motion frequency 0.2Hz~5Hz and be designated as f
low~f
hi, real pulse shape f
sl~f
sh:
Tl=f
low*M/(f
V/2)
Th=f
hi*M/(f
V/2)
Tsl=f
si*N/(f
H/2)
Tsh=f
sh*N/(f
H/2)
In formula, T1 is longitudinal low frequency threshold value, and Th is longitudinal high frequency threshold value, and Tsl is horizontal low frequency threshold value, and Tsh is horizontal high frequency threshold value.The new figure B DCT coefficient of clutter map being formed according to above formula carries out filtering, specifies the DCT coefficient in area to be cleared, and has retained the data in the frequency range of our needs.
Refer to Fig. 6, for what the embodiment of the present invention provided, to carry out IDCT (DCT inverse transformation) through filtered matrix data, convert the clutter map C obtaining.
In the present embodiment, to carrying out idct transform through filtered matrix data, reduction clutter map is realized according to following constructive formula:
Wherein, in formula, p represents the ordinate of transform domain pixel, and q represents the horizontal ordinate of transform domain pixel
0≤p≤M-1,0≤q≤N-1; a
pa
qit is one group of coefficient of dct transform.M represents the horizontal ordinate of original image pixel, and n represents the ordinate of original image pixel, and M represents the quantity of original image pixels across point, and N represents the quantity of the longitudinal pixel of original image, and Amn represents the value of the pixel that in original image, coordinate is (m, n).
In clutter map C, represent the time-domain diagram after filtering reduction, from clutter map C, can find out the vestige that has motion at several Range Profiles.According to the clutter map C obtaining, each column data is added up, selected take energy as objective function judges here, calculate its covariance (or AC energy), the peak value that is greater than maximal value 80% can be described as obviously larger, will meet the point location target location of this condition.Except covariance, can also each kurtosis being listed as calculate as objective function.Kurtosis representative be the intensity that one group of data distributes, when data sample is very large, or when unknown, should meet normal distribution, its kurtosis is 3; If change because some reason causes the distribution of data, kurtosis must depart from 3, so we also can judge whether target by calculating its kurtosis.The calculating formula of kurtosis can adopt the iterative algorithm of simplification:
Wherein: T is the sampling period, and μ gets 60*T, x conventionally
i(t) represent the value of the file pixel on the rear image of i.e. reduction on respective distances picture,
represent kurtosis.
Refer to Fig. 7, the pattern that 4 windows are carried out to Conjoint Analysis providing for the embodiment of the present invention.What at Fig. 7, represent is that each Range Profile is asked to the result after variance, is specially first each window is calculated, and then the value of 4 windows is stitched together and obtains last criterion.For the sake of assurance, can repeatedly calculate a plurality of such windows, then result be carried out to Monte Carlo analysis, choose and occur the several of maximum probability, as net result.From Fig. 7, show, the method that the embodiment of the present invention provides can effectively be extracted motion and range information thereof.
Refer to Fig. 8, the final decision diagram providing for the embodiment of the present invention.In Fig. 8, horizontal ordinate represents Range Profile, and ordinate represents on corresponding Range Profile, is greater than the number of the point of maximum variance value 80%.From Fig. 5, clearly find out and have 3 moving targets, its position is respectively 0.3 meter, 1 meter, 3.4 meters, coincide with experiment situation.
The above, it is only preferred embodiment of the present invention, not the present invention is done to any pro forma restriction, although the present invention discloses as above with preferred embodiment, yet not in order to limit the present invention, any those skilled in the art, do not departing within the scope of technical solution of the present invention, when can utilizing the technology contents of above-mentioned announcement to make a little change or being modified to the equivalent embodiment of equivalent variations, in every case be not depart from technical solution of the present invention content, any simple modification of above embodiment being done according to technical spirit of the present invention, equivalent variations and modification, all still belong in the scope of technical solution of the present invention.