- Email: [email protected]
Frequency-Domain GTLS Identification Combined with Time-Frequency Filtering for Flight Flutter Modal Parameter Identification
V ol. 19
No. 1
CH INES E JO U RN A L O F AER ON A U TICS
February 2006
Frequency Domain GTLS Identification Combined with Time Frequency Filtering for Flight Flutter Modal Parameter Identification T ANG Wei, SHI Zhong ke, LI Hong chao ( Coll ege of A utomat ion , Northw ester n Polytechnical Uni versity , X i an Abstract:
710072 , China)
T he aim of this paper is to present a new method for flig ht flutter modal parameter identifi
cat ion in noisy environment. T his method employs a time frequency ( T F ) filter to reduce the noise be fore identificatio n, w hich depends on the lo calization pr operty of sweep ex citation in T F domain. T hen, a gener alized total least square ( G T LS) identification algo rithm based o n stochastic framew ork is applied to t he enhanced data. System identificatio n w ith no isy data is transformed into a generalized total least square problem, and the solution is car ried out by the g eneralized sing ular value decomposition ( GSVD) to avoid the intensive nonlinear optimization co mputation. A nearly max imum likelihood pro perty can be achiev ed by optimally! weig hted g eneralized tot al least square. F inally, t he efficiency of the method is illustrated by means of flight test data. Key words:
flutter ; identification; time frequency filtering; GT L S
基于时频 域滤波及频域广义整体最小二乘辨识的飞机颤振模态参数辨识. 唐炜, 史忠科, 李洪超. 中国航空学报( 英文版) , 2006, 19( 1) : 44- 51. 摘
要: 提出一套适用于 噪声环 境的 飞机颤 振模 态参 数辨识 方法。为 减小 噪声 对辨 识结 果的 影
响, 首先设计了一种针对扫频激励的 时频滤 波器, 利用扫 频信号 及其响 应在时 频域分布 较为集 中 的特点, 有效去除噪声, 提高了试验数据的信噪比。为进一步提高 辨识精度, 提出了 一种基于 随机 模型的频域广义最 小二乘辨识算法。将 噪声条件 下的系 统辨识 问题转 化为广 义整体最 小二乘 问 题, 并采用线性的广义奇异值分解 求解模 型系数, 避免 了非线 性优化 的复杂 计算。通过 优化加 权 项, 获得了接近极大似然估计的辨识效果。最后, 通过试飞试验数据验证了方法的有效性。 关键词: 颤振; 辨识; 时频域滤波; 广义整体最小二乘 文章编号: 1000 9361( 2006) 01 0044 08
中图分类号: T N911. 7; N945. 14
文献标识码: A
New or modif ied aircraft of ten necessitates flig ht flut ter test s to verify the saf et y marg ins and
curat e results, w hen applied for flut ter modal anal ysis[ 2] . Recent ly , t he max imum likelihood estima
to prevent the cat ast rophic f lutt er. F light flut ter
t or ( ML E) has been developed and implemented for flutt er modal parameter ident ificat ion success
test ing cont inues to be a challenging research area because of t he concerns w ith cost, t ime and safet y in expanding the envelope of new or modif ied air craft . Hence, the procedures of flight t est s are de sired t o dramatically reduce t he flig ht flut ter t est time, increase the accuracy and reliability of t he [ 1]
modal param eter estimat ion . Aircraft in f light test ing dat a are typically
[ 3]
f ully . However, minimizat ion of m ax imum like lihood nonlinear cost funct ion requires both a pow erful iterat ive solver and a set of st art ing values that are ∀ close enough# t o global minimum. Conse quent ly, t here is no guarant ee at all that the opt i m ization algorit hm converges t o global minimum of ML cost funct ion.
charact erized by high noise level due t o at mospheric turbulence. For t his reason, the t radit ional param
T his paper focuses on invest ig at ing a new met hod f or f light flut ter modal paramet er ident if i
eter identification methods generally produce inac
cat ion in noisy environment. A generalized tot al
R eceived dat e: 2005 06 27; R evision received dat e: 2005 11 23 Foundation it em: N at ional N at ural S cience Foundat ion of China ( 60134010)
Frequency Domain GT LS Identif icat ion Com bined w ith T ime Frequency February 2006
∋ 45 ∋
Filt ering f or Flight Flut ter M odal Paramet er Identif icat ion 2
modal identif icat ion. Comparing w ith t he w avelet
1 - t2 - j 0| t| e e ( 2) 2 T he Fourier spect rum of M orlet w avelet is a
denoisng procedure for impulse ex citat ion discussed in Ref. [ 4] , a more pow erf ul noise reduct ion tech
shift ed Gaussian funct ion, and can be seen as a narrow band filt er w it h central frequency f 0
nique using t ime f requency f ilt er for linear sw eep
!(f ) = 2 e- 2 (f - f 0 ) ( 3) In order t o avoid phase dist ortion in signal re
least square ( GT L S ) estimator combined w ith time f requency f ilt ering is presented for accurate
excit at ion w ill be described in Sect ion 1. T hen, w it h the enhanced dat a, t he GT LS est imat or based on st ochastic f ramework will provide accurate iden t ification. Section 2 gives the detailed descript ion of t he GT LS algorit hm. Finally, t he met hod w ill be validat ed by real flig ht t est dat a.
1
T ime frequency Filtering for Noise Re duct ion Sw eep ( chirp) is t he most popular form of
( t) =
2
2
construction, only real M orlet w avelet is used 2 1 - t2 ( t) = e cos2 f 0 t ( 4) 2 T he localizat ion of M orlet in t ime and frequen cy domain allows for mapping t he signal into T F domain. By adjust ing t he paramet er f 0 , it is possi ble t o generate a family of narrowband filters w ith arbit rary cent ral f requency f i ( i = 1, 2, &, n ) . T he cont inuous w avelet transform maps t he one di
cont inuous excit at ion used in vibrat ion test , and is part icular com mon in flight flut ter test ing. T he
mensional sweep signal int o tw o dimension picture index ed by time ( b ) and frequency ( f ) . A t ime
part icular concent rat ion propert y of sweep sig nal in
frequency representat ion of signal is oft en called a scalogram, w hich is act ually t he pow er spect ral
time frequency ( T F ) domain allows for designing a
densit y | W ( b , f ) | of the signal over t he ( b , f )
filt er in T F f or response. 1. 1 Time frequency analysis and CWT T here have been many T F analysis tech
plane. 1. 2 TF filtering for sweep response
niques, such as spectrograms ( short time F ourier
T he t ypical form of linear sw eep signal is
transforms) , scalograms ( w avelets) and bilinear T F distribut ions in the Cohen! s class. Som e of
( 5) f 0 t + 1 r t2 2 where f 0 is the init ial f requency, r is t he frequency
them, such as spectrograms have a f ixed T F reso lut ion, the bilinear T F dist ributions have a high resolut ion but have crosst erms for mult icomponent signal. T he linear t echniques ( f or ex am ple w avelet transform) don! t have crossterms, and also have a g ood resolut ion by adjust ing dilat ion.
Hence,
w avelet transform is select ed as t he T F analysis tool in t his sect ion.
x ( t ) = A ( t ) cos
2
change rate. F ig. 1 show s an ex ample of scalogram of a 0 10 H z sw eep ex citation signal and m easured ac celeromet er response in flight . Obviously, t he sw eep show s a w ell concentrated propert y in T F, it localizes on a small reg ion. Applying sweep as an input to linear t im e invariant ( LT I) syst em, t he
F or a signal x ( t ) , the cont inuous w avelet transform can be defined as W( a, b) =
1 a
+ %
∃ x(t) - %
*
t - b dt ( 1) a
w here b is a t ranslation indicat ing t he localit y, a is a dilat ion or scale paramet er. ( t ) is a basic w avelet and
*
( t ) is the complex conjugat e of
( t) . One of t he most widely used funct ion in T F analysis is the Morlet w avelet[ 5] ,
( a)
Scalogram of sweep excitation sig nal
∋ 46 ∋
T A NG Wei, SHI Zho ng ke, L I Hong chao
CJA
D( b i , f j ) =
1
if | We ( b i , f j ) | > t 0
0
otherw ise
i = 1, &, N ; j = 1, &, M where D is the mask mat rix f rom w avelet coeff i cient We( b i , f j ) of sw eep ex citation by t hreshold ing , t 0 is t hreshold value. T hen, a mask operat ion is performed on wavelet coeff icients of noisy response, W ^ r ( b i , f j ) = D( b i , f j ) ∋ Wr ( b i , f j ) ( b)
Scalogram of accelerometer response
Fig . 1 Scalograms of flight flutter testing data
where W r denot es t he coefficients mat rix of noisy response and W ^ r represent s the coeff icients m at rix
response usually holds t he propert y and concen
aft er mask operat ion.
t rat es along sw eep line in T F. Hence, t he darker
1. 4
sloping area in ( b) mainly represent s t he true re sponse; the others in T F represent noise. If it is
T F domain t o t ime domain. H owever, t he f ilter
possible to locate t he region where t he chirp signal concent rat es in T F plane, a mask operat ion can be preformed to ext ract the t rue response. T he t ime frequency analysis sugg ests the f ollowing basic fil t ering procedure for response of sw eep ex citat ion, as shown in Fig. 2. ( 1) Perform t ime f requency analysis for sw eep and it! s response.
Response reconstruction Inverse t ransf orms w ill map t he sig nal from
bank reconstruction for discret e wavelet t ransform ( DWT ) can not be applied to CWT . A novel re construction in f requency domain is employed in this sect ion. Before int roducing t he met hod, express t he reconst ruct ion problem into a L S form, the f ilt ered response ^s is t hat minimizes
( 2) Remove noise component s in scalogram
(H^s - D ∋ * Hf (F = min (Hs - D ∋ * Hf (F = s
using mask operation. ( 3 ) Reconstruct t he response using clean
where operator H denot es t he w avelet t ransf orms.
scalogram.
f is the noisy response. Arithmet ic operator ( *∋ )
min ( Ws - D∋ * Wf ( F s
denot es the matrices element by element mult iplica tion. If H is a matrix, the normal solut ion of LS is s = H+ D ∋ * Hf
( 6)
where ∀ + # denot es pseudo inverse. How ever, t he comput at ion of CWT w ith convolut ion can not be represent ed by mat rix directly, + %
∃ x(t) !
1 a
Wx ( b, f 0 ) =
- %
* f0
b- t dt = a
x ( t ) * ! *f 0 , a ( t )
Fig . 2 T F filtering pr ocedure in T F plane
( 7)
Instead of Eq. ( 7) , CWT can be ex pressed in 1. 3
Mask operation Since t he input sw eep signal is known before,
frequency domain, Wx ( , f 0 ) = X( ) ∋ ! *f 0 , a ( )
( 8)
it s pat tern in T F domain may help in designing mask in T F domain for f ilt ering noise. T he com
where Wx (
mon met hod for mask desig n is t hresholding t he
Fourier t ransforms of Wx ( b, f 0 ) , x ( t ) and ! *f 0 , a
modulus of coef ficients of clean sw eep input in T F domain . T he procedure is g iven as
( t ) , respect ively.
, f 0 ) , X ( ) and !
* f , a( 0
) are t he
Write in matrix form in f requency domain,
Frequency Domain GT LS Identif icat ion Com bined w ith T ime Frequency February 2006
∋ 47 ∋
Filt ering f or Flight Flut ter M odal Paramet er Identif icat ion
F Wx = ∀ ∋ X
shown in F ig. 3. T he m easured input Um and out put Ym are dist urbed w it h noise. U0 and Y0 re spect ively are the t rue sw eep and its accelerometer response f ree noise, respectively. M u and M y are
w here Wx (
1,
f 1)
&
f N)
# &
#
F Wx = Wx (
1,
X = diag[ X( *
! f 1, a ( ∀ =
1) ,
# !
* f
N, a
(
Wx (
&
M,
M)]
*
! f 1, a (
& 1)
M,
f 1)
#
&, X(
&
1)
Wx (
the measurem ent noise in input and output, respec t ively. N g denot es t he at mospheric turbulence.
fN)
;
M
)
# ! *f
N, a
(
M)
Hence, LS can be w ritt en as in term of t he F ourier t ransf orm msin (F Ws - F ^ Wf (F = min ( ∀ ∋S- F ^ Wf (F s w here S = diag [ S (
1) ,
&, S (
M)
], S (
i)
Fig. 3
is
the discrete Fourier transforms ( DF T ) of estimat ed signal s . F ^ Wf denot es t he Fourier spect ra matrix of coeff icient s aft er m ask operation. T he minimum problem can be solved by t he
T he measurement s Um and Y m are related t he exact values U 0 and Y 0 by Um ( j ) = U0 ( j ) + N U ( j ) Ym ( j ) = Y 0 ( j ) + N Y ( j ) where N U ( j ) = MU ( j )
follow ing equat ion ∀ ∋s = F ^ Wf
( 9)
T he est imat ed solution is ^ = ∀+ F S ^ Wf T hen t he diagonal elem ent s of S ^ are inversely
Stochastic model of flight flutter test
N Y ( j ) = M Y ( j ) + H 0 ( j ) N g( j ) T
Def ine ∃Z = [ N U ( j k ) N Y ( j k ) ] , t he above equation can be w ritt en as Z m ( j k ) = Z 0 ( j k ) + ∃Z( j k ) ( 10) where T
transformed by IDF T ( inverse discret e Fourier tr
Z 0 = [ U0( j
k)
Y0( j
ansf orm ) to get t he t ime domain reconst ruct ion
Z m = [ Um( j
k)
Ym ( j k ) ] T
^s ( t) .
k) ]
It is reasonable t o make the following assump t ion : Assumpt ion 1. T he noise ∃Z ( j k ) is zero [ 6]
2
Generalized T otal Least Squares Estima tor for Ident ificat ion
T he common approach for modal ident if icat ion st art s from frequency response funct ions ( FRF ) ,
mean, complex normally distribut ed wit h covari ance m at rix E( ∃Z ( j k ) ∃Z H ( j k ) ) = Cz ( j k ) = %2U ( j
w hich are derived using a non parametric FRF est i
2 %Y U ( j
mator such as H 1 or H v . How ever, t he perf or mance of FRF estimator is sensitive to noise. Even w it h t he enhanced dat a, t he t radit ional frequency ident ification algorit hms perform poorly. For t his reason, the GT LS identificat ion algorit hm based on st ochast ic framew ork is applied, w hich improve t he accuracy of the paramet er est imation. 2. 1 Descri ption of flight flutter test with stochastic framework T he stochast ic model of flight f lutt er t est is
k)
E( ∃Z ( j
k)
k) T
%2U Y ( j 2 %Y (
j
k) k)
∃Z ( j l ) ) = 0
k, l
Assumpt ion 2. T he noise at each dif ference frequency is independent ly distributed, E ( ∃Z ( j 2. 2
k)
H
∃Z ( j l ) ) = 0,
k ) l
Parametric model A rational transfer funct ion model of order n/
d is used for frequency domain estimat or, N ( j , &) H 0 ( j , p) = D ( j , ∋) =
∋ 48 ∋
T A NG Wei, SHI Zho ng ke, L I Hong chao
&n ( j ) n + &+ &1 ( j ) + &0 d ∋d ( j ) + &+ ∋1 ( j ) + ∋0 w it h t he real coef ficients p = [ ∋0 & ∋d &
CJA
where A may cont ain errors in all elements.
( 11)
T he generalized tot al least squares ( GT L S)
&0
solut ion to the est imat ion in Eq. ( 14) is given
&n ] . Define p r ( r = 1, 2, &, d ) to be the poles of
by
[ 7]
arg min (( A - A ^ ) C- 1 (2F
the transfer funct ion, t he corresponding modal fre
( 15)
^ ,p A
quency and damping rat io can be obt ained as Im( p r ) Re( p r ) fr = , (r = 2 | pr | Obviously , t he modal parameters are derived
and subjected to A ^ p = 0, p T p = 1, w here C is a square root of t he colum n covariance matrix of A: H
H
C C= E{ ∃A ∃A} . A ^ is t he estimate of A. Eliminat ion of A ^ in Eq. ( 15) g ives the e
from the transfer function model, hence it is im
qualient cost funct ion minimized by GT L S estima
port ant to focus on estimat ing coef ficients of t rans
t or,
f er funct ion. 2. 3 GTLS in frequency domain
p T A H Ap T H p C Cp
arg pmin
If there are no model errors, t he out put and
subject to p T p = 1 ( 16)
input obey t he following system relat ion Y0 ( j ) = H 0 ( j ) U0 ( j ) U sing Eq. ( 11) , t he model equation is readily
p T A H Ap = ), and Eq. ( 16) p T CH Cp can be t ransformed into t he follow ing equation, Define
obtained,
min
Y0 ( j ) D( j , ∋) - U 0 ( j ) N ( j , &) = 0
AT Ap = )CT Cp
T he above equat ion can be w ritt en in m at rix
T he solut ion p t o the equat ion can be obtained by generalized sing ular value decom posit ion ( GSVD) of matrix pair ( A and C) [ 8] . T he vect or
form A0 p = 0
( 12)
in right m at rix corresponding t o minimum general
w here A0 = [
p Tp = 1
subject to
a H01
&
aH02
aH0F ] H ,
a 0k =
Z T0 ( j k )
ized singular value is the solut ion. Comparing w ith ML , this w ill avoid t he intensive nonlinear iterat ive
S( j k ) ,
S( j k ) = block diag(- [ 1, &, ( j ) ] , [ 1, &, ( j ) ] ) n
d
optimizat ion and relies on t he init ial point .
Put t ing the noisy value Eq. ( 10) into Eq. ( 12)
In order t o simplif y Eq. ( 15) , define
def ine the noisy matrix A, A = A 0 + ∃A
∗( j
( 13)
= a kp = Y( j
k)
U( j
w here ∃A = [ ∃aH1
∃aH2
2 WML ( j
& ∃aHF ] H,
k)
N (j
p ) = E { ∃∗( j
k,
E { ∃akpp
T
∃∋k = ∃Z ( j k ) S( j k ) . %2Y
U sing Eq. ( 13) , t he paramet er estimation can be formulat ed as a tot al least squares, looking for a
p)
+ k,
D( j
k,
p)
k)
∃a Hk }
%2U
k,
∃∗( j
p) k)
H
}=
=
N (j
p) D ( j
k, k,
p)
2
-
p) H )
T he cost f unction for minimizat ion in Eq. ( 15) can be rewritt en as
( 14)
F
∗
K GT LS ( p, Z m ) =
k,
2Re( %2U YN ( j
solut ion of Ap = 0
D( j
2
H
k)
∗( j
k)
2
k= 1 F
∗
W ML ( j
k)
2
=
k= 1 F
∗
Y( j
k)
D( j
k,
p ) - U( j
2
+ %U N ( j
2
k,
p)
k) N
(j
k,
p)
2
k= 1 F
∗
k= 1
2
{ %Y D ( j
k,
p)
2
( 17) 2
H
- 2Re( % UYN ( j k , p ) D( j k , p ) ) }
Frequency Domain GT LS Identif icat ion Com bined w ith T ime Frequency February 2006
2. 4
∋ 49 ∋
Filt ering f or Flight Flut ter M odal Paramet er Identif icat ion
Sample noise variance estimation
be used. T he parameter est imated f rom previous
Before using Eq. ( 17) t o est imat e the coeff i cients, it is needed t o est imat e the noise covariance and t he covariance f rom sample data. T he sample noise variance is usually derived from a set of data according t o st at ist ical definit ion, w hich is hardly
step is used to updat e the nex t st ep, and t he cost funct ion is w rit ten as F
K WGTLS ( pi , Z m ) =
∗
W-ML2 ( j k , p i- 1 ) | ∗( j k , pi ) | 2
∗
WML( j k , p i- 1 ) WML( j k, pi )
k= 1 F
-2
2
k= 1
to implement for limit ed t im e in f light flut ter test. F or this reason, a method to est imate the variance
( 19)
from one set of measurements is present ed. T he in
where K W GTLS ( pi , Z m ) denotes the cost funct ion
put sweep signal for excit at ion is know n before,
at i t h iterat ion step. p i- 1 is t he est imat ion derived
w hich is free noise. T o simplif y t he problem, only the error on the measured out put sequences will be
from previous step.
considered w ith %2U = 0, %2UY = 0.
reduced to M L cost funct ion
T he noise variance on t he out put is %2Y ( ) = E { N Y ( ) N HY ( ) } = H
K WGT LS, i = K ML / F
H
where K M L is the cost funct ion of M L est imat or.
E { N Y ( ) [ Ym ( ) - Y0 ( ) ] } =
T he f act t hat WGT L S and ML cost funct ion are e
H
E { N Y ( ) Ym ( ) } = H
quivalent indicat es t hat t he optimally! w eighted
H
E { Ym Y m - Y 0 Y m} = H
WGT L S has nearly M L property. Hence, t he
H
E ( Y m Ym ) - H 0 ( j ) E ( U0 Ym ) T he dat a are divided int o M segm ent s f or H 1 est imator, and H 0 is replaced by est imat ion H ^ 0, M
^%2Y (
∗
)=
Ym ( j )
Considering t he convergences of est imator, as the it eration i + % : p i - 1 = pi . Eq. ( 19) is readily
WGT L S can be used to generate the st art ing values that are ∀ close enough# t o t he act ual minimum for ML .
)- H ^ 0( j ) ∋
YHm ( j
3
i= 1 M
∗
H
U0 ( j ) Ym ( j ) / M = S Y Y( j ) -
In this sect ion, t he GT LS est imat or combin ing w ith time frequency filtering will be applied t o
i= 1
H ^ 0 ( j ) SUY ( j ) 2. 5
M easurement Results
Weighted GTLS estimator Due t o the equal w eighting over all t he fre
quencies GT L S est imator overemphasize t he high frequency. T he efficiency may be poor in some
flig ht test ing dat a. T he test w as performed using linear sweep ex citat ion vary ing f rom 1 H z t o 20 Hz, sampled at 256 Hz. T he dat a w as measured using tw o channels corresponding t o the force and the accelerat ion response, respectively. T he t rans
case. Adding an appropriat e f requency dependent w eight ing is the key solution to improve its eff i
f er f unction is to be est imated in t he band [ 8 Hz,
ciency [ 9] . A possible w ay to im prove t he eff iciency
14 H z] , because t he dominant modes of aircraft
consists in providing each row of A by a w eighting
most ly dist ribute in t he frequency band. F ig. 4 shows the com parison of original noisy
derived f rom M L w eighting . T his results in f ollow ing cost funct ion F
∗
K W GTL S ( p, Z m) =
W -ML2 ( j
k,
p ) | ∗( j
k)
|2
k= 1 F
∗ W -ML2 ( j
k,
p ) W 2M L( j
k)
k= 1
W -ML2 (
Note t hat
( 18) j k , p ) depends on t he pa
ramet er p, in analog y w ith iterat ive linear least [ 10]
squares met hod
, and an it erative procedure can
response and cleaned response by T ime frequency filt ering.
T hen,
f requency response funct ions
( FRF) est im at ed by cleaned data ( solid line) and original dat a ( dott ed line ) are show n in Fig. 5. Not e t he resonance third peak at 10. 6 H z, w hich is dism issed by FRF derived from t he original noisy dat a. Comparing w it h the noisy data, t he improve ment is obvious w it h the t ime f requency filtering dat a in ident ifying modal peaks.
∋ 50 ∋
T A NG Wei, SHI Zho ng ke, L I Hong chao
( a) orig inal noisy response
( a)
CJA
clean F RF ( dotted line) and the estimated transfer function by W LS ( solid line)
( b) cleaned response Fig . 4
Orig inal noisy response and cleaned r esponse by time frequency filter
( b)
clean FRF ( do tted line) and the estimated transfer function by weig hted GT L S ( solid line)
Fig . 6
Comparisons of est imated transfer function be tween W LS and weighted GT LS in frequency domain
4 Conclusion In this paper, a GT LS estimator in frequency domain combined w it h t ime frequency filtering is presented for improvement of modal paramet er est i mat ion. T he real ex ample results show that t he fil t er in t ime f requency domain can reduce t he non Fig . 5
FRF derived from clean data ( solid line ) and noisy data ( dotted line)
A w eight ed GT L S estimat or is implement ed t o model t he F RF derived from clean dat a w ith a ra
stationary noise in testing dat a signif icantly, and GT L S based on stochast ic model allows us t o obt ain an accurate est imat ion of transfer funct ion in modal paramet er ident ificat ion. References
t ional form n = 21, d = 21. T he w eighted least square est imator ( WLS) has been added for com parison. T he result is show n in F ig. 6. Ag reement is consistently bett er betw een clean FRF and
[ 1]
Brenner J M , Lind C R , V oracek F D. O verview of recent f light flut t er t est ing research at NA SA dryden [ R ] . A IAA 97 1023, N A SA TM 4793, 1997.
[ 2]
Ghiringhell L G , Lanz M , M ant egazza P. A comparison of
WGT L S identification result s, as compared to t he
methods used for t he ident if icat ion of f lut t er f rom experimen
WL S t ransf er funct ion identification.
t al data[ J] . Journal of Sound and V ibrat ion, 1987, 119( 1 ) : 39- 51.
Frequency Domain GT LS Identif icat ion Com bined w ith T ime Frequency February 2006 [ 3]
Filt ering f or Flight Flut ter M odal Paramet er Identif icat ion
Guillaume P, Pint elon R, Schoukens J. D escription of a para
ratio of tw o complex polynomials [ J] . IEE E Transaction on
met ric maximum likelihood est imator in the frequen cy domain
A utomat ic Cont rol, 1963, 9( 1) : 56- 58.
f or multi input , mult i output syst em and its applicat ion t o
Biographies:
f light flut ter analysis[ J] . M echanical Syst em and Signal Pro
TANG wei Born in 1977, male, a native o f XiangF an, HuBei, he received his B. S. and M . S. in 2000 and 2003 fro m N orthwestern Polytechnical U niversity respectively. He is pr esently a Ph. D. candidate of automat ic contr ol theory and automation engineering. His research inter est is system identification and signal processing. T el: ( 029) 88494465, E mail: attracker @ 163. com SHI Zhong ke Born in 1956, he received B. S. from Northwestern Polytechnical U niv ersity in 1981. He receiv ed his M . S. and Ph. D. in 1988 and 1994 r espectiv ely. He has published over hundr ed scientific papers in various periodi cals. He is currently a professor at Northwestern Polytechni cal U niversity of automation. T el: ( 029) 88494465, 88495823, E mail: zkeshi@ nwpu. edu. cn LI Hong chao Bo rn in 1965, he received B. S. and M . S. fro m N orthwestern Poly technical U niversity in 1985 and 1988 respectively . He is presently a Ph. D. candidate of sys tem engineering. His research interest is robust control. T el: ( 029) 88494465, E mail: hongcli@ 163. com
cessing, 1990, 4( 5) : 405- 416. [ 4]
T ang W , S hi Z K . Wavelet denoising of flight f lut t er t est ing dat a f or improvement of parameter ident ificat ion[ J] . Chinese Journal of A eronautics, 2005, 18( 1) : 72- 77.
[ 5]
Daubechies I. Ten l ect ures on wavelet s[ M ] . Philadelphia: Societ y for Indust rial an d Applied M at hemat ics, 1992.
[ 6]
V anderst een G , van Hamme H, Pint elon R . G eneral f rame w ork for asympt ot ic properties of generaliz ed w eighted non linear least squares est imat ors w ith det erminist ic and stochas t ic w eight ing[ J ] . IEEE T ransact ion on A ut omat ic Cont rol, 1996, 41( 10) : 1501- 1507.
[ 7]
van Huf f el S, V andewalle J. T he t ot al least squares problem : computational aspect s and analysis, in front iers in applied mathemat ics, V ol. 9 [ M ] . Philadelphia: Societ y for Indus t rial and A pplied M at hematics, 1991.
[ 8]
Gene H G , van Loan C. M at rix computations[ M ] . T hird Edit ion. Baltimore: Johns Hopkins U niversity Press, 1996.
[ 9]
Pint elon R, Guillaume P, R olain Y , et al . Paramet ric iden t ification of transf er f unct ions in t he frequency domain , a survey[ J] . IEEE Transaction on A utomat ic Control, 1994, 39( 11) : 2245- 2260.
[ 10]
∋ 51 ∋
Sanat hanan K C, Koemer J. Transf er funct ion syt hesis as a
No. 1
CH INES E JO U RN A L O F AER ON A U TICS
February 2006
Frequency Domain GTLS Identification Combined with Time Frequency Filtering for Flight Flutter Modal Parameter Identification T ANG Wei, SHI Zhong ke, LI Hong chao ( Coll ege of A utomat ion , Northw ester n Polytechnical Uni versity , X i an Abstract:
710072 , China)
T he aim of this paper is to present a new method for flig ht flutter modal parameter identifi
cat ion in noisy environment. T his method employs a time frequency ( T F ) filter to reduce the noise be fore identificatio n, w hich depends on the lo calization pr operty of sweep ex citation in T F domain. T hen, a gener alized total least square ( G T LS) identification algo rithm based o n stochastic framew ork is applied to t he enhanced data. System identificatio n w ith no isy data is transformed into a generalized total least square problem, and the solution is car ried out by the g eneralized sing ular value decomposition ( GSVD) to avoid the intensive nonlinear optimization co mputation. A nearly max imum likelihood pro perty can be achiev ed by optimally! weig hted g eneralized tot al least square. F inally, t he efficiency of the method is illustrated by means of flight test data. Key words:
flutter ; identification; time frequency filtering; GT L S
基于时频 域滤波及频域广义整体最小二乘辨识的飞机颤振模态参数辨识. 唐炜, 史忠科, 李洪超. 中国航空学报( 英文版) , 2006, 19( 1) : 44- 51. 摘
要: 提出一套适用于 噪声环 境的 飞机颤 振模 态参 数辨识 方法。为 减小 噪声 对辨 识结 果的 影
响, 首先设计了一种针对扫频激励的 时频滤 波器, 利用扫 频信号 及其响 应在时 频域分布 较为集 中 的特点, 有效去除噪声, 提高了试验数据的信噪比。为进一步提高 辨识精度, 提出了 一种基于 随机 模型的频域广义最 小二乘辨识算法。将 噪声条件 下的系 统辨识 问题转 化为广 义整体最 小二乘 问 题, 并采用线性的广义奇异值分解 求解模 型系数, 避免 了非线 性优化 的复杂 计算。通过 优化加 权 项, 获得了接近极大似然估计的辨识效果。最后, 通过试飞试验数据验证了方法的有效性。 关键词: 颤振; 辨识; 时频域滤波; 广义整体最小二乘 文章编号: 1000 9361( 2006) 01 0044 08
中图分类号: T N911. 7; N945. 14
文献标识码: A
New or modif ied aircraft of ten necessitates flig ht flut ter test s to verify the saf et y marg ins and
curat e results, w hen applied for flut ter modal anal ysis[ 2] . Recent ly , t he max imum likelihood estima
to prevent the cat ast rophic f lutt er. F light flut ter
t or ( ML E) has been developed and implemented for flutt er modal parameter ident ificat ion success
test ing cont inues to be a challenging research area because of t he concerns w ith cost, t ime and safet y in expanding the envelope of new or modif ied air craft . Hence, the procedures of flight t est s are de sired t o dramatically reduce t he flig ht flut ter t est time, increase the accuracy and reliability of t he [ 1]
modal param eter estimat ion . Aircraft in f light test ing dat a are typically
[ 3]
f ully . However, minimizat ion of m ax imum like lihood nonlinear cost funct ion requires both a pow erful iterat ive solver and a set of st art ing values that are ∀ close enough# t o global minimum. Conse quent ly, t here is no guarant ee at all that the opt i m ization algorit hm converges t o global minimum of ML cost funct ion.
charact erized by high noise level due t o at mospheric turbulence. For t his reason, the t radit ional param
T his paper focuses on invest ig at ing a new met hod f or f light flut ter modal paramet er ident if i
eter identification methods generally produce inac
cat ion in noisy environment. A generalized tot al
R eceived dat e: 2005 06 27; R evision received dat e: 2005 11 23 Foundation it em: N at ional N at ural S cience Foundat ion of China ( 60134010)
Frequency Domain GT LS Identif icat ion Com bined w ith T ime Frequency February 2006
∋ 45 ∋
Filt ering f or Flight Flut ter M odal Paramet er Identif icat ion 2
modal identif icat ion. Comparing w ith t he w avelet
1 - t2 - j 0| t| e e ( 2) 2 T he Fourier spect rum of M orlet w avelet is a
denoisng procedure for impulse ex citat ion discussed in Ref. [ 4] , a more pow erf ul noise reduct ion tech
shift ed Gaussian funct ion, and can be seen as a narrow band filt er w it h central frequency f 0
nique using t ime f requency f ilt er for linear sw eep
!(f ) = 2 e- 2 (f - f 0 ) ( 3) In order t o avoid phase dist ortion in signal re
least square ( GT L S ) estimator combined w ith time f requency f ilt ering is presented for accurate
excit at ion w ill be described in Sect ion 1. T hen, w it h the enhanced dat a, t he GT LS est imat or based on st ochastic f ramework will provide accurate iden t ification. Section 2 gives the detailed descript ion of t he GT LS algorit hm. Finally, t he met hod w ill be validat ed by real flig ht t est dat a.
1
T ime frequency Filtering for Noise Re duct ion Sw eep ( chirp) is t he most popular form of
( t) =
2
2
construction, only real M orlet w avelet is used 2 1 - t2 ( t) = e cos2 f 0 t ( 4) 2 T he localizat ion of M orlet in t ime and frequen cy domain allows for mapping t he signal into T F domain. By adjust ing t he paramet er f 0 , it is possi ble t o generate a family of narrowband filters w ith arbit rary cent ral f requency f i ( i = 1, 2, &, n ) . T he cont inuous w avelet transform maps t he one di
cont inuous excit at ion used in vibrat ion test , and is part icular com mon in flight flut ter test ing. T he
mensional sweep signal int o tw o dimension picture index ed by time ( b ) and frequency ( f ) . A t ime
part icular concent rat ion propert y of sweep sig nal in
frequency representat ion of signal is oft en called a scalogram, w hich is act ually t he pow er spect ral
time frequency ( T F ) domain allows for designing a
densit y | W ( b , f ) | of the signal over t he ( b , f )
filt er in T F f or response. 1. 1 Time frequency analysis and CWT T here have been many T F analysis tech
plane. 1. 2 TF filtering for sweep response
niques, such as spectrograms ( short time F ourier
T he t ypical form of linear sw eep signal is
transforms) , scalograms ( w avelets) and bilinear T F distribut ions in the Cohen! s class. Som e of
( 5) f 0 t + 1 r t2 2 where f 0 is the init ial f requency, r is t he frequency
them, such as spectrograms have a f ixed T F reso lut ion, the bilinear T F dist ributions have a high resolut ion but have crosst erms for mult icomponent signal. T he linear t echniques ( f or ex am ple w avelet transform) don! t have crossterms, and also have a g ood resolut ion by adjust ing dilat ion.
Hence,
w avelet transform is select ed as t he T F analysis tool in t his sect ion.
x ( t ) = A ( t ) cos
2
change rate. F ig. 1 show s an ex ample of scalogram of a 0 10 H z sw eep ex citation signal and m easured ac celeromet er response in flight . Obviously, t he sw eep show s a w ell concentrated propert y in T F, it localizes on a small reg ion. Applying sweep as an input to linear t im e invariant ( LT I) syst em, t he
F or a signal x ( t ) , the cont inuous w avelet transform can be defined as W( a, b) =
1 a
+ %
∃ x(t) - %
*
t - b dt ( 1) a
w here b is a t ranslation indicat ing t he localit y, a is a dilat ion or scale paramet er. ( t ) is a basic w avelet and
*
( t ) is the complex conjugat e of
( t) . One of t he most widely used funct ion in T F analysis is the Morlet w avelet[ 5] ,
( a)
Scalogram of sweep excitation sig nal
∋ 46 ∋
T A NG Wei, SHI Zho ng ke, L I Hong chao
CJA
D( b i , f j ) =
1
if | We ( b i , f j ) | > t 0
0
otherw ise
i = 1, &, N ; j = 1, &, M where D is the mask mat rix f rom w avelet coeff i cient We( b i , f j ) of sw eep ex citation by t hreshold ing , t 0 is t hreshold value. T hen, a mask operat ion is performed on wavelet coeff icients of noisy response, W ^ r ( b i , f j ) = D( b i , f j ) ∋ Wr ( b i , f j ) ( b)
Scalogram of accelerometer response
Fig . 1 Scalograms of flight flutter testing data
where W r denot es t he coefficients mat rix of noisy response and W ^ r represent s the coeff icients m at rix
response usually holds t he propert y and concen
aft er mask operat ion.
t rat es along sw eep line in T F. Hence, t he darker
1. 4
sloping area in ( b) mainly represent s t he true re sponse; the others in T F represent noise. If it is
T F domain t o t ime domain. H owever, t he f ilter
possible to locate t he region where t he chirp signal concent rat es in T F plane, a mask operat ion can be preformed to ext ract the t rue response. T he t ime frequency analysis sugg ests the f ollowing basic fil t ering procedure for response of sw eep ex citat ion, as shown in Fig. 2. ( 1) Perform t ime f requency analysis for sw eep and it! s response.
Response reconstruction Inverse t ransf orms w ill map t he sig nal from
bank reconstruction for discret e wavelet t ransform ( DWT ) can not be applied to CWT . A novel re construction in f requency domain is employed in this sect ion. Before int roducing t he met hod, express t he reconst ruct ion problem into a L S form, the f ilt ered response ^s is t hat minimizes
( 2) Remove noise component s in scalogram
(H^s - D ∋ * Hf (F = min (Hs - D ∋ * Hf (F = s
using mask operation. ( 3 ) Reconstruct t he response using clean
where operator H denot es t he w avelet t ransf orms.
scalogram.
f is the noisy response. Arithmet ic operator ( *∋ )
min ( Ws - D∋ * Wf ( F s
denot es the matrices element by element mult iplica tion. If H is a matrix, the normal solut ion of LS is s = H+ D ∋ * Hf
( 6)
where ∀ + # denot es pseudo inverse. How ever, t he comput at ion of CWT w ith convolut ion can not be represent ed by mat rix directly, + %
∃ x(t) !
1 a
Wx ( b, f 0 ) =
- %
* f0
b- t dt = a
x ( t ) * ! *f 0 , a ( t )
Fig . 2 T F filtering pr ocedure in T F plane
( 7)
Instead of Eq. ( 7) , CWT can be ex pressed in 1. 3
Mask operation Since t he input sw eep signal is known before,
frequency domain, Wx ( , f 0 ) = X( ) ∋ ! *f 0 , a ( )
( 8)
it s pat tern in T F domain may help in designing mask in T F domain for f ilt ering noise. T he com
where Wx (
mon met hod for mask desig n is t hresholding t he
Fourier t ransforms of Wx ( b, f 0 ) , x ( t ) and ! *f 0 , a
modulus of coef ficients of clean sw eep input in T F domain . T he procedure is g iven as
( t ) , respect ively.
, f 0 ) , X ( ) and !
* f , a( 0
) are t he
Write in matrix form in f requency domain,
Frequency Domain GT LS Identif icat ion Com bined w ith T ime Frequency February 2006
∋ 47 ∋
Filt ering f or Flight Flut ter M odal Paramet er Identif icat ion
F Wx = ∀ ∋ X
shown in F ig. 3. T he m easured input Um and out put Ym are dist urbed w it h noise. U0 and Y0 re spect ively are the t rue sw eep and its accelerometer response f ree noise, respectively. M u and M y are
w here Wx (
1,
f 1)
&
f N)
# &
#
F Wx = Wx (
1,
X = diag[ X( *
! f 1, a ( ∀ =
1) ,
# !
* f
N, a
(
Wx (
&
M,
M)]
*
! f 1, a (
& 1)
M,
f 1)
#
&, X(
&
1)
Wx (
the measurem ent noise in input and output, respec t ively. N g denot es t he at mospheric turbulence.
fN)
;
M
)
# ! *f
N, a
(
M)
Hence, LS can be w ritt en as in term of t he F ourier t ransf orm msin (F Ws - F ^ Wf (F = min ( ∀ ∋S- F ^ Wf (F s w here S = diag [ S (
1) ,
&, S (
M)
], S (
i)
Fig. 3
is
the discrete Fourier transforms ( DF T ) of estimat ed signal s . F ^ Wf denot es t he Fourier spect ra matrix of coeff icient s aft er m ask operation. T he minimum problem can be solved by t he
T he measurement s Um and Y m are related t he exact values U 0 and Y 0 by Um ( j ) = U0 ( j ) + N U ( j ) Ym ( j ) = Y 0 ( j ) + N Y ( j ) where N U ( j ) = MU ( j )
follow ing equat ion ∀ ∋s = F ^ Wf
( 9)
T he est imat ed solution is ^ = ∀+ F S ^ Wf T hen t he diagonal elem ent s of S ^ are inversely
Stochastic model of flight flutter test
N Y ( j ) = M Y ( j ) + H 0 ( j ) N g( j ) T
Def ine ∃Z = [ N U ( j k ) N Y ( j k ) ] , t he above equation can be w ritt en as Z m ( j k ) = Z 0 ( j k ) + ∃Z( j k ) ( 10) where T
transformed by IDF T ( inverse discret e Fourier tr
Z 0 = [ U0( j
k)
Y0( j
ansf orm ) to get t he t ime domain reconst ruct ion
Z m = [ Um( j
k)
Ym ( j k ) ] T
^s ( t) .
k) ]
It is reasonable t o make the following assump t ion : Assumpt ion 1. T he noise ∃Z ( j k ) is zero [ 6]
2
Generalized T otal Least Squares Estima tor for Ident ificat ion
T he common approach for modal ident if icat ion st art s from frequency response funct ions ( FRF ) ,
mean, complex normally distribut ed wit h covari ance m at rix E( ∃Z ( j k ) ∃Z H ( j k ) ) = Cz ( j k ) = %2U ( j
w hich are derived using a non parametric FRF est i
2 %Y U ( j
mator such as H 1 or H v . How ever, t he perf or mance of FRF estimator is sensitive to noise. Even w it h t he enhanced dat a, t he t radit ional frequency ident ification algorit hms perform poorly. For t his reason, the GT LS identificat ion algorit hm based on st ochast ic framew ork is applied, w hich improve t he accuracy of the paramet er est imation. 2. 1 Descri ption of flight flutter test with stochastic framework T he stochast ic model of flight f lutt er t est is
k)
E( ∃Z ( j
k)
k) T
%2U Y ( j 2 %Y (
j
k) k)
∃Z ( j l ) ) = 0
k, l
Assumpt ion 2. T he noise at each dif ference frequency is independent ly distributed, E ( ∃Z ( j 2. 2
k)
H
∃Z ( j l ) ) = 0,
k ) l
Parametric model A rational transfer funct ion model of order n/
d is used for frequency domain estimat or, N ( j , &) H 0 ( j , p) = D ( j , ∋) =
∋ 48 ∋
T A NG Wei, SHI Zho ng ke, L I Hong chao
&n ( j ) n + &+ &1 ( j ) + &0 d ∋d ( j ) + &+ ∋1 ( j ) + ∋0 w it h t he real coef ficients p = [ ∋0 & ∋d &
CJA
where A may cont ain errors in all elements.
( 11)
T he generalized tot al least squares ( GT L S)
&0
solut ion to the est imat ion in Eq. ( 14) is given
&n ] . Define p r ( r = 1, 2, &, d ) to be the poles of
by
[ 7]
arg min (( A - A ^ ) C- 1 (2F
the transfer funct ion, t he corresponding modal fre
( 15)
^ ,p A
quency and damping rat io can be obt ained as Im( p r ) Re( p r ) fr = , (r = 2 | pr | Obviously , t he modal parameters are derived
and subjected to A ^ p = 0, p T p = 1, w here C is a square root of t he colum n covariance matrix of A: H
H
C C= E{ ∃A ∃A} . A ^ is t he estimate of A. Eliminat ion of A ^ in Eq. ( 15) g ives the e
from the transfer function model, hence it is im
qualient cost funct ion minimized by GT L S estima
port ant to focus on estimat ing coef ficients of t rans
t or,
f er funct ion. 2. 3 GTLS in frequency domain
p T A H Ap T H p C Cp
arg pmin
If there are no model errors, t he out put and
subject to p T p = 1 ( 16)
input obey t he following system relat ion Y0 ( j ) = H 0 ( j ) U0 ( j ) U sing Eq. ( 11) , t he model equation is readily
p T A H Ap = ), and Eq. ( 16) p T CH Cp can be t ransformed into t he follow ing equation, Define
obtained,
min
Y0 ( j ) D( j , ∋) - U 0 ( j ) N ( j , &) = 0
AT Ap = )CT Cp
T he above equat ion can be w ritt en in m at rix
T he solut ion p t o the equat ion can be obtained by generalized sing ular value decom posit ion ( GSVD) of matrix pair ( A and C) [ 8] . T he vect or
form A0 p = 0
( 12)
in right m at rix corresponding t o minimum general
w here A0 = [
p Tp = 1
subject to
a H01
&
aH02
aH0F ] H ,
a 0k =
Z T0 ( j k )
ized singular value is the solut ion. Comparing w ith ML , this w ill avoid t he intensive nonlinear iterat ive
S( j k ) ,
S( j k ) = block diag(- [ 1, &, ( j ) ] , [ 1, &, ( j ) ] ) n
d
optimizat ion and relies on t he init ial point .
Put t ing the noisy value Eq. ( 10) into Eq. ( 12)
In order t o simplif y Eq. ( 15) , define
def ine the noisy matrix A, A = A 0 + ∃A
∗( j
( 13)
= a kp = Y( j
k)
U( j
w here ∃A = [ ∃aH1
∃aH2
2 WML ( j
& ∃aHF ] H,
k)
N (j
p ) = E { ∃∗( j
k,
E { ∃akpp
T
∃∋k = ∃Z ( j k ) S( j k ) . %2Y
U sing Eq. ( 13) , t he paramet er estimation can be formulat ed as a tot al least squares, looking for a
p)
+ k,
D( j
k,
p)
k)
∃a Hk }
%2U
k,
∃∗( j
p) k)
H
}=
=
N (j
p) D ( j
k, k,
p)
2
-
p) H )
T he cost f unction for minimizat ion in Eq. ( 15) can be rewritt en as
( 14)
F
∗
K GT LS ( p, Z m ) =
k,
2Re( %2U YN ( j
solut ion of Ap = 0
D( j
2
H
k)
∗( j
k)
2
k= 1 F
∗
W ML ( j
k)
2
=
k= 1 F
∗
Y( j
k)
D( j
k,
p ) - U( j
2
+ %U N ( j
2
k,
p)
k) N
(j
k,
p)
2
k= 1 F
∗
k= 1
2
{ %Y D ( j
k,
p)
2
( 17) 2
H
- 2Re( % UYN ( j k , p ) D( j k , p ) ) }
Frequency Domain GT LS Identif icat ion Com bined w ith T ime Frequency February 2006
2. 4
∋ 49 ∋
Filt ering f or Flight Flut ter M odal Paramet er Identif icat ion
Sample noise variance estimation
be used. T he parameter est imated f rom previous
Before using Eq. ( 17) t o est imat e the coeff i cients, it is needed t o est imat e the noise covariance and t he covariance f rom sample data. T he sample noise variance is usually derived from a set of data according t o st at ist ical definit ion, w hich is hardly
step is used to updat e the nex t st ep, and t he cost funct ion is w rit ten as F
K WGTLS ( pi , Z m ) =
∗
W-ML2 ( j k , p i- 1 ) | ∗( j k , pi ) | 2
∗
WML( j k , p i- 1 ) WML( j k, pi )
k= 1 F
-2
2
k= 1
to implement for limit ed t im e in f light flut ter test. F or this reason, a method to est imate the variance
( 19)
from one set of measurements is present ed. T he in
where K W GTLS ( pi , Z m ) denotes the cost funct ion
put sweep signal for excit at ion is know n before,
at i t h iterat ion step. p i- 1 is t he est imat ion derived
w hich is free noise. T o simplif y t he problem, only the error on the measured out put sequences will be
from previous step.
considered w ith %2U = 0, %2UY = 0.
reduced to M L cost funct ion
T he noise variance on t he out put is %2Y ( ) = E { N Y ( ) N HY ( ) } = H
K WGT LS, i = K ML / F
H
where K M L is the cost funct ion of M L est imat or.
E { N Y ( ) [ Ym ( ) - Y0 ( ) ] } =
T he f act t hat WGT L S and ML cost funct ion are e
H
E { N Y ( ) Ym ( ) } = H
quivalent indicat es t hat t he optimally! w eighted
H
E { Ym Y m - Y 0 Y m} = H
WGT L S has nearly M L property. Hence, t he
H
E ( Y m Ym ) - H 0 ( j ) E ( U0 Ym ) T he dat a are divided int o M segm ent s f or H 1 est imator, and H 0 is replaced by est imat ion H ^ 0, M
^%2Y (
∗
)=
Ym ( j )
Considering t he convergences of est imator, as the it eration i + % : p i - 1 = pi . Eq. ( 19) is readily
WGT L S can be used to generate the st art ing values that are ∀ close enough# t o t he act ual minimum for ML .
)- H ^ 0( j ) ∋
YHm ( j
3
i= 1 M
∗
H
U0 ( j ) Ym ( j ) / M = S Y Y( j ) -
In this sect ion, t he GT LS est imat or combin ing w ith time frequency filtering will be applied t o
i= 1
H ^ 0 ( j ) SUY ( j ) 2. 5
M easurement Results
Weighted GTLS estimator Due t o the equal w eighting over all t he fre
quencies GT L S est imator overemphasize t he high frequency. T he efficiency may be poor in some
flig ht test ing dat a. T he test w as performed using linear sweep ex citat ion vary ing f rom 1 H z t o 20 Hz, sampled at 256 Hz. T he dat a w as measured using tw o channels corresponding t o the force and the accelerat ion response, respectively. T he t rans
case. Adding an appropriat e f requency dependent w eight ing is the key solution to improve its eff i
f er f unction is to be est imated in t he band [ 8 Hz,
ciency [ 9] . A possible w ay to im prove t he eff iciency
14 H z] , because t he dominant modes of aircraft
consists in providing each row of A by a w eighting
most ly dist ribute in t he frequency band. F ig. 4 shows the com parison of original noisy
derived f rom M L w eighting . T his results in f ollow ing cost funct ion F
∗
K W GTL S ( p, Z m) =
W -ML2 ( j
k,
p ) | ∗( j
k)
|2
k= 1 F
∗ W -ML2 ( j
k,
p ) W 2M L( j
k)
k= 1
W -ML2 (
Note t hat
( 18) j k , p ) depends on t he pa
ramet er p, in analog y w ith iterat ive linear least [ 10]
squares met hod
, and an it erative procedure can
response and cleaned response by T ime frequency filt ering.
T hen,
f requency response funct ions
( FRF) est im at ed by cleaned data ( solid line) and original dat a ( dott ed line ) are show n in Fig. 5. Not e t he resonance third peak at 10. 6 H z, w hich is dism issed by FRF derived from t he original noisy dat a. Comparing w it h the noisy data, t he improve ment is obvious w it h the t ime f requency filtering dat a in ident ifying modal peaks.
∋ 50 ∋
T A NG Wei, SHI Zho ng ke, L I Hong chao
( a) orig inal noisy response
( a)
CJA
clean F RF ( dotted line) and the estimated transfer function by W LS ( solid line)
( b) cleaned response Fig . 4
Orig inal noisy response and cleaned r esponse by time frequency filter
( b)
clean FRF ( do tted line) and the estimated transfer function by weig hted GT L S ( solid line)
Fig . 6
Comparisons of est imated transfer function be tween W LS and weighted GT LS in frequency domain
4 Conclusion In this paper, a GT LS estimator in frequency domain combined w it h t ime frequency filtering is presented for improvement of modal paramet er est i mat ion. T he real ex ample results show that t he fil t er in t ime f requency domain can reduce t he non Fig . 5
FRF derived from clean data ( solid line ) and noisy data ( dotted line)
A w eight ed GT L S estimat or is implement ed t o model t he F RF derived from clean dat a w ith a ra
stationary noise in testing dat a signif icantly, and GT L S based on stochast ic model allows us t o obt ain an accurate est imat ion of transfer funct ion in modal paramet er ident ificat ion. References
t ional form n = 21, d = 21. T he w eighted least square est imator ( WLS) has been added for com parison. T he result is show n in F ig. 6. Ag reement is consistently bett er betw een clean FRF and
[ 1]
Brenner J M , Lind C R , V oracek F D. O verview of recent f light flut t er t est ing research at NA SA dryden [ R ] . A IAA 97 1023, N A SA TM 4793, 1997.
[ 2]
Ghiringhell L G , Lanz M , M ant egazza P. A comparison of
WGT L S identification result s, as compared to t he
methods used for t he ident if icat ion of f lut t er f rom experimen
WL S t ransf er funct ion identification.
t al data[ J] . Journal of Sound and V ibrat ion, 1987, 119( 1 ) : 39- 51.
Frequency Domain GT LS Identif icat ion Com bined w ith T ime Frequency February 2006 [ 3]
Filt ering f or Flight Flut ter M odal Paramet er Identif icat ion
Guillaume P, Pint elon R, Schoukens J. D escription of a para
ratio of tw o complex polynomials [ J] . IEE E Transaction on
met ric maximum likelihood est imator in the frequen cy domain
A utomat ic Cont rol, 1963, 9( 1) : 56- 58.
f or multi input , mult i output syst em and its applicat ion t o
Biographies:
f light flut ter analysis[ J] . M echanical Syst em and Signal Pro
TANG wei Born in 1977, male, a native o f XiangF an, HuBei, he received his B. S. and M . S. in 2000 and 2003 fro m N orthwestern Polytechnical U niversity respectively. He is pr esently a Ph. D. candidate of automat ic contr ol theory and automation engineering. His research inter est is system identification and signal processing. T el: ( 029) 88494465, E mail: attracker @ 163. com SHI Zhong ke Born in 1956, he received B. S. from Northwestern Polytechnical U niv ersity in 1981. He receiv ed his M . S. and Ph. D. in 1988 and 1994 r espectiv ely. He has published over hundr ed scientific papers in various periodi cals. He is currently a professor at Northwestern Polytechni cal U niversity of automation. T el: ( 029) 88494465, 88495823, E mail: zkeshi@ nwpu. edu. cn LI Hong chao Bo rn in 1965, he received B. S. and M . S. fro m N orthwestern Poly technical U niversity in 1985 and 1988 respectively . He is presently a Ph. D. candidate of sys tem engineering. His research interest is robust control. T el: ( 029) 88494465, E mail: hongcli@ 163. com
cessing, 1990, 4( 5) : 405- 416. [ 4]
T ang W , S hi Z K . Wavelet denoising of flight f lut t er t est ing dat a f or improvement of parameter ident ificat ion[ J] . Chinese Journal of A eronautics, 2005, 18( 1) : 72- 77.
[ 5]
Daubechies I. Ten l ect ures on wavelet s[ M ] . Philadelphia: Societ y for Indust rial an d Applied M at hemat ics, 1992.
[ 6]
V anderst een G , van Hamme H, Pint elon R . G eneral f rame w ork for asympt ot ic properties of generaliz ed w eighted non linear least squares est imat ors w ith det erminist ic and stochas t ic w eight ing[ J ] . IEEE T ransact ion on A ut omat ic Cont rol, 1996, 41( 10) : 1501- 1507.
[ 7]
van Huf f el S, V andewalle J. T he t ot al least squares problem : computational aspect s and analysis, in front iers in applied mathemat ics, V ol. 9 [ M ] . Philadelphia: Societ y for Indus t rial and A pplied M at hematics, 1991.
[ 8]
Gene H G , van Loan C. M at rix computations[ M ] . T hird Edit ion. Baltimore: Johns Hopkins U niversity Press, 1996.
[ 9]
Pint elon R, Guillaume P, R olain Y , et al . Paramet ric iden t ification of transf er f unct ions in t he frequency domain , a survey[ J] . IEEE Transaction on A utomat ic Control, 1994, 39( 11) : 2245- 2260.
[ 10]
∋ 51 ∋
Sanat hanan K C, Koemer J. Transf er funct ion syt hesis as a