CN106294284A - Rainfall KPT Scatter computing accelerated method based on frequency angle two-dimensional mixing interpolation - Google Patents

Abstract

The invention discloses rainfall KPT Scatter computing accelerated method based on frequency angle two-dimensional mixing interpolation, the method is combined into model with cubic spline interpolation and Stoer Bulirsch.RCS (RCS) comprises frequency and angle information simultaneously, the frequency at wave band interested repeated and angle band is needed to obtain RCS, the core of the present invention is exactly at the widest angle band application Interpolatory Splines of faradic current, at broadband band application Stoer Bulirsch.In order to effectively calculate Electromagnetic Scattering Characteristics, improve the operation efficiency of matrix-vector multiplication in improved Electric Field Integral Equation, apply sparse matrix/canonical grid (SM/ CG).The present invention saves a large amount of operation time when carrying out frequency and angle scanning, and effectively raises the arithmetic speed of rainfall population bandwidth scattering.

Description

Rainfall KPT Scatter computing accelerated method based on frequency angle two-dimensional mixing interpolation

Technical field

The invention belongs to Radar Technology field, particularly relate to a kind of rainfall particle based on frequency angle two-dimensional mixing interpolation Scattering computing accelerated method.

Background technology

The concept of radar originated from for 20 beginnings of the century, became army because having round-the-clock, the remote distant sensing function of round-the-clock The electronic equipment that thing field is indispensable.Along with developing rapidly of scientific and technological level, for different application background radar system also Appear widely in the every field of national economy.Precipitation radar is as an important branch of weather radar, to precipitation, ice The weather phenomena such as hail are monitored and forecast aspect plays an important role, and have obtained the great attention of national governments.From rainfall Em polarization characteristic angle considers, rainfall particle is the basis of rainfall detection to the scattering of radar beam.But, precipitation particles Scattering properties is the most complicated, not only relevant with the size of precipitation particles, shape and phase, is also largely dependent upon particle and exists The factors such as the orientation in space, orientation.Therefore, the scattering properties analyzing rainfall particle theoretically has important researching value.

In recent years, about in the remote sensing application of environment, to the rainfall population electromagnetic wave phase interaction with complexity Characterization research is very important.The common method solving complicated rainfall population scattering is dependent on digital technology, such as, Moment method (MOM).Moment method (MOM) can provide the precise solution that complex dielectrics scatters, but is solving integral equation Time computing extremely complex and need the biggest memory space.In order to alleviate these bottlenecks, there has been proposed some quickly, effectively Method, such as: conjugate gradient fast Fourier transform (CG-FFT) method, Fast multipole method (FMM) etc. reduce in computer Deposit the demand with amount of calculation, improve the efficiency of moment method.Zeng Liang and Chen Zhihao professors etc. are analyzing random rough surface scattering problems Time, product impedance matrix being divided into strong interaction matrix and weak interaction matrix with waiting to ask vector, greatly reduce meter Calculate complexity and storage complexity.In recent years, Fast numerical method sparse matrix/regular grid (SM/CG) method quilt Propose, be applied in the scattering analysis of three-dimensional rainfall population.

RCS (RCS) comprises frequency and angle information simultaneously.In many actual application, prediction in ideal There is frequency domain and spatial domain in mono-static RCS target simultaneously.Although computation complexity and internal memory can be reduced by SM/CG method Requirement, but we need for repeat each frequency in frequency range interested or angle calculation, so as at the widest angle Frequency domain wave band obtains RCS.In order to solve this difficulty, it is proposed that many interpolating methods improve efficiency, as Asymptotic Waveform is estimated (AWE) method and parameter estimation based on model (MBPE) method.But these methods should use the most highly difficult, AWE method uses Time, sometimes it is difficult to impedance matrix and the derivative of faradic current vector obtained.The where the shoe pinches of MBPE method is to need matrix Invert and obtain the coefficient of alternative model, cause wasting the plenty of time when calculating frequency and angle scanning.

Summary of the invention

The problem existed for prior art, the present invention provides a kind of rainfall grain based on frequency angle two-dimensional mixing interpolation Son scattering computing accelerated method, saves the plenty of time when calculating frequency and angle scanning and effectively raises rainfall grain The arithmetic speed of subgroup bandwidth scattering.

