Modal parameter identification of a CMUT membrane using response data only

Capacitive micromachined ultrasonic transducers (CMUTs) are microelectromechanical systems used for the generation of ultrasounds. The fundamental element of the transducer is a clamped thin metallized membrane that vibrates under voltage variations. To control such oscillations and to optimize its dynamic response it is necessary to know the modal parameters of the membrane such as resonance frequency, damping and stiffness coefficients. The purpose of this work is to identify these parameters using only the time data obtained from the membrane center displacement. Dynamic measurements are conducted in time domain and we use two methods to identify the modal parameters: a subspace method based on an innovation model of the state-space representation and the continuous wavelet transform method based on the use of the ridge of the wavelet transform of the displacement. Experimental results are presented showing the effectiveness of these two procedures in modal parameter identification.


Introduction
Capacitive micromachined ultrasonic transducers (CMUTs) introduced in 1994 constitute a good alternative to conventional piezoelectric transducers in terms of broader bandwidth, better transduction efficiency and the possibility of on-chip integration with electronics components.Basically, a CMUT cell is a capacitor structure composed of a thin flexible membrane, metallized or not metallized, and a bottom fixed plate.If an alternating voltage, applied between the membrane and the back plate, is superimposed on the bias voltage the modulation of the electrostatic force moves the membrane which generates ultrasonic waves.Fig. 1 shows a schematic cross-section of a CMUT elementary cell.A CMUT cell uses the deformation of a suspended thin membrane to transmit ultrasonic waves and the dynamic analysis of a CMUT has been studied in time domain using a FEM package [1,2].They determined the normal displacement at the center of the membrane by finite element analysis and by experimental tests.The idea of generating acoustic waves by the electrostatic attraction force between the plates of a capacitor is very old and a tutorial can be obtained in [3].Dynamic and acoustic modelization of capacitive micromachined ultrasonic transducers have been proposed in [4] where the acoustic pressure emitted by a CMUT is estimated by taking into account explicitly the dynamics of each membrane of the network.

Fig. 1. A conceptual cross-section of a capacitive micromachined ultrasonic transducer
CMUTs are fabricated using silicon processing routes, leading to a thin (0.5 -2 µm thick) silicon nitride or polysilicon membrane above a small gap (typically about one micron in depth) to result in a capacitive structure.In general, such CMUTs tend to be less than 0.2 mm in diameter.We note R the membrane radius, tm the membrane thickness and tg the gap (or cavity) thickness.
A good design of a capacitive ultrasonic sensor requires a large displacement of the membrane and an optimum energy coupling between the membrane and the air that are achieved when the membrane vibrates at its mechanical resonance frequency.Stiffness is an important property and is vital to define both the sensitivity of CMUTs sensors and the dynamic behavior of the device that is strongly dependent on material characteristics and on clamping conditions of the membrane.The quality factor is also fundamental, as it provides immediate quantification of dissipations, which are generally addressable to air damping and internal friction in the material.In the paper the membrane displacement of a CMUT cell is reduced to a single degree of freedom and only the first vibration mode is considered.Higher order modes are neglected in our model and one objective of the paper is to identify the parameters of the first vibration mode from output data only.
Numerous papers have been presented on system parameter identification from measured data.Ho and Kalman [5] introduced the principle of minimal realization theory in which a state-space model has an alternative representation with a sequence of Markov parameters.Zeiger and Mc Ewem [6] proposed a concept of combining singular value decomposition (SVD) and minimal realization for parameter identification.Based on the developments of SVD and minimal realization algorithm, the eigensystem realization algorithm (ERA) was proposed by Juang [7] to identify parameters of a system from the noisy input-output measurements.To reduce the noise effect Juang [7] proposed the eigensystem realization algorithm with data correlation (ERA/DC).James et al. [8] developed the natural excitation technique (NExT) based on the principle that a correlation function measured under natural excitation can be expressed as a sum of decaying sinusoids.Each decaying sinusoid has an eigenfrequency and a damping ratio that is identical to the one of the structural mode.A time domainbased stochastic subspace identification (SSI) method has been proposed by Van Overschee and De Moor [9] and is based on the state-space model of input/output signals and parameter identification is realized by extracting system matrices from input/output signals or output signals only.The key concept of SSI, which is very time consuming, is the projection of the row space of the future time data into the row space of the past time data.The system parameter identification is also possible in a timefrequency domain using the continuous wavelet transform (CWT) of signals [10].The CWT is a timefrequency representation and applied to the free response of a system allows the identification of eigenfrequencies, damping ratios and mode shapes.Indeed, the CWT of a signal results in a complex valued form, whose modulus and phase are related to the parameters of the analysed system [11][12][13][14][15].
In this paper, the dynamic characterization of the metallized membrane is aimed at defining its behavior in working conditions, as well as its performances in terms of resonance frequency, stiffness coefficient, damping factor, quality factor and half-power bandwidth.Such coefficients are identified using two methods: the subspace method derived from a state-space innovation model and the continuous wavelet transform (CWT) method based on the use of the ridge of the CWT.These methods extract the modal parameters from structural response data only, using the time displacement of the membrane center.The procedures presented in the paper can identify closely eigenfrequencies that cannot be identified by the traditional Fourier transform.
The paper is organized as follows: the second section is devoted to the modelization of a vibrating system and to the generation of a discrete time state-space model.An innovation representation with the estimation of the transition matrix is derived in the third section.The continuous wavelet transform is presented in section 4. Validity of the modal parameter identification procedures is presented in the fifth section with simulated and an experimental test in laboratory.The paper is briefly concluded in section 6.

