CN102638318B - Method applicable to correcting equiamplitude vibration error of aerofoil conformal array - Google Patents

Abstract

The invention discloses a method for correcting the constant-amplitude vibration error of a wing conformal array, which mainly solves the problems of time variation and 2π ambiguity in the process of correcting the constant-amplitude vibration error of the wing conformal array in the prior art. The implementation process is as follows: first, assuming that the phase center errors of adjacent snapshots are the same, the signal subspaces at each time are estimated respectively; then, according to the estimated signal subspaces, the cost function of array element position vector estimation is established; using Newton iteration The method calculates the cost function and obtains the array element position vector at each moment; uses the curve fitting method to fit the array element position vector, and the difference between the fitted value and the ideal value of the array element position vector is the phase center error; finally, Perform phase center error correction on the received data according to the obtained phase center error value. The invention overcomes the time-varying phenomenon and the 2π fuzzy phenomenon of the wing conformal array phase center error, and can be applied to the correction of the constant amplitude vibration error of the wing conformal array.

Description

Be applicable to the continuous vibration error calibration method of the conformal array of wing

Technical field

The invention belongs to signal process field, relate to the phase center error calibration method of even linear array, can be used for the correction of the phase center error that the conformal array antenna of wing causes due to continuous vibration in aircraft flight process.

Background technology

Aircraft is in the middle of flight course, and wing and air-flow interact and can produce chatter phenomenon, and in the time that speed reaches a certain specific threshold, the amplitude of the flutter of aerofoil being caused by disturbance just remains unchanged, and claims that this speed is flutter speed.Generally, think that now the zitterbewegung of wing is simple harmonic oscillation.Change because the continuous vibration of wing causes with it conformal array antenna phase center, thereby produced phase center error.

The phase center error that wing continuous vibration causes has following characteristics:

1) airfoil variable is larger.Bay is far away apart from wing root, and phase center error is larger, easily causes that 2 π are fuzzy;

2) data dependence variation.Due to flutter of aerofoil, phase center error can slowly change along with the time, thereby caused the correlation variation between sample data.

At present, conventional phase center error calibration method, as add auxiliary array element, add the accurate known information source of direction etc., mainly for the constant situation of error, do not consider that flutter of aerofoil causes the time dependent situation of phase center error, and because airfoil variable is larger, easily cause that 2 π are fuzzy, adopt conventional method phase calibration errors of centration cannot obtain correct result.

Summary of the invention

The object of the invention is to the deficiency for conventional phase center error calibration method, the bearing calibration of the conformal array continuous vibration of a kind of wing phase center error has been proposed, with eliminate 2 π that cause due to flutter of aerofoil fuzzy and time change effect, effectively improve the correction accuracy of phase center error.

The technical scheme that realizes the object of the invention is that the phase center error change feature comparatively slowly first causing for flutter of aerofoil, supposes that the phase center error of adjacent snap is identical, estimates respectively the signal subspace in each moment; Then,, according to the signal subspace of estimating to obtain, adopt the method for subspace fitting to carry out the rough estimate of error; Finally, adopt the method fit error curve of curve, carry out the accurate estimation of error.Specific implementation step is as follows:

(1) suppose that the phase center error of adjacent n snap is identical, n is odd number, utilizes following formula to obtain the k time snap to receive the covariance matrix of data

$R_X^k=\frac{1}{n}\underset{l=k-\frac{n-1}{2}}{\overset{k+\frac{n-1}{2}}{\mathrm{\&Sigma;}}}x\left(l\right)x{\left(l\right)}^{H}k=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,L,---<1>$

Wherein, L is fast umber of beats, and x (l) represents the l time sample data, () ^hrepresenting matrix conjugate transpose; According to

estimate the signal subspace in this moment;

(2), according to the signal subspace of estimating to obtain, adopt the method based on subspace fitting to set up the cost function that element position vector is estimated:

Wherein,

represent the estimated value of element position vector, U _srepresent signal subspace,

represent array manifold matrix,

⁺representing matrix generalized inverse; Arg () represents to get plural argument, represent to get the transition formula evaluation minimum that makes in bracket

value, || || ₂represent to get 2 norms;