The present invention is directed to what above-mentioned technical problem was mainly addressed by following technical proposals:

Rainfall KPT Scatter computing accelerated method based on frequency angle two-dimensional mixing interpolation, comprises the steps:

(1) single order of all sampled points in addition to two boundary points is tried to achieve at angle domain application cubic spline interpolation Derivative；

(2) at wideband band domain application Stoer-Bulirsch algorithm, rational interpolating function is obtained by calculating sampled point；

(3) use two-dimensional mixing adaptive interpolation method that angle domain and the frequency domain of target wave band are carried out two-dimensional adaptive Sampling, and sampled result is carried out structure division from coarse to fine.

Further, described step (1) is tried to achieve in addition to two boundary points at angle domain application cubic spline interpolation The first derivative of all sampled points specifically include:

Assume sampling node in angle domainWhole angular range is divided into multiple interval；? Each intervalInBe oneThe three rank multinomials that place launches, by Hermite interpolation theory,Insert Value formula is:

WhereinRepresent sampled point,In order to calculate angular rangeInValue, need ArriveWith WithFirst derivative can obtain under zero natural boundary conditions, then obtain it by formula (2) The first derivative of his sampled point；

Wherein h_iRepresentWithBetween interval,

Further, described step (2), at wideband band domain application Stoer-Bulirsch algorithm, obtains by calculating sampled point To rational interpolating function particularly as follows:

2a) useKnown sampled data is to indicate Manage function, and known sampled data is (f_i,P(f_i)) i=1 ..., k；Wherein, A (f), B (f) are that number of times is less than respectively The multinomial of p and q, P (f_i) it is sampled point f_iCorresponding sampled value；

2b) make R_i,1=P (f_i) i=1 ..., k, R (f)=R_1,k, apply Stoer-Bulirsch algorithm, according to formula (3)- (5) recursive resolve takes the post as Jie's rational function f interpolation R (f) at an arbitrary position:

R_j,k=((f-f_j)R_j+1,k-1+(f_j+k-f)R_j,k-1)/(f_j+k-f_j) (3)

R_j,k=(f_j+k-f_j)/((f-f_j)/R_j+1,k-1+(f_j+k-f)/R_j,k-1) (4)

R_j,k=R_j+1,k-2+(f_j+k-f_j)/((f-f_j)/(R_j+1,k-1-R_j+1,k-2)+(f_j+k-f)/(R_j,k-1-R_j+1,k-2)) (5)。

5, further, described step (3) use two-dimensional mixing adaptive interpolation method to the angle domain of target wave band and Frequency domain carries out two-dimensional adaptive sampling, and sampled result is carried out structure from coarse to fine divide particularly as follows:

3a) take four sampled point (f₁,φ₁),(f₁,φ_n),(f_m,φ₁) and (f_m,φ_n) do linear interpolation, then obtain I (f₁,φ₁),I(f₁,φ_n),I(f_m,φ₁) and I (f_m,φ_n)；

3b) assume in one plane have m × n sampled point, behavior m, be classified as n, obtain approximating thunder by linear interpolation Reach cross section rcs_a1 (f, φ), obtain radar cross section rcs_a2 (f, φ) by cubic spline interpolation approximation；

3c) in each intervalChoose the point that in these two curve of error, mean error is maximum Position i.e. φ_tj, obtain line of sampling；

Maximum error ME in 3d): as convergence error CE (φ), more than step 3c), now sampled point is best, otherwise will Interpolating functionCurve add in sampling interval；

3e) at each interval [f_si,f_si+1] (i=1,2 ..., m), pass through Stoer-Bulirsch algorithm when ME occurs Obtain line of sampling；

3f) when CE (f) occurring more than step 3e) in ME time, retain now corresponding sampled point, otherwise by interpolating functionCurve add in sampling interval；

3g) if each intervalMeet step 3d) in condition, and at [f_si,f_si+1] interval satisfied step Rapid 3f) the middle condition occurred, then sampling terminates.

The invention has the beneficial effects as follows:

1. application Stoer-Bulirsch algorithm calculating sampled point obtains rational interpolating function, and need not solution matrix Inverse, so it can obtain rational interpolating function by many sampled points, avoid singularity problem simultaneously, significantly save meter The operation time of calculation machine.