Modeling a vibrating mechanical system
Consider a multivariate continuous time structural dynamic system described by a positive definite diagonal mass matrix Mo, a symmetric semi-definite viscous damping matrix Co and a symmetric positive definite stiffness matrix Ko.These matrices are (n'xn').The system is represented by the linear differential equation and the observation equation [7] Mo   (t) + Co   (t) + Ko  (t) =  (t) (eq. 1) where  (t) is the (n'x1) vector of displacement of degrees of freedom;  (t) is the (n'x1) unmeasured excitation vector; y(t) is the (mx1) output observation vector; the matrix L (mxn') specifies which points of the system are observed, namely where the sensors are located and v(t) is an (mx1) additive noise disturbance, a vector of measurement errors simulated by a gaussian white noise.The modal characteristics (  ,   )of the system are solutions of: where represents a determinant.
The system described by (eq. 1) and (eq.2) is equivalent to the following continuous time state-space model : (eq. 4) y(t) = C x(t) + v(t) (eq.5) where x(t) is the vector of dimension 2n' = n and Ã ,

B
, C are the (2n'x2n'), (2n'xn') ,(mx2n') matrices After sampling with constant period  t and transformation of the 2n' first order differential equations (4) into a discrete time equation, we obtain the following discrete time state-space model and the discrete time observation equation [7,15] xk+1 = A xk + uk (eq.8) yk = C xk + vk (eq.9) where xk represents the unobserved state vector of dimension n; A = e Ã t  is the (nxn) discrete time state-space matrix or discrete time transition matrix and uk is given by uk = e B t s ds The excitation is not measured and can be considered (a priori) as a centred random process.Note that uk and vk are mutually uncorrelated.The matrix C is (mxn) and is called the observation matrix and yk is an (mx1) vector of observations.
The eigenvalues  and eigenvectors   of the discrete time transition matrix A are related to the modal characteristics (3) by [15]  and the global modal parameters: the natural frequencies fi and damping ratios  i of the vibrating system can be determined by fi= for i = 1,2,. ..n' Since the state vector xk of dimension n = 2n' is not observable we estimate it and transform (eq.8) and (eq.9) into a state-space innovation form, before we can obtain the transition matrix A.

State-space innovation model
We estimate xk by its orthogonal projection onto the subspace spanned by the stacked vector of past data vectors noted ] to denote this orthogonal projection.We obtain zk+1 using properties of orthogonal projections given by M. Aoki [16]  we have : (eq. 16) Note that by assumption E[uk y k  Using the definition of ek and (eq.9), we obtain the discrete time observation equation which contains the innovation vector The set of equations (eq.8) and (eq.9) is now replaced by zk+1 = A zk + B ek (eq.18) yk = C zk + ek (eq.19) Such representation is called state-space innovation representation of the data generating process (eq.8)-(eq.9).The state-space innovation model in (eq.18)-(eq.19) is a suitable representation for the synthesis of temporal samples.
Denote the covariance matrix of the stacked vector . This matrix is theoretically infinite-dimensional, although in any actual implementation we cut off the stacked vector at yk-p for some positive integer p.We use the hat to denote an estimate based on finite amount of data.Thus R ˆis (mpxmp).We develop the state vector zk is called the controllability matrix.This is a reduced form expression for the state vector.At the moment equation ( 21) is not operational since we do not know the matrix K, which will be shown to be determined by a suitable factorisation of the block Hankel matrix defined later.Since K it is also a matrix with a theoretically infinite number of columns (eq. 21) where we define the covariance matrix between zk and yk-1 as G = E[zk yk-1'].
To estimate K it is necessary to introduce the stacked vector of future data vectors noted ' , then from equations ( 19) and ( 20) we obtain where the matrix O is the observability matrix of the model ( 19), (20) and has the form which theoretically has an infinite number of rows and F is a lower block triangular matrix formed by the matrices 0m , Im , A , B and C. Then the covariance matrix between the future and the past is The matrix H, called the block Hankel matrix, is formed of (mxm) covariance matrices Ri=E[yk+i y ' k] as submatrices i=1,2... and is factored as the product of two matrices O and K.Note that sample covariance matrices calculated from data yk, k=1,.., A second factorisation of the estimated block Hankel matrix is obtained using the singular value where Û and  V contain the singular vectors of  H and (eq. 30) We have formed the two factors and with generalised inverses Other factorisations are possible, however the representation that corresponds to the choice of (equ.31) and (eq.32) is called balanced since the gramians are equal, we have The unknown transition matrix A may now be found in terms of the known estimated covariance matrices  R i with i = 0, 1, . .and the singular value decomposition of n H ˆ .The procedure is described below.Multiply (eq.18) from the right by z ' k and take expectation to obtain Dropping the last term since ek and zk are uncorrelated theoretically, we obtain From (eq. 20) we have zk = K R--1  1 k y and we obtain where From (eq. 36) we get and finally the transition matrix is given by Note that the controllability matrix K can be estimated from (eq. 26) using the generalised inverse of With estimates of the transition matrix A in hand we compute its eigenvalues and we can identify eigenfrequencies and damping factors of the vibrating system following (eq.12) and (eq.13).However, all the subspace modal identification algorithms have a serious problem of model order determination.When extracting physical or structural modes, subspace algorithms always generate spurious or computational modes to account for unwanted effects such as noise, leakage, residuals, nonlinearity's … For these reasons, the assumed number of modes, or model order, is incremented over a wide range of values and we plot the stability diagram [15].The stability diagram involves tracking the estimates of eigenfrequencies and damping factors as a function of model order.As the model order is increased, more and more modal frequencies and damping ratios are estimated, hopefully, the estimates of the physical modal parameters are stabilized using a criterion based on the modal coherence of measured modes and identified modes [15].Using this criterion we detect and remove the spurious modes.A numerical example and an experimental tests in laboratory showing the effectiveness of the method in modal parameter identification are presented in section 5.

Definition and properties
The continuous wavelet transform (CWT) gives time and frequency information about the analyzed data.For all functions y(t) satisfying the condition dt y(t) 2       <  , which implies that y(t) decays to zero at   , the wavelet transform of y(t) is defined as [10][11][12][13][14][15] where ψ(t) is an analyzing function called mother wavelet, a is the dilatation or scale parameter defining the analysing window stretching and b is the translation parameter localising the wavelet function in the time domain.The CWT represents the correlation between the signal y(t) and a scaled version of the function ψ(t) and the idea of the CWT is to decompose a signal y(t) into wavelet coefficients using the basis of wavelet functions.The decomposition is obtained locally at different time windows and frequency bands.Once the mother wavelet is chosen, the centre of the time window controlled by the translation parameter b while the length of the frequency band is controlled by the dilatation parameter a.Hence, one can examine the signal at different time windows and frequency bands by controlling translation and dilatation.

Modal parameter identification by the CWT
The CWT gives time and frequency information about the analysed data.Any function ψ (t) can be used as an analysing wavelet when it satisfies the admissibility condition [10], this condition being necessary to obtain the inverse of the CWT.There are a number of different complex and real valued functions used as analysing wavelets, but the Morlet wavelet may be the most commonly considered and will be used in this communication.The Morlet wavelet is defined in the time domain as and the maximum is located at a ω = o ω .
Consider a modulated signal in amplitude and frequency: where /db is the instantaneous frequency of the signal.The dilatation parameter is Now, consider the impulse response of a viscously damped single degree of freedom system where is the damped angular frequency and ξ the damping factor.From (eq. 48) we can deduce The phase function is a straight line whose slope represents the damped angular frequency.The log of the wavelet modulus cross-section is again a straight line whose slope is the decay rate n ω ξ .Once the instantaneous frequency and the decay rate have been estimated we can identify the natural frequency and the damping ratio of the vibrating system.

Simulated data
To prove the effectiveness of the identification procedure based on the subspace analysis we consider a two-DOF system with very closely spaced modes.The parameters of the system are f1 = 2 MHz , f2 = 2,01 MHz, 1 ξ = 0,01 and 2 ξ = 0,02.Fig. 2 (a) shows the free response of the system where a Gaussian white noise has been added: the generated data were corrupted by a random noise.The sampling frequency is 5 MHz and only 151 time samples are used in the simulation.The signal duration is 30 s μ .We convert to the frequency domain this time response by taking the discrete Fourier transform of the noisy signal.It is impossible to identify the two frequencies components by using the FFT as shown in Fig. 2(b) and (c), where the power spectral density has been plotted.In our identification procedure by subspaces, we plot the stabilization diagram on eigenfrequencies and damping factors.Fig. 3 shows stabilization diagrams with the modal coherence indicator: spurious modes have been eliminated and only physical modes are present.We apply now the continuous wavelet transform to our numerical test.Fig. 4 presents the continuous wavelet transform modulus and its ridge for the two-DOF system.The CWT is a time-frequency representation and the two modes appear after 18 s μ of signal recording.The identified first mode is plot in red dash line and the identified second mode is plot in magenta dot-dash-line.The two vertical blue lines correspond to the limitations of edge effects.The identified modal parameters for the two dof system 2 ξ ˆ= 0,49.A very satisfactory estimation on eigenfrequencies are obtained, however damping factors are not well identified.This is due to the fact that the two modes are very close and only 151 time samples are used in the analysis.It is shown in [14] that it is necessary to have a long recording time and an important number of time samples to identify correctly closely spaced modes.Furthermore, the vibrating system is very quickly attenuated: the time response decreases rapidly and, it is shown in [14], the continuous wavelet transform gives good results if the time response is not attenuated abruptly.The limitations of the continuous wavelet transform concern the low number of time samples used in the identification procedure and the fast decay of the time response.

Fig. 4. Continuous wavelet transform modulus and its ridge for the two-DOF system
An experimental test in laboratory is presented in the following section.

Modal parameter identification of a CMUT membrane
A CMUT cell is basically a metallized membrane suspended over an electrode and by applying an alternating voltage to the parallel plate capacitor, an electrostatic force will cause the deflection of the membrane and the generation of ultrasounds.The main purpose is to find the first resonance frequency of the membrane, the effectiveness of a CMUT cell being analyzed in terms of resonance frequency, stiffness and damping coefficients.It can be reminded that vibration amplitude of the membrane is maximized at the mechanical resonance frequency.A CMUT is commonly fabricated by means of the surface micromachining technology, using standard integrated circuits techniques.Several processes have been reported in the literature to fabricate CMUTs, using different materials and film deposition techniques.In our experimental test the CMUTs have been fabricated using the anodic bonding technology of a fixed thickness monocrystalline silicon layer of a SOI wafer (Silicon On Insulator technology) on a borosilicate glass substrate.The process evolution for each technological step and the fabrication runs are not discussed, but details can be found in Bellaredj [17].Two photographs representing a circular CMUT cell and an array of CMUTs are presented in Fig. 5.The bonding technology is a particularly promising technique for CMUTs fabrication because it is easy and reliable to control and reduces process cost and time.

Fig. 5. A circular CMUT cell and an array of CMUTs
The dynamic measurements of a metallized membrane are conducted in the time domain by means of a POLYTEC laser vibrometer.The vibrometer is composed of a fiber-coupled vibrometer sensor head OFV-534 and an OFV-5000 controller.The controller has a DD300 decoder being able to measure displacements up to 150 nm peak to peak at 24 MHz.In our experimental test the external force acting in the membrane is a step force and Fig. 6 shows the time response of the membrane center, where the sampling frequency of signal is 500 MHz and the duration is 20 s μ .Only a part of this time response of the structure is used in the identification procedure (the damped response).In this paper a single degree of freedom of a spring-mass-damper model is used to study the membrane behavior.The model is constituted by the following parameters: the metallized membrane mass m concentrated in its center, the damping coefficient c and the stiffness coefficient k.Table 1 shows the dimensions and material properties of the CMUT cell and the first theoretical resonance frequency of the metallized membrane can be computed by the relation given in Rossi [17]

. Dimensions and properties of the CMUT cell
The metallized membrane mass is m=  π R 2 tm = 6,61x10 -11 kg.The membrane stiffness is given by k=m (2 π f ) 2 , the damping coefficient is c=4 π m f ξ , the quality factor is Q=1/2 ξ and the half- power bandwidth of the system is B Δ = f/Q, where f is the resonance frequency of the membrane and ξ the damping factor.These two modal parameters are obtained by an average over the orders of the stabilization diagram presented in Fig. 7. To our knowledge, it is the first time that the subspace algorithm is applied to a microelectromechanical system and in particular to a CMUT cell.Fig. 8 shows the modulus of the CWT and the identified ridge obtained from the membrane center displacement.From the ridge of the CWT we obtain Fig. 9 and the instantaneous modal parameters of the membrane using the slope of straight lines.To our knowledge, it is the first time that the CWT algorithm is applied to the membrane of a CMUT cell.Tab. 3. Identified parameters of the metallized membrane using the continuous wavelet transform method Fig. 10 shows the comparison between the measured time response of the metallized membrane and the reconstructed time response obtained using the identified modal parameters obtained by the subspace method.In this time response we have only taken in consideration the damped response portion between 6,5 μs and 20 μs .Gm and β = c/(2m) can also be evaluated directly by interpolating the measured response of the membrane and in this case we can obtain the value of the damping coefficient c from β and m.In [20]   nonlinearities are added in the model to analysis the performances of the CMUT for several DC voltages.

Conclusion
In modal parameter identification the key point is to determine the relationship between the system parameters and the measured dynamical data and a particular attention has been paid to modal parameter identification when only output measurements are used.Such parameters can then be used to improve the ultrasonic generating process by a CMUT, for model validation and for fault detection in CMUTs.In this paper, we have proposed two methods for modal parameter identification: the subspace method based in the state-space innovation model and the continuous wavelet transform method based in the use of the ridge of the CWT.Numerical and experimental results have shown the effectiveness of the procedures presented in the paper for modal parameter identification.However, the CWT has some limitations: if the number of data points is insufficient or if the vibrating system is quickly attenuated closely modes are not well identified in terms of damping factors.

1 ]- 1 -
is called the innovation component in the data vector yk .It is an (mx1) stochastic vector with zero mean and covariance matrix Q = E[ek e'k].Note that ek is uncorrelated with y kand hence with zk by construction, and that the elements of {ek} are serially uncorrelated.zk is an (nx1) vector of unobservable states.The (nxm) matrix B is called the Kalman gain.Since E[xk+1 e k ] = E[xk+1 ek'](E[ek ek']) -1 ek (eq.15) . The singular values of  H are contained in descending order in the diagonal matrix  S .When estimating a model with n states, only the n largest singular values are used.Setting the smallest singular values to zero constitutes the set of exclusion restrictions used to identify the model.Denote the approximation to the block Hankel matrix created by enforcing these exclusion restrictions by n H ˆ .The S.V.D. of n the above the estimated singular values smaller than the cut-off value n s ˆ are all replaced by zero, and only the columns vectors of Û and  V corresponding to the included singular values are retained as The model specification problem has been made an approximation problem.The two factorizations of the block matrix n H ˆ are n

ω
is the oscillating frequency of the analysing wavelet.The Morlet wavelet is basically a sinusoidal function, oscillating at the frequency o ω and modulated by a Gaussian envelope of unit width.In the frequency domain the dilated Morlet wavelet becomes

Fig. 2 .
Fig. 2. Time response for the analyzed two-DOF system (a) and its spectrum (b) and (c)

Fig. 3 . 1 f ˆ = 2 MHz; 2 f 1 ξ 2 ξ
Fig. 3. Stabilization diagram on eigenfrequencies and damping factors for the two-DOF system Our procedure can separate closely spaced modes and the mean values on identified eigenfrequencies and damping factors obtained by an average over the orders of the stabilization diagram are 1 f ˆ = 2 MHz;

Fig. 6 .
Fig. 6.Time response of the membrane center to a step force actuation

Fig. 7 .
Fig. 7. Stabilization diagram on resonance frequency and damping factor for the experimental CMUT membrane

Fig. 10 . 7 .
Fig. 10.Comparison between the measured time response of the membrane (in blue) and the reconstituted time response (in red)

Table 2 and
Table3show the identified metallized membrane parameters using respectively the subspace method and the CWT method.These two methods give very similar results.