(3) utilize Newton iteration method to calculate above-mentioned cost function, obtain the element position vector in each moment, realize the rough estimate to phase center error;

(4) adopt curve-fitting method to carry out matching to the element position vector obtaining in step (3), the match value of this element position vector and the difference of ideal value are phase center error, realize the accurate estimation of phase center error;

(5) according to the phase center error of estimating to obtain, carry out the correction of phase center error.

The present invention compared with prior art, has the following advantages:

1) the present invention adopts Newton iteration method direct estimation element position vector, solves the element position vector that obtains and the difference of desirable element position vector and is phase center error, and 2 π that avoided conventional method phase calibration errors of centration to produce are fuzzy;

2) the present invention adopts the method for curve to carry out matching to array element position vector, and the time realization that has reduced phase center error resembles the impact on error correction precision.

Accompanying drawing explanation

Fig. 1 is flow chart of the present invention;

Fig. 2 is the analogous diagram of element position correcting action of the present invention with signal to noise ratio situation of change;

Fig. 3 is the analogous diagram of element position correcting action of the present invention with snap situation of change.

Embodiment

With reference to Fig. 1, performing step of the present invention is as follows:

Step 1 is estimated the signal subspace of each moment data.

(1a) suppose that the phase center error of adjacent n snap is identical, n is odd number, utilizes following formula to obtain the k time snap to receive the covariance matrix of data

$R_X^k=\frac{1}{n}\underset{l=k-\frac{n-1}{2}}{\overset{k+\frac{n-1}{2}}{\mathrm{\&Sigma;}}}x\left(l\right)x{\left(l\right)}^{H}k=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,L,---<3>$

Wherein, L is fast umber of beats, and x (l) represents the l time sample data, () ^hrepresenting matrix conjugate transpose;

(1b) covariance matrix to the k time snap reception data

carry out feature decomposition, be decomposed into the multiply accumulating form of characteristic value and characteristic vector, and characteristic value is arranged by order from big to small:

$R_X^k=\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_i{v}_i{v}_i^H,---<4>$

${\mathrm{\&lambda;}}_1\&GreaterEqual;{\mathrm{\&lambda;}}_2\&GreaterEqual;\&CenterDot;\&CenterDot;\&CenterDot;\&GreaterEqual;{\mathrm{\&lambda;}}_P>>{\mathrm{\&lambda;}}_{P+1}={\mathrm{\&lambda;}}_{P+2}=\&CenterDot;\&CenterDot;\&CenterDot;={\mathrm{\&lambda;}}_M={\mathrm{\&sigma;}}_n^2$

Wherein, M is array number, λ _ifor covariance matrix

characteristic value, v _ibe and eigenvalue λ _icorresponding characteristic vector, i=1,2 ... M; P is information source number,

for white noise power;

(1c) defined feature value λ ₁, λ ₂..., λ _pthe subspace that characteristic of correspondence vector is opened is signal subspace U _s:

U _s＝[v ₁，v ₂，…，v _P]。<5>

Step 2 adopts subspace fitting method to set up the cost function that element position vector is estimated.

(2a) establish

for element position vector:

Wherein, x _n, y _nbe respectively n abscissa and the ordinate after array element deformation, n=1,2 ..., M, () ^trepresenting matrix transposition;

(2b), by the following formula of coordinate figure substitution of element position vector, obtain the spatial domain steering vector a (θ of i information source _i):

$a\left({\mathrm{\&theta;}}_i\right)={[e^{j\frac{2\mathrm{\&pi;}}{\mathrm{\&lambda;}}{\mathrm{\&tau;}}_1\left({\mathrm{\&theta;}}_i\right)},e^{j\frac{2\mathrm{\&pi;}}{\mathrm{\&lambda;}}{\mathrm{\&tau;}}_2\left({\mathrm{\&theta;}}_i\right)},\&CenterDot;\&CenterDot;\&CenterDot;,e^{j\frac{2\mathrm{\&pi;}}{\mathrm{\&lambda;}}{\mathrm{\&tau;}}_M\left({\mathrm{\&theta;}}_i\right)}]}^{T}i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,P,<7>$