2. cubic spline interpolation obtains the first derivative rather than repeatedly of each sampling node by another way Solve system of linear equations, greatly reduce scanning cost when angle domain scans.

Accompanying drawing explanation

In order to be illustrated more clearly that the embodiment of the present invention or technical scheme of the prior art, below will be to embodiment or existing In having technology to describe, the required accompanying drawing used is briefly described, it should be apparent that, the accompanying drawing in describing below is only this Some embodiments of invention, for those of ordinary skill in the art, on the premise of not paying creative work, it is also possible to Other accompanying drawing is obtained according to these accompanying drawings.

Fig. 1 is the structural representation of rainfall KPT Scatter computing accelerated method based on frequency angle two-dimensional mixing interpolation；

Fig. 2 is the shape approximate diagram of raindrop particle；

Fig. 3 is the cross section of different size rainfall particle；

Fig. 4 (a) is the schematic diagram of model I；

Fig. 4 (b) is model I RCS curve under vertical polarization；

Fig. 5 (a) is the schematic diagram of model II；

Fig. 5 (b) is model II RCS curve under vertical polarization；

Fig. 6 (a) is the schematic diagram of model III；

Fig. 6 (b) is model III RCS curve under vertical polarization.

Detailed description of the invention

Below in conjunction with the accompanying drawings the preferred embodiments of the present invention are described in detail, so that advantages and features of the invention energy It is easier to be readily appreciated by one skilled in the art, thus protection scope of the present invention is made apparent clear and definite defining.

The main thought of the present invention is: the faradic current angle domain at wave band interested (i.e. target wave band) applies three samples Bar algorithm；Stoer-Bulirsch algorithm is applied in broadband rainfall population Analysis of Electromagnetic Character；By two-dimensional mixing certainly Adapt to interpolation technique, use by roughly to fine hierarchy as a process repeatedly, uneven to generate one group Sampling node model, improves operation efficiency.

Rainfall KPT Scatter computing accelerated method based on frequency angle two-dimensional mixing interpolation, including:

1, at angle domain application cubic spline interpolation, it is not necessary to repeatedly solve system of linear equations to obtain each sampled point Derivative, can directly try to achieve the first derivative of all sampled points in addition to two boundary points；Particularly as follows: assume in angle domain Middle sampling nodeWhole angular range is divided into multiple interval；In each intervalIn Be oneThe three rank multinomials that place launches, by Hermite interpolation theory,Formula for interpolation be:

WhereinRepresent sampled point,In order to calculate angular rangeInValue, need ArriveWith WithFirst derivative can obtain under zero natural boundary conditions, then obtain it by formula (2) The first derivative of his sampled point；

Wherein h_iRepresentWithBetween interval,

2, at wideband band domain application Stoer-Bulirsch algorithm, rational interpolating function is obtained by calculating sampled point；Tool Body is:

First, useKnown sampled data is for representing Rational function, and known sampled data is (f_i,P(f_i)) i=1 ..., k；Wherein, A (f), B (f) are that number of times does not surpasses respectively Cross the multinomial of p and q, P (f_i) it is sampled point f_iCorresponding sampled value；

Then, R is made_i,1=P (f_i) i=1 ..., k, R (f)=R_1,k, apply Stoer-Bulirsch algorithm, according to formula (3)-(5) recursive resolve takes the post as Jie's rational function f interpolation R (f) at an arbitrary position；The method is made without matrix inversion fortune Calculate, thus fundamentally avoiding singularity problem.

R_j,k=((f-f_j)R_j+1,k-1+(f_j+k-f)R_j,k-1)/(f_j+k-f_j) (3)

R_j,k=(f_j+k-f_j)/((f-f_j)/R_j+1,k-1+(f_j+k-f)/R_j,k-1) (4)

R_j,k=R_j+1,k-2+(f_j+k-f_j)/((f-f_j)/(R_j+1,k-1-R_j+1,k-2)+(f_j+k-f)/(R_j,k-1-R_j+1,k-2)) (5)。

Wherein, the concrete calculation process of application Stoer-Bulirsch algorithm is as follows:

