Linear elasticity
In order to construct a meaningful, yet computationally feasible computing scheme, material properties of the tissue should be properly modeled. In general, the biological tissue is a nonrigid object that deforms over time and has very complicated material properties in terms of the underlying constitutive laws [43]. However, in this paper, we only study the quasistatic ultrasound elastography, which is to utilize a steady state quasistatic excitation on the tissues using the ultrasound transducer. Since the tissue is compressed slightly during quasistatic elastography, constitute law of linear elasticity is still valid [2]. Hence, for computational simplicity, in our current two dimensional implementation, the linear isotropic continuum material is adopted to describe the mechanical behavior of the tissue during the quasistatic elastography, and the relationship of the stress σ and the strain ε obeys the Hooke’s law [44]:
\begin{array}{l}\sigma =S\epsilon \end{array}
(1)
where S is the strainstress matrix.
For the two dimensional case discussed in this paper, let u(x, y) and v(x, y) be the displacement along the x and yaxis of a point, the infinitesimal strain tensor ε of the point is:
\begin{array}{l}\epsilon =\left[\begin{array}{c}\frac{\mathrm{\partial u}}{\mathrm{\partial x}}\\ \frac{\mathrm{\partial v}}{\mathrm{\partial y}}\\ \frac{\mathrm{\partial u}}{\mathrm{\partial y}}+\frac{\mathrm{\partial v}}{\mathrm{\partial x}}\end{array}\right]\end{array}
(2)
In twodimensional situation, under the plane strain condition [44], strainstress matrix S can be derived as:
\begin{array}{l}S=\frac{E}{(1+\nu )(12\nu )}\left[\begin{array}{ccc}1\nu & \nu & 0\\ \nu & 1\nu & 0\\ 0& 0& \frac{12\nu}{2}\end{array}\right]\end{array}
(3)
S is a materialdependent matrix, E stands for the Young’s modulus, and ν stands for the Poisson’s ratio. Here, the Young’s modulus E and the Poisson’s ratio ν are two materialspecific parameters. This fact that the internal stress caused by the deformation is a function of the displacement vector and the material parameters is quite clear from these equations. Ideally, the problem could be tackled in three dimension to avoid the throughplane motion effect, but in this paper we only study the twodimensional problem because the images of quasistatic ultrasound elastography are two dimensional data collected by the linear array of ultrasonic probe. Further, a linear material model is used to illustrate the basic ideas and rationales of a joint estimation strategy to recover the distribution of real elasticity for quasistatic ultrasound elastography, but more realistic models can be easily adopted into current framework in the future works.
Stochastic finite element method
In the past decades, the deterministic finite element method has been able to provide an effective and convenient platform for biomechanical studies in the past decades [45, 46]. However, it does not have the capability to process the uncertainties of material properties, and kinematic observations, i.e., measurements of the displacement. Especially the ultrasound imagingderived data are usually corrupted by noises of various nature. It is thus necessary to develop a strategy, the stochastic finite element method (SFEM), which has been used widely for structural dynamic analysis in probability analysis frameworks [47, 48], to process the data from the ultrasound elastography. The main difference between SFEM and finite element method (FEM) is that material properties are deterministic variables in FEM, but are described by random fields, possibly with known prior statistics in SFEM.
In order to build our implementation of SFEM for the quasistatic ultrasound elastography, a Delaunay triangulated finite element mesh is constructed at the first image frame before compression. An isoparametric formulation defined in a natural coordinate system is used, where, for trinodal linear element, the basis functions are linear functions of the nodal coordinates [44]. With mesh constructed, assuming that the material parameters E (Young’s modulus) and ν (Poisson’s ratio) are temporally constant but varying spatially, we can have the following dynamic governing equation [44]:
\begin{array}{l}M\mathit{\xdc}+C\stackrel{\u0307}{U}+\mathit{\text{KU}}=R\end{array}
(4)
where M, C and K are the mass, damping and system stiffness matrices, R is the load vector, and U is the displacement vector. In this work, U consists of u(x, y) and v(x, y). The calculation of M, K and C matrices are the standard FEM procedure, which can be read in [44]. Because the tissue density can be generally considered to be uniform over the region of interest, M is a known function of the material density and is temporally and spatially constant. K is a function of the constitutive law of material, which is linear elasticity here, and is related to the materialspecific Young’s modulus and Poisson’s ratio, which are considered as constants temporarily in this study. In our framework, these two local material parameters are treated as random variables with known a priori statistics, and will be estimated along with the motion parameters. The damping matrix C is frequency dependent, and we assume small proportional Rayleigh damping with C = α M + β K in our implementation. In Rayleigh damping α is the mass proportional Rayleigh damping coefficient and β is the stiffness proportional Rayleigh damping coefficient [44]. In practice, it is difficult to determine the damping parameters because they are frequency dependent. We assume Raleigh damping for the very low damping biological tissue during quasistatic elastography, and fix the two weighting constants to be 1%.
State space representation
In order to employ our filtering strategy to integrate biomechanical model, the dynamics equation (Equation (4)) needs to be transformed into a statespace representation of the continuoustime linear stochastic system. First the kinematic vector x(k) and the material parameter vector θ are defined as:
\begin{array}{l}x\left(t\right)=\left[\begin{array}{c}U\left(t\right)\\ \stackrel{\u0307}{U}\left(t\right)\phantom{\rule{1em}{0ex}}\end{array}\right]\end{array}
(5)
\begin{array}{l}\phantom{\rule{9.0pt}{0ex}}\theta =\phantom{\rule{1em}{0ex}}\left[\begin{array}{c}E\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\\ \nu \phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\end{array}\phantom{\rule{1em}{0ex}}\right]\end{array}
(6)
Where the kinematic vector x(t) is consisted of displacement U(t) and velocity \stackrel{\u0307}{U}\left(t\right) information, and the material parameter vector θ is consisted of Young’s modulus E and Poisson’s ratio ν. In general, tumor inclusion or tissue blocked from its blood nutrients is stiffer than normal tissue, which mostly reflects in the variation of E. Since benign and cancerous tumors usually have distinguishing elastic properties (i.e., different Young’s modulus value), the value of Poisson’s ratio can be fixed in the following implementation. For example, as one of incompressible materials, the Poisson’s ratio of tissue can be set close to 0.5.
The state space representation of Equation (4) becomes
\begin{array}{l}\stackrel{\u0307}{x}\left(t\right)={A}_{c}x\left(t\right)+{B}_{c}W\left(t\right)+n\left(t\right)\end{array}
(7)
where n(t) is the process noise which is an additive, zeromean, white noise (E[n(t)] = 0; E[n(t)n(s)^{′}] = Q_{
n
}(t)δ_{
ts
}, where Q_{
n
} is the process noise covariance). Process noise cannot be ignored if modeldata discrepancy exists. The system matrices A_{
c
}, B_{
c
} and input forces W(t) are:
\begin{array}{lll}\hfill {A}_{c}& =& \left[\begin{array}{cc}0& I\\ {M}^{1}K& {M}^{1}C\end{array}\right]\\ \hfill {B}_{c}& =& \left[\begin{array}{cc}0& 0\\ 0& {M}^{1}\end{array}\right]\\ \hfill W\left(t\right)& =& \left[\begin{array}{c}0\\ R\left(t\right)\end{array}\right]\end{array}
An associated measurement equation, which describes the observations y(t) extracted from the ultrasonic RF data, can be expressed in the form:
\begin{array}{l}y\left(t\right)=\mathit{\text{Hx}}\left(t\right)+e\left(t\right)\end{array}
(8)
where e(t) is the measurement noise which is additive, zero mean, and white (E[e(t)] = 0;E[e(t)e(s)^{′}] = R_{
e
}(t)δ_{
ts
}, where R_{
e
} is the measurement noise covariance), independent of n(t). Although the noise in the displacement measurement extracted from ultrasonic RF data might not be white Gaussian noise, the property of the noise from ultrasonic RF data is rarely available due to the complicated collecting process of ultrasonic RF data. Hence our best guess is to assume it is white Gaussian noise, which is one common treatment in engineering field [49–51]. Furthermore, we found that the performance of this assumption is satisfied in the following experiments. H is the measurement matrix which should be specified by the relation between state vector x(t) and measurement vector y(t).
In order to run the dynamic equation (Equation (4)) in computer, it should be transformed into a discrete statespace equation, typically seen in control and estimation literature [50, 52]. We discretize Equation (7) and (8) over a constant time interval T. Since the interval T is always a known constant, we can replace kT with k in following equations
\begin{array}{l}x(k+1)=\mathit{\text{Ax}}\left(k\right)+\mathit{\text{Bw}}\left(k\right)+n\left(k\right),\end{array}
(9)
where
\begin{array}{lcr}A& =& {e}^{{A}_{c}T},\end{array}
(10)
\begin{array}{lcr}B& =& {A}_{c}^{1}({e}^{{A}_{c}T}I){B}_{c}\end{array}
(11)
The associated discrete measurement equation is:
\begin{array}{l}y\left(k\right)=\mathit{\text{Dx}}\left(k\right)+e\left(k\right)\end{array}
(12)
where y(k) is the measurement vector contained the displacement extracted from ultrasound RF data. In most quasistatic elastography, only the axial components of displacement vector are extracted. Therefore, the corresponding place of D will be set to 1 or 0 according to the available measurement data. n(k) and e(k) are process noise and the measurement noise respectively.
In order to perform the estimation of the material parameter distribution, the state vector x is augmented by the material parameter vector θ to form the new state vector z:
\begin{array}{l}z={[x,\theta ]}^{T}\end{array}
(13)
Accordingly, the new augmented state equation is:
\begin{array}{l}z(k+1)=f\left(z\right(k),w(k\left)\right)+{v}_{s}\left(k\right)\end{array}
(14)
with
\begin{array}{l}f\left(z\right(k),w(k\left)\right)=\left[\begin{array}{c}A\left(\theta \right)x\left(k\right)+B\left(\theta \right)w\left(k\right))\\ \theta \end{array}\right]\end{array}
(15)
\begin{array}{l}{v}_{s}\left(k\right)=\left[\begin{array}{c}v\left(k\right)\\ 0\end{array}\right]\end{array}
(16)
Where v(k) is the process noise, introduced by the model mismatch. The new augmented measurement equation is derived:
\begin{array}{l}y\left(k\right)=h\left(z\right(k\left)\right)+{e}_{o}\left(k\right),\end{array}
(17)
with
\begin{array}{l}h\left(z\right(k\left)\right)=[D,0]\left[\begin{array}{c}x\left(k\right)\\ \theta \left(k\right)\end{array}\right]\end{array}
(18)
\begin{array}{l}{e}_{o}\left(k\right)=e\left(k\right)\end{array}
(19)
Using previous assumptions of noises, the augmented noises are also Gaussian with distribution:
\begin{array}{l}{v}_{s}\left(k\right)\sim N(0,{Q}_{s}),\phantom{\rule{1em}{0ex}}\mathit{\text{where}}\phantom{\rule{1em}{0ex}}{Q}_{s}=\left[\begin{array}{cc}{Q}_{v}& 0\\ 0& 0\end{array}\right]\end{array}
(20)
\begin{array}{l}{e}_{o}\left(k\right)\sim N(0,{R}_{o}),\phantom{\rule{1em}{0ex}}\mathit{\text{where}}\phantom{\rule{1em}{0ex}}{R}_{o}={R}_{e}\end{array}
(21)
Extended kalman filter (EKF) for joint state estimation
The reconstruction of quasistatic elastographic image now can be defined as calculating the optimal estimates of the nonlinear statespace estimation problem represented by the Equation (14) and (17) above. This leads to solve this nonlinear estimation problem using the EKF framework, which is to linearize the augmented state (Equation (14)) at each time step and perform a recursive procedure with natural block structure to estimate the material parameter. EKF operates like Kalman filter, adopting the same form of feedback control in estimation: the filter makes an estimation of the state at one moment and then obtains feedback in form of measurements. The operations of EKF can be divided into two steps: time update equations and measurement update equations. Time update equations are responsible for projecting forward the current state and error covariance estimates to obtain the priori estimates for the next time step, while the measurement update equations are responsible for the feedback with the posterior data.
Initializing the filter with \hat{z}{}^{}\left(0\right)={\hat{z}}_{0} and P\left(0\right)=\sum _{0}, the EKF runs sequentially as follows:

(1)
compute the prior estimation of the state, from t−1 to t:
\begin{array}{l}\hat{z}{}^{}\left(t\right)=f\left(\hat{z}\right(t1),w(t\left)\right),\end{array}
(22)

(2)
project the error covariance from t−1 to t:
\begin{array}{l}P{}^{}\left(t\right)={F}_{t}P(t1){F}_{t}^{T},\end{array}
(23)

(3)
compute the Kalman gain at t:
\begin{array}{l}G\left(t\right)={P}^{}\left(t\right){H}^{T}{\left(H{P}^{}\right(t){H}^{T}+{R}_{0})}^{},\end{array}
(24)

(4)
compute posterior estimation of the state with the measurement at t:
\begin{array}{l}\hat{z}\left(t\right)=\hat{z}{}^{}\left(t\right)+G\left(z\right(t)h(\hat{z}{}^{}\left(t\right)\left)\right),\end{array}
(25)

(5)
update the error covariance at t:
\begin{array}{l}P\left(t\right)=(IG(t\left)\right)\left(z\right(t)h(\hat{z}{}^{}\left(t\right)\left)\right),\end{array}
(26)

(6)
repeat (1)(5) to convergence. The needed quantities for Equations (22)–(26) are defined by:
\begin{array}{l}{F}_{t}=\frac{\partial}{\partial}f\left(z\right(t),w(t\left)\right){}_{z=\hat{z}}=\left[\begin{array}{cc}A\left(\hat{\theta}\right)& {M}_{t}\\ 0& I\end{array}\right]\end{array}
(27)
\begin{array}{l}H\left(\phantom{\rule{0.3em}{0ex}}\hat{z}\right(t\left)\right)=\frac{\partial}{\partial}h\left(z\right(t\left)\right){}_{z=\hat{z}}=\left[\begin{array}{cc}D& 0\end{array}\right]\end{array}
(28)
\begin{array}{l}{M}_{t}=\frac{\partial}{\partial \theta}\left(A\right(\theta \left)\hat{x}\right(t)+B(\theta \left)w\right(t\left)\right){}_{\theta =\hat{\theta}}\end{array}
(29)
In our implementation, the initial error covariance matrix P(0) consists of four parts:
\begin{array}{l}P\left(0\right)=\left[\begin{array}{cc}{P}_{1}\left(0\right)& {P}_{2}\left(0\right)\\ {P}_{2}^{T}\left(0\right)& {P}_{3}\left(0\right)\end{array}\right]\end{array}
(30)
P_{1}(0), P_{3}(0) are the submatrices inside the initial error covariance matrix related to the trustworthiness of the kinematic functions and material parameters respectively, and P_{2}(0) is the kinematicsmaterial correlation submatrix, which is a zero matrix in this paper.
Summary of the algorithmic flow
As described above, the system vector z, which consists of kinematic vector x and material vector θ, is updated through the EKF filtering framework, but in this implementation we only have one set of data from quasistatic elastography, i.e., one frame of displacement measurement acquired in deformed configuration. In Figure 1, the flow chart of our filtering framework is illustrated. If the difference between the estimate of current time step and the estimate of previous time step is small enough, the operation of EKF will be terminated and the final estimate will be considered as close enough to the optimal estimation. Therefore, through our filtering framework, a quantification of the internal material properties distribution jointly with kinematic state could be recovered from the ultrasonic RF data collected from quasistatic elastography.