Wherein, θ _irepresent the arrival bearing of i information source, i=1,2 ..., P, P is information source number τ _n(θ _i)=[sin θ _i, cos θ _i] [x _n, y _n] ^t, n=1,2 ..., M, e ⁽⁾represent the exponential function take natural logrithm e the end of as, () ^trepresenting matrix transposition, λ is carrier wavelength, j represents imaginary unit;

(2c) by spatial domain steering vector a (θ _i), obtain array manifold matrix

for:

Wherein, () ^trepresenting matrix transposition, θ=[θ ₁, θ ₂..., θ _p] ^trepresent arrival bearing's vector;

(2d) according to array manifold matrix

the cost function that obtains the estimation of element position vector is:

，<9>

Wherein,

represent the estimated value of element position vector, arg () represents to get plural argument,

represent to get the transition formula evaluation minimum that makes in bracket

value, || || ₂represent to get 2 norms, () ⁺represent to ask Generalized Inverse Matrix; I _mrepresent M rank unit matrix, matrix P _awith C be intermediate variable, w is weight matrix, and Tr () represents to ask matrix trace, for cost matrix:

Weight matrix W obtains by following formula:

$W={({\mathrm{\&Lambda;}}_s-{\widehat{\mathrm{\&sigma;}}}_n^2)}^2{\mathrm{\&Lambda;}}_s^{-1},---<11>$

Wherein, () ²represent squared operation, () ^-1the representing matrix operation of inverting, Λ _sfor intermediate variable, Λ _s=diag (λ ₁, λ ₂..., λ _p), diag () represents vector to form diagonal matrix,

represent noise power estimation value:

${\widehat{\mathrm{\&sigma;}}}_n^2=\frac{1}{M-P}\mathrm{Tr}\left(C{R}_X^k\right).---<12>$

Step 3 utilizes Newton iteration method to calculate element position vector estimate cost function, obtains the element position vector in each moment, realizes the rough estimate to phase center error.

(3a) hypothesis array element initial position vector for desirable element position vector;

(3b) establishing k is iterations, makes k=0, by array element initial position vector

substitution formula <3>, adopts Newton iterative to calculate the element position vector of the 1st iteration

wherein,

represent the element position vector that the k+1 time iteration obtains,

represent the element position vector that the k time iteration obtains, the step factor that β is iteration,

represent cost matrix

the k time iterative value of second dervative,

represent cost matrix

the k time iterative value of first derivative, () ^-1the representing matrix operation of inverting;

provided by following formula:

Wherein, Re () represents to get real part operation, () ^trepresenting matrix matrix transpose operation,

for intermediate variable, provided by following formula:

$\tilde{J}=\left[\begin{array}{cc}\tilde{I}& \\ & \tilde{I}\end{array}\right],---<15>$

for M × (M-1) rank matrix:

$\tilde{I}=\left[\begin{array}{c}{0}^T\\ I_{M-1}\end{array}\right],---<16>$

M is array number, ' 0 ' expression null matrix, I _m-1for M-1 rank unit matrix;

B is intermediate variable, obtains by following formula:

Wherein,

i _mrepresent M rank unit matrix, the operation of ' ⊙ ' representing matrix dot product, H _xX, H _xY, H _yX, H _yYfor intermediate variable, obtain by following formula:

，<18>

() ^hrepresenting matrix conjugate transpose, w is weight matrix;

with

for intermediate variable, provided by following formula respectively:

，<19>

Wherein, f ₀represent carrier frequency, c represents the light velocity, and j is imaginary unit, Λ _sinand Λ _cosjust be respectively, cosine diagonal matrix:

Λ _sin＝diag([sinθ ₁，sinθ ₂，…，sinθ _P] ^T)

Λ _cos＝diag([cosθ ₁，cosθ ₂，…，cosθ _P] ^T)

θ _ibe the arrival bearing of i information source, i=1,2 ..., P, P is number of source, diag () represents vector to form diagonal matrix;

obtain by following formula:

Wherein, vecd () represents that the diagonal element of getting matrix forms column vector;

(3c) make k=k+1, by the k time iterative value of element position vector

substitution formula <3>, the k+1 time iterative value of calculating element position vector