1st step: calculate the actual value EM (f of broadband two ends f1, f2₁)、EM(f₂), and by f₁、EM(f₁)、f₂、EM(f₂) Join in set of samples

2nd step: utilize formula (3) and formula (4) as recurrence formula, calculate each frequency-of-interest point f in broadband respectively_t's Value P₁(f_t) and P₂(f_t), and find out maximum difference point f₃And difference value error12.

$error12=|P_1\left(f_3\right)-P_2\left(f_3\right)|=\underset{i}{max}|P_1\left(f_{ti}\right)-P_2\left(f_{ti}\right)|---\left(6\right)$

3rd step: if maximum error12 is less than error threshold eps set, then terminate to calculate, forward the 6th step to.As Really error12 is bigger than eps, then calculate f₃Actual value EM (the f at place₃), and by f₃With EM (f₃) join in set of samples.

4th step: utilize formula (3) and formula (5) as recurrence formula, calculate each frequency-of-interest point f in broadband respectively_t's Value P₁(f_t) and P₃(f_t), and find out maximum difference point f₄And difference value error13.

$error13=|P_1\left(f_4\right)-P_3\left(f_4\right)|=\underset{i}{max}|P_1\left(f_{ti}\right)-P_3\left(f_{ti}\right)|---\left(7\right)$

5th step: if maximum error13 is less than error threshold eps set, then terminate to calculate, forward the 6th step to.As Really error13 is bigger than eps, then calculate f₄Actual value EM (the f at place₄), and by f₄With EM (f₄) join in set of samples, then turn Enter the 2nd step and carry out the calculating of next sampled point.

6th step: according to existing data formula (3) modeling in sampling, clicks on row interpolation and calculates each frequency-of-interest.

3, use two-dimensional mixing adaptive interpolation method that angle domain and the frequency domain of target wave band are carried out two-dimensional adaptive Sampling, and sampled result is carried out structure division from coarse to fine；Particularly as follows:

1st step: take four sampled point (f₁,φ₁),(f₁,φ_n),(f_m,φ₁) and (f_m,φ_n) do linear interpolation.Then To I (f₁,φ₁),I(f₁,φ_n),I(f_m,φ₁) and I (f_m,φ_n)。

2nd step: assume in one plane there be m × n sampled point, behavior m, it is classified as n, is approximated by linear interpolation Radar cross section rcs_a1 (f, φ), obtains radar cross section rcs_a2 (f, φ) by cubic spline interpolation approximation.

3rd step: in each intervalChoose mean error in these two curve of error maximum The position of point i.e. φ_tj, obtain line of sampling, such as

4th step: as convergence error CE (φ), more than maximum error ME in step 3, now sampled point is best, otherwise will Interpolating functionCurve add in sampling interval.

5th step: at each interval [f_si,f_si+1] (i=1,2 ..., m), calculated by Stoer-Bulirsch when ME occurs Method obtains line of sampling.

6th step: when occurring that CE (f) is more than the ME in step 5, now sampled point is best, otherwise by interpolating functionCurve add in sampling interval.

7th step: if each intervalMeet the condition in step 4, and at [f_si,f_si+1] interval satisfied step The condition occurred in rapid 6, the most adaptively sampled process just finishes.

Shown below is several the numerical example this method is illustrated.

In order to verify that S-B (Stoer-Bulirsch) adaptively sampled method is calculating having in broadband rainfall population Effect property, this section gives several the numerical example and is illustrated.The iterative algorithm of solution matrix equation uses GPBi-CG, repeatedly The initial value in generation is all set to zero, and greatest iteration step number is set to 10000, and iteration precision is 10^-5.In this section, all examples calculate Working environment is Intel Core (TM) II Quad processor, 3.50GB internal memory, 2.83GHz dominant frequency, single precision.

The random size in the most individual position of model arbitrary spheroidal particle group

At 0.3*0.3*0.3m³In space, random distribution the spherical rainfall particle of 1000 radius 0-2mm, unknown quantity Number is 35308.The dielectric constant of raindrop particle changes along with frequency, between 4+0.04i and 3.5245+0.08755i.

The random size in II.500 position of model arbitrary elliposoidal population

