Project topics and materials
Get Complete Project Material File(s) Now! »
State of the art in inverse problems. The problem of deconvolution
In physics, an inverse problem appears when determining the input or the inter-nal characteristics of a physical system from external measurements. This kind of problems appears in many physical systems where the output depends on the prop-erties of the system and an input. Trying to infer the external applied force from measured displacements, or the structure of the soil from known electromagnetic inputs and outputs are two examples of inverse problems. The main characteristic of inverse problems is that the small errors introduced in the measurements can affect considerably the solution. In most of the problems, the analytic solution is not available, and the direct resolution of the problem can lead to a solution that is far from the exact one. By definition, inverse problems suffer of lack of infor-mation making impossible the reconstruction of the input, and numerical methods must be considered to find an approximation of the solution. Even if perfect mea-sures (i.e. noiseless) could be done, one would still need to perform computations in infinite precision to avoid noise amplification. Therefore, approaches based on each problem characteristics must be considered. From the mathematical analysis to the numerical solution of the problem, numerous authors have contributed to the analysis of inverse problems.
The mathematical standard form for inverse problems in many physical sys-tems can be formulated as the Fredholm integral equation of the first kind, which is written as [44]: u(x) = h(x; y)f(y)dy; (1.1)
where h(x; y) is called the kernel, and f(y) and u(x) are the input and the output functions respectively. The function f(y) is to be determined, h(x; y) is known in the interval a x; y b and u(x) is known in the interval a x b. Consider a physical system, e.g. a tomograph, which can be modelled and considered the kernel. The obtained images will be then the outputs u(x), which will be probably obtained with a certain amount of errors due to the precision of the machine, de-fects, noise, etc. One wants to find the input f(y), i.e. the function that represents the body organ, with the best possible accuracy.
A general idea of the inverse problems can be obtained from the consideration of the mapping between the input space Si and the output space So: small perturba-tions in the object space lead to small perturbations in the image space. However, in the inverse mapping small perturbations in the image space can lead to big per-turbations in the object space [49], leading to solutions far from expected. This feature is shown graphically in figure 1.1. Therefore, computation of f(y) could lead to an erroneous representation of the body organ.
A classification can be made in the inverse problems field [35]. When the known variables are the input and the output the problem is classified as a system identification problem. When the known variables are the kernel and the output, the problem is classified as a deconvolution problem (the word unfolding can be also found in some textbooks). The previous definitions are summarized in table 1.1, extracted from [35].
The use of system identification methods started to be important since 1960s. The augment of data collection and measurement techniques allowed new possi-bilities for physicians. New techniques were developed taking profit of the data to obtain models for unknown systems. These models considers the system as black boxes that relates inputs and outputs. The identification of this black box is the objective of the system identification techniques [57]. [86, 87] are two examples where system identification techniques are used to estimate unknown parameters of the models. The exact composition and materials of the earth layers is not usually accessible, or there not exist a model on a complex physic process. System identi-fication techniques use the collected data, inputs and outputs, to try to identify the hidden physical system or the unknown parameters of the model.
Deconvolution problems appear in physical systems where the kernel is known and the problem can be modelled by the convolution theorem. Wave or heat equa-tions allow to model with accuracy dynamics or heat transfer in structures, and are usually applied to compute displacements and temperatures.
Both system identification and deconvolution share numerical methods for res-olution procedures. When the kernel is a convolution operator, the properties of the integral allow to interchange the kernel and the input roles, which makes no difference between both fields. In this work, despite this equivalence, the inverse problem will be referred as a deconvolution problem. The reason can be found on the structure of this thesis: a direct problem will be first solved by using the convolution operator, so for the sake of simplicity, is called deconvolution to the computation of the inverse problem which is treated subsequently, even if it is the kernel what will be computed.
The convolution operation appears in many physical systems. It is a special case of the Fredholm equation of the first kind where the kernel is a convolution operator. The integral in this case can be written as: u(t) = h(t )f( )d ; (1.2) where the time parameter t 2 It is introduced and a range of [0; T ] is considered. For the range [ ; 0] the input is zero, so it is the output, and causality holds. This equation is also known as Volterra equation of the first kind.
One of the most known problems in the community, used as an example of deconvolution problems, is the geomagnetic prospecting problem [44]. Consider that a deposit is situated under the surface of the earth in a plane stratum at a known distance h from the surface, as shown in figure 1.2.
This deposit emits a magnetic signal whose intensity determines its volume, then its economic value. Measures of the vertical component J(x) of the magnetic field intensity at the surface S(x; 0) are available. The field is created by a mass whose distribution I(x) is required. The value of J(x) in the surface S due to an infinitesimal part ds of the mass is:
sin I(s)ds = hI(s)ds : (1.3)
(p (h2 + (x s)2) 3
h2 + (x s)2)2
The deposit occupies a space from x1 to x2 , then the total contribution is:
Js(x) = x2 (1.4)
Zx1 (h2 + (x s)2)23 :
I(s)ds
Obtaining I(s) from (1.4) is a deconvolution problem. Many other physical systems require the resolution of the deconvolution problem, among them one can find astronomy [73], biology [71] or tomography [74].
Deconvolution is a non well-posed problem in the sense of Hadamard definition
[41]. Hadamard defined the term well-posed problem when it fulfils the following conditions:
The solution exists.
The solution is unique.
The solution depends continuously on the data.
When a problem does not fulfill at least one of the conditions, the problem is called ill-posed problem. Inverse problems are generally ill-posed problems, and consequently, the deconvolution problem is generally ill-posed. Furthermore, it is usually ill-conditioned, which means that small variations in the data produce big variations in the solution, so special theoretical and numerical approaches must be applied to solve the problem. The convolution problems are usually described as u~(t) = h(t )f( )d ; (1.5) where u~(t) = u(t) + e. This means that an error has been added to the output, and it is generally unknown. These errors come principally from measurement errors and truncation (identified as noise in signal treatment). The introduction of these errors can cause unacceptable deviations in the solution, even if they are negligible from a physical point of view.
The theory about the Fredholm equations of the first kind states that the inte-gration of f with h is a smoothing operation, what means that high frequencies are damped. In fact, the convolution of f and h maximizes the differentiability of the functions f and h [36], in other words, the function u is smoother than f. By following this property, the deconvolution operation will undo the smoothing trans-formation. If in the convolution operation the high frequencies are damped, in the deconvolution the high frequencies are expected to be augmented. Thus, the rela-tion between the errors and the high frequencies is a key point in the deconvolution procedure.
Discrete convolution and deconvolution. Properties of the Toeplitz matrix
The integral of convolution (1.2) can be written in discrete form as:
Xi (1.6)
u(k) = h(k i)f(i) =0
The equivalent version in matrix form is: u = Hf; (1.7) where H is the Toeplitz matrix of the impulse response h, where the entries of the matrix have been added until a time step n, and the corresponding values in the range of time [ ; 0] are set to zero. The Toeplitz matrix takes the form:
2 h(2) h(1) ::: 0 0 3
6 h(1) 0 ::: 0 0 7
H = … … … … … : (1.8)
This Toeplitz matrix is a n n matrix, where only 2n 1 elements are different and there are placed with a characteristic form: each element of h occupies one diagonal of the matrix. This feature makes Toeplitz matrices ill-conditioned, but the repetitive structure makes possible to derive appropriate methods to im-prove the resolution of the problems involving Toeplitz matrices. One of the main characteristics is the persymmetry, which means that they are symmetric across the anti-diagonal.
The previous Toeplitz matrix is a special case of the more general Toeplitz matrix definition [63]: 3
2 h(2) h(1) ::: h( n 2) h( n 1)
6 h(1) h( 1) ::: h( n 1) h( n) 7
H = … … … … … : (1.9)
A Toeplitz matrix is then symmetric when h(j) = h( j) for all j.
The direct resolution f = H 1u by inverting the Toeplitz matrix H is usually precluded because it is usually ill-conditioned [43]. The relative error of the so-lution u is related to the condition number of H. How large a condition number is acceptable depends on the resolution accuracy or the tolerance associated to the problem.
Let us consider the discrete form of the deconvolution problem to obtain f, which can be easily deduced from equation (1.6):
u(k) P k 1 u(i + 1)h(k i + 1)
f(k) = i=1 : (1.10)
h(1) h(1)
It can be observed the importance of the errors in measured u. Errors in the in-put u are inherent to discrete calculus, and they are reintroduced in the subsequent computations, leading to an unstable resolution. The equation (1.10) shows that the importance of the errors in the solution grows proportionally to the factorial of the number of computations k.
Numerical solution methods
The existence of inherent errors in the inverse problems can be considered as a lack of information, and this lack of information can lead to a lack of uniqueness of the solution. The lack of information can never be avoided, so the approaches are based on the addition of information to the problem: statistical knowledge of the noise, cost function energy control, etc.
Different techniques exist in the literature to solve the lack of information, which try in some way to add information or constraint in the resolution process to reduce the subspace of possible solutions. Ivanov, Tikhonov and Phillips [50] pro-posed a minimization problem, which is actually known as regularization method. Iterative variants of the regularization methods, where the solution is sought in an iterative procedure until a certain prescribed error is achieved, also exist [47]. The Truncated Singular Value Decomposition can also be used to deal with the ill-conditioning of the Toeplitz matrices [63]. Another kind of approaches comes from the statistical analysis of the problem by introducing the theory of probabili-ties to cover the lack of information [48]. There has been also attempts to solve the deconvolution problem in real time, where some methods can be found in [54][55].
The resolution of the problem in the frequency domain is based in the Fourier’s theory [8]. With the development of the Fast Fourier Transform, the computational cost of the transformation between the time domain and the frequency domain was considerably reduced. Consequently, a wide family of methods and applications has emerged in the last decades taking profit of its benefits. Signal processing, seismography, tomography, etc. are some of the fields where it is applied [21].
The resolution of the problem in the frequency domain eliminates the bad con-ditioning which arises when inverting the Toeplitz matrix form of H. The equa-tions are transformed to the frequency domain, where the convolution becomes a multiplication [49]:
u^ ^^ (1.11) = hf;
where the symbol (^) denotes that the function lives in the frequency domain, so u^ ^ ^ 2 C. ^ ; h and f To obtain f, a simple division is made, which means that the deconvolution in the frequency domain is transformed into a division between the functions u^ and h as follows:
^ u^ (1.12)
f = ^ : h
The existence of a solution depends on the r.h.s. of the equation: it must be well-defined, which is not always the case. Consider the case when h is zero or near zero for one or more values of !: then the solution in frequency will have singularities and the inverse Fourier transform may not exist.
Fourier analysis gives a useful view about the deconvolution problem. Remem-ber that the acquired data is usually obtained in discrete form, and the analytic solution is not available. Errors in the data is unavoidable as the precision of the machines is not infinite. If no errors were added to the output data, the deconvo-lution would give back the exact searched function f, but usually is not the case. Consider then the equation (1.12) and add the noise errors u^ in the acquired data:
^ u^ + u^ u^ u^ (1.13)
f = = + :
^ ^ ^
h h h
Physically motivated transfer functions tend usually to zero when the frequency tends to infinity. On the other hand, transfer functions related to noise usually show a nearly flat spectrum (e.g. white noise has equal intensity at all frequencies). The division of the high frequency zone of the data spectrum by the high frequency zone of the transfer function spectrum will amplify the effect of the noise, and will dramatically affect the solution.
Truncated SVD
The singular value decomposition is a factorization technique for matrices. Over the field of real numbers, the SVD decomposition is defined as: H = UDVT; (1.14) where U is the matrix of the left-singular vectors, D is a pseudo-diagonal matrix containing the singular values and V is a matrix that contains the right-singular vectors. The singular vectors are stored by columns in U and V, so their columns are orthonormal: UTU = VTV = I; (1.15) where I is the identity matrix. The decomposition can be also written as:
Xi (1.16)
H =i i i; =1
where n is the rank of H.
Then, the right and left singular vectors are stored in matrices U and V: U = ( 1; : : : ; n); V = ( 1; : : : ; n); (1.17) and the right and singular vectors are orthogonal following: iT j = iT j = ij: (1.18)
The singular values reveal some important properties of H:
The singular values of H decay with a certain slope until it reaches the ma-chine precision times 1, where the slope tends to be horizontal.
The previous statement means that the condition number is related to the machine precision.
There is no gap in the singular values spectrum, and it usually follows a harmonic or a geometric progression.
The number of sign changes of the singular vectors is inversely proportional to the value of the corresponding singular value.
The number of non-zero singular values is the rank of H, and the condition number of H can be computed from: cond(H) = 1 ; (1.19) where 1 is the largest singular value and n is the minimum non-zero singular value.
The SVD permits to stablish an important relation between the matrix H and the singular values and singular vectors: H i = i i; (1.20) where the L2 norm of H i is: jjH ijj = i: (1.21)
These properties can be exploited to analyse the problem of deconvolution.
Consider then the problem of the equation (1.7): Hf = u: (1.22)
By using (1.20), the terms u and f can be expanded as: f = ( iT f) i ; u = ( iT u) i; (1.23) i=1 =1 and substituting (1.23) in (1.20), is straight to write the equation: Xi ( iT u) i: (1.24) Hf = i =1
Finally, the approximated solution can be written as a decomposition in singu-lar values and singular right and left vectors as: f = i i i: (1.25) i=1
This decomposition allows to compute an approximation of the solution by truncating the sum. The truncation can be made by controlling the singular values and the corresponding left and right eigenvectors. In fact, the solution is governed by the relation between Ti u and the singular values. The Discrete Picard Con-dition [64] states that the quantity j Ti uj decays faster than the singular values i until they arrive to the level of the machine precision. In this point, the round noise starts to govern over the SVD components, and the information obtained in the decomposition is no longer useful. It is evident that the terms Ti u and the sin-gular values must be monitored to avoid the introduction of spurious terms on the approximation of the solution. If the quantity Ti u does not decay faster than the singular values, there is not any function square integrable f that is a solution of the problem such that: jf(t)j2dt < 1; (1.26) and any effort solving the problem should be avoided.
Therefore, useful information can be extracted from the SVD decomposition of matrix H. The first singular values and singular vectors still carry the most important information: the effort must be then in the separation of the valuable information and the spurious information, by applying a truncation:
~ N iT u
f = Xi i; (1.27)
=1 i
where N < n. This approximation is known as Truncated Singular Values De-composition. The election of the truncation number N is quite simple: a SVD decomposition must be done on matrix H, and the quantities Ti u must be com-puted. The number N corresponds to the change of the slope of the quantities Ti u. There exists methods that reduce the computational cost of the previous method [65].
The main problem in the TSVD method is the computational cost: the decom-position of large matrices can be unaffordable for the actual computational power.
Regularization method
At the same time, Ivanov and Tikhonov published similar techniques to solve the Fredholm equation of the first kind [44][49], and their theories can be also applied to the deconvolution problem. These methods are known as Tikhonov regulariza-tion method or just regularization method. Those techniques use a least squares approximation with a prescribed energy weighted by a parameter to resolve the problem. As the solution is not unique, a least squares technique is applied to find a unique solution, which will be the best approximation to the solution in the con-strained problem. In the discrete form of the regularization method, f is sought in order to minimize the functional :
(f) = jjHf ujj2; (1.28)
with a prescribed energy for f:
E2(f) = jjfjj2 E2: (1.29)
The bound E2 of f must be in the order of the known u. A bound for u is not usually a simple guess. Consider a structural dynamic problem, where the structure characteristics are known. In most of the engineering problems, an order of magnitude of the displacements is usually known, so an order of magnitude of the displacements is also guessed. This added information is what can in some way substitute the lack of information.
This problem can be solved by applying the method of Lagrange multipliers, where the condition is introduced in the same functional of the problem as: = jjHf ujj22 + 2jjfjj22: (1.30) where is an arbitrary positive real number. The main problem of the method is the estimation of the parameter . Great values of this parameter will lead to a solution which may differ significantly from the real solution, while very small values will lead to the ill-conditioned problem, where the solution exists but can be far from what expected. Several methods exist on the research of the optimal value of [61][62][60]. In the middle of the extreme values, there exists infinite possible solution for the problem.
The election of the regularization parameter is still a challenge for the research [67]. The different possible solutions are related to the information one has about the problem. The distribution of the error, or the assumptions that can be made about the correlation with the solution leads to different methods. The discrepancy principle searches the regularization parameter as a function of the error norm, which must be known or assumed known. The cross validation method consid-ers that the noise is uncorrelated with the solution, and therefore some previous information is needed. One of the most known methods is the L-curve criterion [34][66], which consists in computing the regularized solution for several different values of the parameter > 0. The plot of the f versus the norm of the residual is a curve in a more or less L shape, where a minimum can be found along the curve. [34][66] recommend using the corresponding value to the minimum of the plot, which is a good choice in the balance between the solutions influenced by the ill-conditioning of the problem and the solutions influenced by the noise.
Table of contents :
Introduction
1 State of the art in inverse problems. The problem of deconvolution
1.1 Discrete convolution and deconvolution. Properties of the Toeplitz matrix
1.2 Numerical solution methods
1.3 Numerical methods for Toeplitz matrices
2 Monitoring of displacements
2.1 Structural dynamics
2.1.1 Model problem
2.1.2 Time integration methods
2.1.3 Modal method
2.1.4 Harmonic Analysis
2.1.5 Example
2.2 Model order reduction techniques
2.2.1 Proper orthogonal decomposition
2.2.2 Reduced Basis Method
2.2.3 Proper Generalized Decomposition
2.3 Generalized transfer function
2.4 Generalized impulse response
2.5 Generalized displacements
2.5.1 Generalized displacements in frequency space
2.5.2 Generalized displacements in time domain
2.6 Results
2.7 Fractional damping
2.7.1 Fractional damping as a parameter
2.7.2 Results
3 Monitoring of forces
3.1 Formulation of the problem
3.2 Generalized inverse impulse response
3.2.1 Dual problem. Flexibility method
3.2.2 Training
3.2.3 Avoiding regularization with the separated representation
3.2.4 Computation of the generalized inverse impulse response
3.3 Results
4 Nonlinear applications of the Generalized Impulse Response
4.1 Generalized impulse response application in nonlinear problems
4.2 Nonlinear external applied force
4.3 Nonlinear stiffness
4.4 Numerical examples
4.4.1 Nonlinear external applied force
4.4.2 Nonlinear stiffness
Conclusions
Bibliography
Appendix