If (3d) λ ₀for carrier wavelength, iteration stopping,

be the element position vector that estimation obtains; Otherwise continue to carry out (3c).

Step 4 adopts curve-fitting method to carry out matching to array element position vector, and the match value of this element position vector and the difference of ideal value are phase center error, realizes the accurate estimation of phase center error.

(4a) establish

represent the element position vector being obtained by the n time snap data estimation:

Wherein,

represent respectively abscissa and the ordinate value of i element position, i=2 ..., M;

(4b) make matching vector d be:

d＝[1，2，…，L] ^T，<23>

If d _irepresent i the element of vector d, i=1,2 ..., L, with three rank fitting of a polynomial error curves, H represents observation matrix, H is:

$H=\left[\begin{array}{cccc}{d}_1^3& d_1^2& d_1& 1\\ d_2^3& d_2^2& d_2& 1\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ d_L^3& d_L^2& d_L& 1\end{array}\right];---<24>$

(4c) establish I _xm, I _ymrepresent abscissa and the ordinate value of m array element error vector,

abscissa and the ordinate value of m the array element error vector that expression matching obtains, m=1,2 ..., M,

Order

$I_{\mathrm{Xm}}={[x_m^1,x_m^2,\&CenterDot;\&CenterDot;\&CenterDot;,x_m^L]}^T$ ，<25>

$I_{\mathrm{Ym}}={[y_m^1,y_m^2,\&CenterDot;\&CenterDot;\&CenterDot;,y_m^L]}^{T}m=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,M$

Obtain according to curve-fitting method:

${\hat{I}}_{\mathrm{Xm}}={\left(H^{T}H\right)}^{-1}{H}^T{I}_{\mathrm{Xm}}$ <26>

${\hat{I}}_{\mathrm{Ym}}={\left(H^{T}H\right)}^{-1}{H}^T{I}_{\mathrm{Ym}}$

Respectively with

with

represent vector

with

n element, m=1,2 ..., M, n=1,2 ..., L; Order

represent to estimate the element position vector of the n time snap of m array element obtaining,

for:

${\hat{I}}_{\mathrm{mn}}=[{\hat{I}}_{\mathrm{Xm}}^n,{\hat{I}}_{\mathrm{Ym}}^n];---<27>$

(4d) establish Q _mnbe the element position vector ideal value of the n time snap of m array element, m=1,2 ..., M, n=1,2 ..., L, the n time snap of m array element phase center error estimate be:

${\mathrm{\&phi;}}_{\mathrm{mn}}=Q_{\mathrm{mn}}-{\hat{I}}_{\mathrm{mn}}.---<28>$

Step 5, according to the phase center error of estimating to obtain, is carried out the correction of phase center error.

(5a) by estimate the n time snap of m array element of obtaining phase center error estimate φ _mnsubstitution formula <19>, obtain the n time snap phase center error vector estimated value Φ _n:

${\mathrm{\&Phi;}}_{\text{n}}={[e^{j\frac{2\mathrm{\&pi;}}{\mathrm{\&lambda;}}{\mathrm{\&phi;}}_{1n}},e^{j\frac{2\mathrm{\&pi;}}{\mathrm{\&lambda;}}{\mathrm{\&phi;}}_{2n}},\&CenterDot;\&CenterDot;\&CenterDot;,e^{j\frac{2\mathrm{\&pi;}}{\mathrm{\&lambda;}}{\mathrm{\&phi;}}_{\mathrm{Mn}}}]}^T,---<29>$

Wherein, m=1,2 ..., M, n=1,2 ..., L, M is array number, L is fast umber of beats;

(5b) according to phase center error vector estimated value Φ _n, the n time fast beat of data x (n) carried out to phase center error correction:

$\hat{x}\left(n\right)={\left[\mathrm{diag}\left({\mathrm{\&Phi;}}_n\right)\right]}^{-1}x\left(n\right),---<30>$

Wherein,

represent the correction result to the n time fast beat of data x (n), diag () represents vector to be formed to diagonal matrix, () ^-1representing matrix is inverted.

Effect of the present invention can be by illustrating the processing of emulated data below:

1 experimental situation and condition