In actual rainfall, due to the double influence by gravity and air drag, most rainfall shape of particle It is not that standard is spherical, but is rendered as the structure of similar ellipsoid, as shown in Figure 2.Wherein, the semi-minor axis semiaxis of ellipsoid is a length of a₃, major semiaxis half axial length is respectively a₁,a₂(a₁=a₂).Then, the shape of rainfall particle just can be by the axle ratio (a of ellipsoid₁/ a₃) describe.

This example devises at 0.3*0.3*0.3m³In space, random distribution the elliposoidal fall of 500 semi-minor axis 0-2mm Rain particle, ellipsoid axle is than change at random between 1 and 2, and unknown quantity number is 34604.The dielectric constant of raindrop particle is along with frequency Change, between 4+0.04i and 3.4+0.1i.In the frequency range of 1-13GHz, every 0.01GHz takes a sampled point, uses SM/CG method point by point scanning need to calculate 1201 Frequency points.

The random size in III.500 position of model meets the rainfall population of Keenan distribution

It is situation during spheroid that a upper example describes rainfall particle, and in true rainfall, the shape of rainfall particle Shape and size have close relationship, Fig. 3 to give the cross sectional representation of different size of raindrop and corresponding equivalent sphere.By Figure is it can be seen that along with the increase of raindrop volume, real rain drop shapes presents the feature becoming closer to oblate ellipsoid.

In this example, it will be assumed that at 0.5*0.5*0.5m³In space, random distribution the ellipse of 500 semi-minor axis 0-2mm Spherical rainfall particle, their axle meets Keenan et al. model than relation, and unknown quantity number is 27498.Raindrop grain The dielectric constant of son changes along with frequency, between 4+0.04i and 3.3495+0.105i.Frequency range at 10-23GHz Interior every 0.01GHz takes a sampled point, need to calculate 1301 Frequency points with SM/CG method point by point scanning.

Fig. 4,5 and 6 sets forth above three example and directly use SM/CG method and self adaptation S-B interpolation method Result of calculation.It can be seen that the result of calculation in the case of three kinds is all coincide good.In example 1, as shown in Figure 4, at 1- With 20MHz as sampling interval in the frequency range of 20GHz, 951 Frequency points need to be calculated when using SM/CG method point by point scanning, And when combining adaptive S-B interpolation method, the sampled point number directly calculated reduces to 224.In example 2,3, such as Fig. 5,6 institutes Showing, frequency sampling is spaced apart 10MHz, is respectively required for calculating 1201 and 1301 Frequency points when using SM/CG method point by point scanning, and When combining adaptive S-B interpolation method, the sampled point number directly calculated reduces to 387 and 310 respectively.

Table 1 gives three above example and the most directly uses SM/CG method and self adaptation S-B interpolation method when calculating Comparison between.There it can be seen that after using self adaptation S-B interpolation method, in the case of three kinds, calculate total frequency band RCS Time needed for value only has original 23.66%, 26.62% and 28.20% respectively.It should be noted that this sampling algorithm has There is self adaptation feature, therefore can automatically select the sample frequency point position being more suitable for and sampled point according to practical situation Number, i.e. RCS curve is more complicated, and when fluctuating more violent, sampling number is more, and RCS curve is relatively simple, samples when changing shallower Count less.Above example all shows, self adaptation S-B interpolation method is a kind of to calculate having of rainfall particle broadband electromagnetic property Efficacious prescriptions method.

The above, one combines Stoer-Bulirsch model two-dimensional adaptive hybrid interpolation technology based on cubic spline Apply the quick calculating in precipitation particles broadband electromagnetic property.Comparing traditional AWE technology, this new method is without Matrix Calculating Inverse, it is to avoid the calculating of first derivative, result of calculation shows, self adaptation S-B interpolation method is a kind of calculating rainfall particle broadband The effective ways of electromagnetic property.

Table 1 directly uses the contrast that SM/CG method calculates and uses the self adaptation S-B method of sampling to calculate

The above, the only detailed description of the invention of the present invention, but protection scope of the present invention is not limited thereto, and any The change expected without creative work or replacement, all should contain within protection scope of the present invention.Therefore, the present invention Protection domain should be as the criterion with the protection domain that claims are limited.

Claims (4)

1. rainfall KPT Scatter computing accelerated method based on frequency angle two-dimensional mixing interpolation, it is characterised in that: including:

(1) try to achieve the single order of all sampled points in addition to two boundary points at angle domain application cubic spline interpolation to lead Number；

(2) at wideband band domain application Stoer-Bulirsch algorithm, rational interpolating function is obtained by calculating sampled point；

(3) use two-dimensional mixing adaptive interpolation method that angle domain and the frequency domain of target wave band are carried out two-dimensional adaptive and adopted Sample, and sampled result is carried out structure division from coarse to fine.

Rainfall KPT Scatter computing accelerated method based on frequency angle two-dimensional mixing interpolation the most according to claim 1, It is characterized in that: described step (1) tries to achieve owning in addition to two boundary points at angle domain application cubic spline interpolation The first derivative of sampled point specifically includes:

Assume sampling node in angle domainWhole angular range is divided into multiple interval；Each IntervalInBe oneThe three rank multinomials that place launches, by Hermite interpolation theory,Interpolation public Formula is:

WhereinRepresent sampled point,In order to calculate angular rangeInValue, need to obtainWith WithFirst derivative can obtain under zero natural boundary conditions, then obtain other by formula (2) The first derivative of sampled point；

Wherein h_iRepresentWithBetween interval,

Rainfall KPT Scatter computing accelerated method based on frequency angle two-dimensional mixing interpolation the most according to claim 1, It is characterized in that: described step (2), at wideband band domain application Stoer-Bulirsch algorithm, obtains reasonable by calculating sampled point Interpolating function particularly as follows:

2a) useKnown sampled data is for representing reasonable letter Count, and known sampled data is (f_i,P(f_i)) i=1 ..., k；Wherein, A (f), B (f) are that number of times is less than p and q respectively Multinomial, P (f_i) it is sampled point f_iCorresponding sampled value；

2b) make R_i,1=P (f_i) i=1 ..., k, R (f)=R_1,k, apply Stoer-Bulirsch algorithm, pass according to formula (3)-(5) Return to solve and take the post as rank rational function f interpolation R (f) at an arbitrary position:

R_j,k=((f-f_j)R_j+1,k-1+(f_j+k-f)R_j,k-1)/(f_j+k-f_j) (3)

R_j,k=(f_j+k-f_j)/((f-f_j)/R_j+1,k-1+(f_j+k-f)/R_j,k-1) (4)

R_j,k=R_j+1,k-2+(f_j+k-f_j)/((f-f_j)/(R_j+1,k-1-R_j+1,k-2)+(f_j+k-f)/(R_j,k-1-R_j+1,k-2)) (5)。

Rainfall KPT Scatter computing accelerated method based on frequency angle two-dimensional mixing interpolation the most according to claim 1, It is characterized in that: described step (3) uses two-dimensional mixing adaptive interpolation method to enter angle domain and the frequency domain of target wave band Row two-dimensional adaptive is sampled, and sampled result is carried out structure from coarse to fine divide particularly as follows:

3a) take four sampled point (f₁,φ₁),(f₁,φ_n),(f_m,φ₁) and (f_m,φ_n) do linear interpolation, then obtain I (f₁, φ₁),I(f₁,φ_n),I(f_m,φ₁) and I (f_m,φ_n)；

3b) assume in one plane have m × n sampled point, behavior m, be classified as n, obtain approximating radar by linear interpolation horizontal Cross section rcs_a1 (f, φ), obtains radar cross section rcs_a2 (f, φ) by cubic spline interpolation approximation；

3c) in each interval(j=1,2 ..., n) choose the position of the point that mean error is maximum in these two curve of error Put namely φ_tj, obtain line of sampling；

Maximum error ME in 3d): as convergence error CE (φ), more than step 3c), now sampled point is best, otherwise by interpolation FunctionCurve add in sampling interval；

3e) at each interval [f_si,f_si+1] (i=1,2 ..., m), obtained by Stoer-Bulirsch algorithm when ME occurs Sampling line；

3f) when CE (f) occurring more than step 3e) in ME time, retain now corresponding sampled point, otherwise by interpolating functionCurve add in sampling interval；

3g) if each intervalMeet step 3d) in condition, and at [f_si,f_si+1] interval meets step The condition occurred in 3f), then sampling terminates.