Summary of the invention
The present invention proposes a kind of weather radar image interpolation method based on Fourier spectrum analysis and device, can highlight because of beam-broadening, resolution decline and cause radar detection less than strong echo center, and algorithm routine is solidified into chip, make fast operation, the chip apparatus of real-time requirement can be met.
For achieving the above object, the technical solution adopted in the present invention is as follows:
Based on the air radar image interpolation method that Fourier spectrum is analyzed, comprise the following steps:
1) on given radar observation direction, with Δ t=1 storehouse for sampling interval, N number of echo strength that obtains that this side up, is designated as:
X
1,X
2,……,X
N(1)
Wherein, X
t, t=1,2 ..., N, represents the echo strength of t position, represents distance for radial data t, and for tangential data, t represents position angle;
2) according to Fourier spectrum analysis, the echo strength X of t position
talso can be expressed as:
Wherein, A
0for X
1, X
2..., X
nthe arithmetic mean of sequence, A
ibe the amplitude of the i-th subharmonic,
be the phase angle of the i-th subharmonic, ω
i=2 π f
ibe the i-th subharmonic circular frequency, t is X
tpositional information, p is harmonic component number, f
iit is the i-th subfrequency;
Formula (2) also can be expressed as:
Wherein Fourier coefficient a
iand b
ibe respectively:
By Fourier coefficient a
iand b
iobtain harmonic amplitude A
iwith phase place
3) a certain elevation angle PPI figure of document border radar has G position angle, is designated as: α
1, α
2..., α
g, each position angle has H range bin, is designated as: r
1, r
2..., r
h, therefore, the data that in a radial direction, H range bin forms are called radial data, are designated as:
the same one-dimensional data formed apart from a upper G orientation is called tangential data, is designated as:
4) utilize formula (4) and (5) tangential to the 1st on G data calculating Fourier coefficient, harmonic amplitude and phase place; 1st tangentially refers to the nearest range bin of distance radar, i.e. r
1position, in the calculation, makes the N in formula (4) equal G, X in same up-to-date style (4)
t, t=1,2 ..., N, gets r
1g tangential data of position
5) by step 4) harmonic amplitude that calculates and phase place substitute into formula (2), and make the t in formula (2) equal α, and calculate the Echo Rating at the 1st tangential upper, orientation α place, wherein, α is the position angle of interpolation point;
6) X successively in modus ponens (4)
t, t=1,2 ..., N is r
2, r
3..., r
hthe G of position tangential data, repeat step 4) and step 5), calculate the 2nd tangentially on, that is, r
2the Echo Rating at position, orientation α place, the 3rd tangential on, i.e. r
3the Echo Rating at position, orientation α place ..., until calculate all H tangential on, Echo Rating when position angle is α, using this H tangential Echo Rating as one group of radial data;
7) formula (4) and (5) are utilized, calculation procedure 6) this group Fourier coefficient of radial data, harmonic amplitude and phase place of obtaining; In the calculation, the N in formula (4) is made to equal H, X in same up-to-date style (4)
t, t=1,2 ..., N, gets step 6) obtain this group radial data;
8) by step 7) harmonic amplitude that calculates and phase place substitute into formula (2), and make the t in formula (2) equal r, calculate position angle and be α and distance is the Echo Rating at r place, namely interpolation point P (r is obtained, Echo Rating α), wherein, r is the distance of interpolation point.
The described step 8 of aforesaid change) in r, the Echo Rating that position angle is α any distance place can be obtained.
The described step 5 of aforesaid change) in α, the Echo Rating that arbitrary orientation angular distance is r place can be obtained.
Based on the device of the air radar image interpolation method that Fourier spectrum is analyzed, comprise PC and FPGA plate;
Described PC is used for pre-service and realizes the communication between FPGA plate;
Described FPGA plate is provided with control module, computing module, memory module, communication adaptation module and FIFO; Wherein, the control that described control module is responsible for the image algorithm treatment scheme of whole FPGA plate performs, and the cooperation between coordination modules is with synchronous;
Described computing module reads original image data from memory module, processes, and result is sent to main frame by communication adaptation module according to the air radar image interpolation method based on Fourier spectrum analysis;
Described memory module adopts the ram in slice of FPGA, and employing VerilogHDL language realizes 8 RAM in a sheet, can store 4096 unit; Described computing module reads data and write data in RAM from memory module RAM;
Communication protocol between described communication adaptation module and main frame adopts RS232, start bit, eight bit data position and a position of rest during transmission, and communication interface adopts the serial ports on FPGA plate;
Establish FIFO between described computing unit and communication adaptation module, FIFO is the data buffer of a first in first out.
The inventive method is when processing the radar return data that the strong echo central feature in convective region is given prominence to, the architectural feature of strong echo area can be highlighted, especially can distance radar remotely interpolation go out because of the low and undetectable strong echo center of radar resolution, for accurately identifying that strong echo center and intensity provide reference further, and can the algorithm routine of the inventive method be solidified into chip, make fast operation, the chip apparatus of real-time requirement can be met.
Embodiment
Now with embodiment, the present invention is described in further detail by reference to the accompanying drawings.
The essence of the inventive method is that spectrum signature is also used for the matching of observation data discrete point, carries out resampling again by spectrum signature further that utilize Fourier spectrum analysis to obtain echo strength change in time and space.Therefore be conducive to showing some because of the low and undetectable strong echo center of radar resolution.
Fourier spectrum analysis is a kind of method conventional in data signal analysis, utilizes frequency, amplitude and phase place to describe the vibrational waveform in time domain or spatial domain, and then analyzes the rule and characteristic of data, and is widely used in meteorology.The present invention regards radar observation data as time series, and e.g., on given radar observation direction, with Δ t=1 storehouse for sampling interval, the maximum distance apart to sampled point has T=460 storehouse, then obtain N=T/ Δ t=460 echo strength this side up, be designated as:
X
1,X
2,……,X
N(1)
Wherein, X
t(t=1,2 ..., N), represent the echo strength of t position, represent distance for radial data t, for tangential data, t represents position angle.
The basic cycle of echo strength equals the length of sequence, i.e. T=N, and ripple corresponding to basic cycle is called first-harmonic, and fundamental frequency is f
1=1/N.I-th subfrequency is designated as f
i=i*f
1=i/N.According to Fourier spectrum analysis, the echo strength X of t position
talso can be expressed as:
Wherein, A
0for X
1, X
2..., X
nthe arithmetic mean of sequence, A
ibe the amplitude of the i-th subharmonic,
be the phase angle of the i-th subharmonic, ω
i=2 π f
ibe the i-th subharmonic circular frequency, t is X
tpositional information (distance, orientation etc. can be represented), p is harmonic component number, p=N/2=230, has 230 harmonic components.
Formula (2) also can be expressed as:
Wherein Fourier coefficient a
iand b
ibe respectively:
Can be calculated by Fourier transform.
Harmonic amplitude A is obtained by Fourier coefficient
iwith phase place
Due to X
tobtain through radar observation discrete sampling with distance continually varying echo strength, therefore can be considered the change of radar echo intensity with distance, so, the harmonic amplitude that formula (5) is calculated and phase place substitute into formula (2), demand perfection portion harmonic component sum (p=N/2=230,230 harmonic components), the X corresponding to any distance can be obtained
t, namely for any distance in this radial direction, corresponding Echo Rating can be provided.In fact, the essence of carrying out " interpolation " like this utilizes Fourier spectrum analysis to obtain the spectrum signature of echo strength change in time and space and is further used for the matching of observation data discrete point, carries out resampling again.
The a certain elevation angle PPI figure of document border radar has G position angle, is designated as: α
1, α
2..., α
g, each position angle has H range bin, is designated as: r
1, r
2..., r
h.Therefore, the data that in a radial direction, H range bin forms are called radial data (H data), are designated as:
the same one-dimensional data formed apart from a upper G orientation is called tangential data (G data), is designated as:
so a radar PPI figure contains G group radial data and the tangential data of H group.
Utilize above-mentioned theory, to the tangential data of the two-dimensional radar echo data of PPI form and radial data first, after carry out one dimension interpolation arithmetic, the Echo Rating of arbitrary orientation and range bin can be obtained on PPI figure.As shown in Figure 1, according to the character of two-dimensional Fourier transform, order that is tangential, radial data interpolation does not affect result, the interpolation sequence that employing is here first tangentially, radial again.As shown in Figure 1, black color dots is the point that radar detection is arrived, and P is any interpolation point, α
1, α
2for adjacent position angle, r
1, r
2, r
3for neighbor distance storehouse, α and r is the azimuth-range of interpolation point P.
The step of the inventive method is as follows:
(1) utilize formula (4) and (5) tangential to the 1st on G data calculating Fourier coefficient, harmonic amplitude and phase place; 1st tangentially refers to the nearest range bin of distance radar, the r namely in Fig. 1
1position, in the calculation, makes the N in formula (4) equal G, X in same up-to-date style (4)
t(t=1,2 ..., N) and get r
1g tangential data of position
(2) harmonic amplitude step (1) calculated and phase place substitute into formula (2), and make the t in formula (2) equal α, calculate the Echo Rating at the 1st tangential upper, orientation α place; Wherein, α is the position angle of interpolation point;
(3) X successively in modus ponens (4)
t(t=1,2 ..., N) and be r
2, r
3..., r
hg tangential data of position, repeat step (1) and step (2), calculate the 2nd the tangentially upper (r namely in Fig. 1
2position), the Echo Rating at orientation α place, the 3rd the tangential upper (r namely in Fig. 1
3position), the Echo Rating at orientation α place ..., until calculate all H tangential on, Echo Rating when position angle is α, using this H tangential Echo Rating as one group of radial data, the dotted line ringlet namely in Fig. 1;
(4) formula (4) and (5) are utilized, this group Fourier coefficient of radial data, harmonic amplitude and phase place that calculation procedure (3) obtains; In the calculation, the N in formula (4) is made to equal H, X in same up-to-date style (4)
t(t=1,2 ..., N) and get this group radial data that step (3) obtains;
(5) harmonic amplitude step (4) calculated and phase place substitute into formula (2), and make the t in formula (2) equal r, calculate position angle and be α and distance is the Echo Rating at r place, namely interpolation point P (r is obtained, Echo Rating α), wherein, r is the distance of interpolation point.
In said process, change r, the Echo Rating that position angle is α any distance place can be obtained.
In said process, change α, the Echo Rating that arbitrary orientation angular distance is r place can be obtained.
Like this, the Echo Rating of arbitrfary point in radar coverage on PPI can be drawn.
In addition, above-mentioned algorithm routine can be solidified into chip, make fast operation, the chip apparatus of real-time requirement can be met, as follows:
Programmable gate array (FPGA) occurs as a kind of semi-custom circuit in special IC (ASIC) field, has both solved the deficiency of custom circuit, has overcome again the shortcoming that original programming device gate circuit number is limited.Therefore become can the optimum development scheme of repeated configuration circuit module for FPGA.FPGA technology and performance have reached the algorithm that can be used for realizing image procossing.So the present invention is used for FPGA technology in the process application of radar image.
Hardware platform build the Resources on Chip utilizing FPGA completely, use VerilogHDL hardware description language to achieve data communication between one and host computer.Because FPGA itself has repeated configuration completely, the hardware platform of this framework can be upgraded and repeated configuration very easily.As shown in Figure 2, PC and FPGA plate is mainly comprised.PC is used for pre-service realize the communication between FPGA plate by its communication software, and pre-service refers to the process such as the reading of radar data and quality control.
FPGA plate is provided with control module, computing module, memory module, communication adaptation module and FIFO (as shown in Fig. 2 dotted line frame).Control module CA (Controlautomation) is in top-level module, and the control being responsible for the image algorithm treatment scheme of whole FPGA plate performs, and the cooperation between coordination modules is with synchronous.
Computing module is responsible for the evaluation work of specific algorithm, reads data and write data from memory module RAM in RAM.Memory module adopts the ram in slice of FPGA, and employing VerilogHDL language realizes 8 RAM in a sheet, can store 4096 unit.Computing module reads original image data from RAM, and the algorithm proposed according to the present invention processes, and result is sent to main frame by communication adaptation module.Communication protocol between communication adaptation module and main frame adopts RS232, and start bit, eight bit data position and a position of rest during transmission, communication interface employs the simple serial ports on FPGA plate.Between computing unit and communication adaptation module, be also provided with FIFO, FIFO is the data buffer of a first in first out.
When the hardware platform of this framework realizes different image processing algorithms, what need design iterations only has computing module.
Now with a specific embodiment, the inventive method is described in further detail.
During 29 days 00 July in 2008 18 points, No. 0808 typhoon " phoenix " logs at Foochow Southern Coast, and now Foochow radar observes the echo PPI obtained scheme with 0.5 ° of elevation angle, as shown in Figure 3.The central point of Fig. 3 is radar site.Can see having in the southeastern direction of radar the contra solem banded echo district that a scope is larger, banded echo has multiple strong echo center.For strong convection region in distant square frame, as Fig. 4 (a) after amplifying, resolution is 1 ° × 1km (" 1 ° " is equivalent to 3km at distance 180km place).In order to carry out interpolation contrast, according to " space weighted average " (see in radar meteorology to the definition of radar reflectivity factor) resolution of Fig. 4 (a) is reduced to 2 ° × 2km, as Fig. 4 (b), Fig. 4 (b) strong echo center, several places because resolution reduces in level and smooth Fig. 4 (a).Suppose that Fig. 4 (b) detects the reflectogram obtained for low resolution radar, according to 1 ° × 1km resolution, interpolation is carried out to Fig. 4 (b), use bilinear interpolation and Fourier's interpolation two kinds of methods respectively, obtain result respectively as Fig. 4 (c) and 4 (d).Although Fig. 4 (c) spatial resolution improves, but the smoothing property because of bilinear interpolation has buried the strong echo center, several places of original existence further, Fig. 4 (d) is then because of the strong echo center display improved in Fig. 4 (b) with the obvious advantage of Fourier's interpolation.Fig. 4 (d) and Fig. 4 (a) closer to.Especially in the selected region A of Fig. 4 dotted line frame, have more strong echo area (in Fig. 4 (a) dark point) in high-resolution radar data, low resolution radar has only detected the strong echo of a fritter (in Fig. 4 (b) dark point).After carrying out interpolation to low resolution radar echo, bilinear interpolation (Fig. 4 (c)) does not only have interpolation to go out actual due echo, and the scope of strong echo central point have also been smaller; After Fourier's interpolation (Fig. 4 (d)), occurred some stronger echoes in region, only its position and Fig. 4 (a) differ from a pixel (i.e. 1km).Same, in the B of region, high-resolution radar has detected two strong echo central points, and low resolution radar does not detect, there is not strong echo central point in the result after bilinear interpolation, Fourier's interpolation makes two strong echo centers reappear, only alternate position spike pixel yet.There is large stretch of weak echo region in Fig. 4 lower right corner, in this sheet weak echo region, Fig. 4 (b), 4 (c) are more consistent with 4 (d), and illustrate that Fourier's interpolation is when processing weak echo region data, its effect and bilinear interpolation are more or less the same.
In order to express this advantage of Fourier's interpolation further, using the echo strength of pixels all in Fig. 4 (a) as " true value ", carry out pixel paired comparisons one by one respectively to Fig. 4 (c) and 4 (d), pixel and Fig. 4 (c) of obtaining more than true value 40dBZ are as shown in table 1 with the statistics of the respective pixel value in 4 (d).
Table 1 bilinear interpolation and Fourier's interpolation true value comparing result
In table 1, data lexical or textual analysis is with the 1st behavior example, true value is the pixel group of 40.5dBZ, 58 are had in Fig. 4 (a), the mean value of the corresponding value of this pixel group in Fig. 4 (d) (i.e. Fourier's interpolation) is 39.4dBZ, standard deviation is 1.9dBZ, and the mean value of the value of this pixel group in Fig. 4 (c) (i.e. bilinear interpolation) is 39.3dBZ, standard deviation is 2.1dBZ, this shows, for the pixel group (containing 58 pixels) that true value is 40.5dBZ, Fourier's interpolated average is closer to true value, this shows that Fourier's method of interpolation is better than bilinear interpolation.From row data each in table 1, along with true value increases (namely real echo strength strengthens), it is more obvious that Fourier's method of interpolation is better than bilinear interpolation.In table 1, last column data show, the pixel of more than true value 40dBZ has 672, average 45.5dBZ, and Fourier's interpolated average is 44.8dBZ, standard deviation is 1.6dBZ, and bilinear interpolation mean value is 44.1dBZ, standard deviation is 1.5dBZ.Therefore, can illustrate, at " stronger echo area ", the result that two kinds of method of interpolation obtain is from no significant difference viewed from standard deviation, but Fourier's interpolated average is closer to true value, again shows that Fourier's method of interpolation is better than bilinear interpolation.This judgement can also be found out from the Fig. 5 (a) made by table 1, in figure, ● represent Fourier's interpolated average, ■ represents bilinear interpolation mean value, solid line represents linear Fourier's interpolated average, dotted line represents Linear Double linear interpolation mean value, and from Fig. 5 (a), Fourier's interpolated average and true value have better consistance (degree of fitting R
2comparatively large, be 0.9839; Fit line slope closer to 1, intercept closer to 0dBZ), and consistance between bilinear interpolation mean value and true value is just weaker, and (degree of fitting is less, is 0.9432; Fit line slope is also less, is 0.7930, and intercept is comparatively large, is 8.0369dBZ).
In like manner, the pixel of true value within the scope of 30-40dBZ (i.e. medium tenacity echo area) in Fig. 4 (a) is carried out contrast with the respective pixel value in Fig. 4 (c) and 4 (d) add up and make Fitting Analysis, can Fig. 5 (b) be obtained; The pixel of true value within the scope of 10-30dBZ (i.e. weak echo region) in Fig. 4 (a) is carried out contrast with the respective pixel value in Fig. 4 (c) and 4 (d) add up and make Fitting Analysis, can obtain Fig. 5 (c) (pixel of below 10dBZ due to number less, and without diastrous weather indicative significance, so it is not contrasted).In figure, ● represent Fourier's interpolated average, ■ represents bilinear interpolation mean value, solid line represents linear Fourier's interpolated average, dotted line represents Linear Double linear interpolation mean value, from loose distribution and the fitting parameter thereof of Fig. 5 (b) and 5 (c), in medium tenacity and weak echo region, especially the medium tenacity around " stronger echo area " and weak echo region, two kinds of interpolation method difference on effect are little.
In sum, Fourier's method of interpolation that the present invention proposes, can highlight strong echo area and strong center thereof better, contributes to identifying the strong convection of its correspondence, diastrous weather position and scope.