Experimental situation: experiment adopts the half-wavelength uniform line-array that an array number is 8, and linear array is installed along aerofoil surface.

Simulated conditions: wing length is 4.2m, tip largest deformation amount is 0.21m, flutter of aerofoil frequency is 11.14Hz, flutter angle is [0.05rad, 0.05rad], and first phase is 0, radar operation wavelength is 1.2m, pulse repetition frequency is 1400Hz, and far field exists two mutual incoherent echo signals, and its arrival bearing is respectively-20 °, 10 °.

2 experiment contents and result

If I _mnbe the element position vector actual value of the n time snap of m array element, m=1,2 ..., M, n=1,2 ..., L, definition element position correcting action is:

$\frac{1}{\mathrm{MN}}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{L}{\mathrm{\&Sigma;}}}{\left|\right|{\hat{I}}_{\mathrm{mn}}-I_{\mathrm{mn}}\left|\right|}_2.---<31>$

Experiment 1, makes the signal to noise ratio of echo signal change, observes under different signal to noise ratio conditions, and when accurately estimating and accurately estimating, the situation that element position estimated bias changes, simulation result is as shown in Figure 2.Wherein, fast umber of beats is 500.

Experiment 2, the fast umber of beats that order receives data changes, observes under different snap conditions, when accurately estimating and accurately estimating, the situation that element position estimated bias changes, simulation result is as shown in Figure 3.Wherein, signal-to-noise ratio settings is 20dB.

As can be seen from Figure 2, under the lower condition of signal to noise ratio, the effect that adopts the inventive method to carry out phase center estimation error is better than does not carry out the accurately effect of estimation, and along with the raising of signal to noise ratio, the element position estimated bias curve of two kinds of methods overlaps gradually.

As can be seen from Figure 3, under the signal to noise ratio 20dB condition of setting, adopt the inventive method to carry out phase center estimation error and be greater than at 350 o'clock at fast umber of beats in this experiment, element position estimated bias tends towards stability; And the method phase estimation deviation of not carrying out accurately estimation is always larger, cannot reach the object of estimating phase center error, can find out from testing 1 result, while accurately estimation, need to be in the time that signal to noise ratio be greater than 30dB, performance could approach by the inventive method, and this condition is very inaccessible in actual applications.

Claims (4)

1. the conformal array continuous vibration of a wing error calibration method, comprises the steps:

(1) suppose that the phase center error of adjacent n snap is identical, n is odd number, utilizes following formula to obtain the k time snap to receive the covariance matrix of data

Wherein, L is fast umber of beats, and x (l) represents the l time sample data, () ^hrepresenting matrix conjugate transpose;

(2) basis

estimated signal subspace:

(2a) covariance matrix to the k time snap reception data

carry out feature decomposition, be decomposed into the multiply accumulating form of characteristic value and characteristic vector, and characteristic value is arranged by order from big to small:

Wherein, M is array number, λ _ifor covariance matrix

characteristic value, v _ibe and eigenvalue λ _icorresponding characteristic vector, i=1,2 ... M; P is information source number,

for white noise power;

(2b) defined feature value λ ₁, λ ₂..., λ _pthe subspace that characteristic of correspondence vector is opened is signal subspace U _s:

U _S＝[v ₁,v ₂,…,v _P];

(3), according to the signal subspace of estimating to obtain, the method for employing based on subspace fitting set up the cost function of estimation error:

Wherein,

represent the estimated value of element position vector, U _srepresent signal subspace,

represent array manifold matrix, representing matrix

generalized inverse; Arg () represents to get plural argument, represent to get the transition formula evaluation minimum that makes in bracket value, || || ₂represent to get 2 norms;

(4) utilize Newton iteration method to calculate above-mentioned cost function, obtain the element position vector in each moment, realize the rough estimate to phase center error;

(5) adopt curve-fitting method to carry out matching to the element position vector obtaining in step (3), the match value of this element position vector and the difference of ideal value are phase center error, realize the accurate estimation of phase center error;

(6) according to the phase center error of estimating to obtain, carry out the correction of phase center error.

2. the conformal array continuous vibration of wing according to claim 1 error calibration method, wherein the described Newton iteration method of utilizing of step (4) is calculated above-mentioned cost function, carries out as follows:

(4a) hypothesis array element initial position vector

for desirable element position vector;

(4b) establishing k is iterations, makes k=0, by array element initial position vector

substitution formula <3>, adopts Newton iterative to calculate the element position vector of the 1st iteration

Wherein,

represent the element position vector that the k+1 time iteration obtains,

represent the element position vector that the k time iteration obtains, the step factor that β is iteration,

represent cost matrix

the k time iterative value of second dervative,

represent cost matrix

the k time iterative value of first derivative, () ^-1the representing matrix operation of inverting;

provided by following formula:

Wherein, Re () represents to get real part operation, () ^trepresenting matrix matrix transpose operation, for intermediate variable, provided by following formula:

for M × (M-1) rank matrix:

M is array number, ' 0' represents null matrix, I _m-1for M-1 rank unit matrix;

B is intermediate variable, obtains by following formula:

Wherein,

i _mrepresent M rank unit matrix, ' the operation of ⊙ ' representing matrix dot product, H _xX, H _xY, H _yX, H _yYfor intermediate variable, obtain by following formula:

() ^hrepresenting matrix conjugate transpose,

w is weight matrix;

with

for intermediate variable, provided by following formula respectively:

Wherein, f ₀represent carrier frequency, c represents the light velocity, and j is imaginary unit, Λ _sinand Λ _cosjust be respectively, cosine diagonal matrix:

Λ _sin＝diag([sinθ ₁,sinθ ₂,…,sinθ _P] ^T)<10>

Λ _cos＝diag([cosθ ₁,cosθ ₂,…,cosθ _P] ^T)，

θ _ibe the arrival bearing of i information source, i=1,2 ..., P, P is number of source, diag () represents vector to form diagonal matrix;

obtain by following formula:

Wherein, vecd () represents that the diagonal element of getting matrix forms column vector;

(4c) make k=k+1, by the k time iterative value of element position vector substitution formula <3>, the k+1 time iterative value of calculating element position vector

If (4d) λ is carrier wavelength, iteration stopping,

be the element position vector that estimation obtains; Otherwise continue to carry out (4c).

3. the conformal array continuous vibration of wing according to claim 1 error calibration method, wherein the described employing curve-fitting method of step (5) carries out matching to the element position vector obtaining in step (3), carries out as follows:

(5a) establish

represent the element position vector being obtained by the n time snap data estimation:

Wherein,

represent respectively abscissa and the ordinate value of i element position, i=2 ..., M;

(5b) make matching vector d be:

d＝[1,2, …,L] ^T， <13>

If d _irepresent i the element of vector d, i=1,2 ..., L, with three rank fitting of a polynomial error curves, H represents observation matrix, H is:

(5c) establish I _xm, I _ymrepresent abscissa and the ordinate value of m array element error vector,

abscissa and the ordinate value of m the array element error vector that expression matching obtains, m=1,2 ..., M,

Order

Obtain according to curve-fitting method:

Respectively with

with

represent vector

with

n element, m=1,2, ..., M, n=1,2 ..., L; Order

represent to estimate the element position vector of the n time snap of m array element obtaining,

for:

(5d) establish Q _mnbe the element position vector ideal value of the n time snap of m array element, m=1,2 ..., M, n=1,2 ..., L, the n time snap of m array element phase center error estimate be:

4. the conformal array continuous vibration of wing according to claim 1 error calibration method, wherein the described phase center error obtaining according to estimation of step (5), carries out the correction of phase center error, carries out as follows:

(5a) by estimate the n time snap of m array element of obtaining phase center error estimate φ _mnsubstitution formula <19>, obtain the n time snap phase center error vector estimated value Φ _n:

Wherein, m=1,2 ..., M, n=1,2 ..., L, M is array number, and L is fast umber of beats, and λ is carrier wavelength;

(5b) according to phase center error vector estimated value Φ _n, the n time fast beat of data x (n) carried out to phase center error correction:

Wherein,

represent the correction result to the n time fast beat of data x (n), diag () represents vector to be formed to diagonal matrix, () ^-1representing matrix is inverted.

Images