Nonlinear inversion of ultrasonic guided waves for in vivo evaluation of cortical bone properties

Xiaojun Song(宋小軍)， Tiandi Fan(樊天地)， Jundong Zeng(曾俊冬)， Qin-Zhen Shi(石勤振)， Qiong Huang(黃瓊)，Meilin Gu(顧美琳)， Petro Moilanen， Yi-Fang Li(李義方)， and Dean Ta(他得安)，4，?

1Center for Biomedical Engineering，F(xiàn)udan University，Shanghai 200433，China

2College of Electronics and Information engineering，Shanghai University of Electric Power，Shanghai 200090，China

3Department of Physics，University of Helsinki，Helsinki，F(xiàn)inland

4Academy for Engineering and Technology，F(xiàn)udan University，Shanghai 200433，China

Keywords: ultrasonic guided waves，inversion algorithm，transverse velocity，bone evaluation

1. Introduction

Osteoporosis can lead to an increased risk of fracture due to the low bone mass and micro-architectural bone deterioration，[1，2]and can affect almost all skeletal sites.[3]Thus， it is important and necessary to perform a bone examination when someone is getting shorter or their upper back begins to curve forward. Nowadays， bone mineral density(BMD)is a widely accepted parameter for assessing bone status，and is tested by dual-energy x-ray absorptiometry(DXA).However，DXA has limitations of ionizing radiation and high cost. Therefore，it is difficult to test BMD frequently and regularly using DXA. In contrast， because of its advantages of being radiation-free， cheap and sensitive to mechanical properties， ultrasound was proposed for the evaluation of bone status.[4–6]

Cortical bone represents 80% of the weight of a human skeleton. The properties of cortical bone， such as cortical thickness and bulk velocities， are important parameters relating to bone status.[7–9]The axial transmission(AT)technique was first used to evaluate cortical bone 30 years ago.[10–13]The ultrasonic signals collected from cortical bone and recorded by the AT technique are mainly of two types. One type is first arriving signals (FASs) and the other type is ultrasonic guided waves (UGWs). FASs can correspond to lateral waves which propagate along the surface of cortical bone，while UGWs， which propagate through the whole thickness of the cortical bone， are more sensitive to bone minerals and architecture.[14–16]Thus，UGWs are very attractive for assessment of cortical bone.[17–22]Nicholsonet al.extracted the A0 Lamb mode， which is sensitive to the thickness of cortical bone.[10，23]Moilanenet al.obtained the cortical thickness in anin vitrostudy from the group velocity of UGWs based on an inversion scheme.[24]Songet al.proposed a blind identification method to extract a single UGW mode， and a mean error of 4.3%for cortical thickness was achieved.[25]Baiet al.used guided waves to assess the fatigue level of cortical bones，and the results showed that the mean phase velocities of Lamb modes A1 and S1 for intact bone were higher than those of fatigue-damaged bones.[26]

Only one parameter was obtained to evaluate the status of the long bones in the above-mentioned studies. However，cortical bone is affected not just by one parameter but by multiple parameters together，for example cortical thickness，bulk velocities， porosity and so on. Therefore， many researchers have begun to focus their attention on extracting multiple parameters for long bones. Pereiraet al.used low frequencies(20 kHz–85 kHz) to investigate intracortical bone and found that the phase velocity and cut-off frequency of UGWs were sensitive to the properties of cortical bone.[27]Minonzioet al.estimated cortical thickness and porosity level and found that these two parameters are sensitive to non-traumatic fractures in post-menopausal women.[28]Tranet al.developed an inversion procedure to estimate the cortical thickness and bulk velocities of one coated bone plate; their results showed that these parameters were accurately obtained.[8]Bochudet al.used two models (a free plate model and a bilayer model) of long bones to obtain cortical thickness and elastic velocities in bothin vitroandin vivoexperiments，and their results suggested that the free plate model was sufficient to inverse the parameters even in thein vivoexperiments.[7]Liet al.investigated the possibility of utilizing a neural network to obtain cortical thickness， transverse velocity and longitudinal velocity for bovine plates;their results demonstrated that the neural network method was able to assess the parameters of bovine cortical bones precisely.[21]

Following these studies， most researchers have limited their evaluation toin vitroorex vivobone-mimicking plates due to the soft tissue layer and irregular bone structure;few of them have identified osteoporosis using UGWs quantitatively.Considering that cortical thickness， longitudinal velocity and transverse velocity are three important parameters for evaluating bone status，we aim to develop a multi-parameter inversion procedure to acquire these three parameters simultaneously.Simulation andin vitroexperiments are used to test the accuracy of the proposed multi-parameter inversion procedure.In vivomeasurements are also carried out to find the quantitative relationship between osteoporosis and UGWs.

This paper is arranged as follows. Theory and methods are reported in Section 2， with details on forward calculation of the dispersion curves of UGWs， experimental dispersion curve extraction and the inversion algorithm. This is followed in Section 3 by the simulations，in vitroexperiments on bone plates andin vivomeasurements. Then the corresponding results and discussion are given in the Sections 4 and 5，respectively. Finally，our conclusions are given in Section 6.

2. Theory and methods

2.1. Computing theoretical dispersion curves

Two bone models are utilized in this research to calculate the theoretical dispersion curves. One is a free solid plate model (cortical bone)， and the corresponding guided modes are Lamb waves.[29，30]The other is a fluid–solid bilayer model(cortical bone coated by soft tissue)，and the corresponding guided modes are leaky waves because the energy of the guided waves leaks from solid to fluid.[8]The solid layer，which models cortical bone，is elastic and isotropic. The soft tissue layer is modeled as an acoustic fluid medium. The theoretical dispersion curves can be depicted by solving an eigenvalue problem

whereK1，K2， andK3are matrices related to global mass and rigidity，kis the wavenumber along the axial direction，i is the imaginary unit，ωis the angular frequency，Pα(α=1，2)is the global pressure andUis the displacement field.[8，31]By solving Eq.(1)，the eigenvaluekand the corresponding eigenvectorV(ω，k) of the UGWs can be obtained for a givenω.Then the theoretical dispersion curves can be calculated by sweepingωin a given range.

2.2. Extraction of experimental dispersion trajectories

In order to inverse cortical thickness，longitudinal velocity and transverse velocity from the UGWs， the dispersion curves for the simulations andin vitroandin vivoexperiments need to be extracted based on the ultrasonic signalsg(x，t). In the first step，g(x，t)collected in the time(t)and spatial(x)domains，was transformed to frequency(f)and wavenumber(k)domains based on the two-dimensional Fourier transform(2DFT)method.The distributions of mode energyH(ω，k)are the corresponding experimental dispersion trajectories，which can be expressed by the following equation:[32]

The 2D-FT method is easily affected by noise and frequency aliasing，and the finite number of receiving signals further limits the wavenumber resolution in long cortical bone.[21]Thus，in a second step， the Burg method was used to obtain a relatively high resolution with the finite sampling points in our case(details of the Burg method are shown in Ref.[33]).

2.3. Multi-parameter inversion procedure

After 2D-FT and the Burg algorithm， the experimental dispersion energy is obtained and saved in the form of a twodimensional matrixA，which is expressed as

wheremandnare the index ranges of frequency and wavenumber andai jis the corresponding energy value at point(i，j). The theoretical dispersion curves in the wavenumber and frequency domains can be written as a matrixB，

wherebi j=1 if the theoretical dispersion curves go through the(i，j)point;otherwisebi j=0.

The non-linear objective functionJis given by

whereTSis the threshold of the experimental dispersion energy. Ifaij ＞TS， then (ai j ＞TS) in Eq. (5) equals 1， otherwise it is 0.wei，jandwti，jare the wavenumbers of experimental dispersion data and the theoretical dispersion curve at (i，j)，respectively， andLis the point number that meets the threshold in the entire matrixA. The minimumJmininverses the optimal parameters (cortical thickness CTh， longitudinal velocityVLand transverse velocityVT)by the best-fit theoretical dispersion curve.

2.4. Finite-difference time-domain method

The simulated ultrasonic signals are calculated by the FDTD method.[34]For stability and precision， the space step Δx(spatial resolution) and time step Δt(temporal resolution)should satisfy the following expressions:

whereλmindenotes the shortest wavelength，dis the dimension of the model andCmaxis the fastest longitudinal bulk velocity in the simulation.

The inversion method was tested on two sets of simulated data. One was for a free plate and the other was for a bilayer model. Perfectly matched layers (PMLs) were on both sides of the models to absorb the reflection. The PMLs comprised an artificial absorbing medium and were able to achieve quasi-perfect absorption of the incident waves so as to reduce boundary reflection. The cortex is simulated as an isotropic and elastic solid layer. The soft tissue is simulated as acoustic fluid. The simulation parameters of cortex and soft tissue are shown in Table 1，and are similar to the parameters of human bones.

Table 1. Simulation parameters of cortex and soft tissue by the FDTD method.

A Gaussian envelope sinusoid was used as the excitation:

wherep(t) is the acoustic pressure of the excitation，Ais a constant value，σ=5.95×10-7andt0=1.5 μs. The center frequencyfwas 1 MHz for simulation 1 (free plate) and 0.5 MHz for simulation 2(bilayer model)，with the same values being used in thein vitroandin vivoexperiments，respectively.

3. Experiments

3.1. In vitro experiments

Four bovine long bones were bought in a supermarket and tested in thein vitroexperiments. After removing the soft tissue， the bovine bones were polished into plates. The thicknesses of the four bovine cortical bone plates were 2.78 mm，3.58 mm，3.83 mm，and 4.00 mm，respectively.

The setup of thein vitroexperiments is shown in Fig.1.A three-cycle Gaussian envelope sinusoidal wave was produced by an arbitrary waveform generator (Agilent 33220A，CO， USA). After amplification， this Gaussian envelope sinusoid signal was used to drive a transmitter with a central frequency of 1 MHz.

After propagating in the bovine cortical bone plate the waves were received by a 64-element array receiver. The pitch size，central frequency and-6 dB bandwidth for the array receiver were 0.675 mm， 1 MHz， and 0.1 MHz–1.5 MHz， respectively. The transmitter and array receiver were aligned on the same side of the bovine cortical bone plate， with 24.5 mm being the closest offset， and coupled by ultrasonic gel，which ensured that the transducers and bovine bones were in very good contact. The received signals were collected by a Vantage-64 system(Verasonics，WA，USA)with a sampling frequency of 31.25 MHz， which was also synchronized with the excitation signals. The 64 signals received for 2.78 mm bovine plate are shown in Fig. 2(a). Each bovine plate was measured 10 times under the same conditions.

Fig.1. Setup for the in vitro experiments.

Fig. 2. Received ultrasound signals: (a) in vitro experiment， (b) in vivo experiment.

3.2. In vivo measurements

Ten volunteers (six normal healthy and four with osteoporosis)who had been tested by DXA were measured at midtibial level with an ultrasonic axial transmission device. Informed consent was obtained from all volunteers who participated in this study， and local ethical committee approval was sought and attained.

Figure 3 shows the setup for thein vivomeasurements which were conducted in a different lab. Thus thein vivoexperimental setup is not the same as thein vitroone. In order to obtain signals with larger amplitudes from the human tibia， two transducers with a central frequency of 0.5 MHz were used as transmitter and receiver， and they were perpendicular to the tibia of the volunteers. Fixing the transmitter，the receiver was moved along the tibia with a closest offset of 22.5 mm. The moving step was 0.5 mm， and 60 receiving signals were scanned under computer control. Sampled at 10 MHz， and averaged 100 times， the received signals were collected by a PC-based digital oscilloscope. The experiment was done three times for each volunteer under the same conditions. A set of 60 ultrasound signals is shown in Fig.2(b).

Fig.3. Setup for in vivo measurements.

4. Results

4.1. Simulation results

Figures 4(a)and 4(b)show the results of simulation 1 and simulation 2. The black trajectories are the extracted dispersion curves by the Burg algorithm. The red lines are the theoretical dispersion curves based on the parameters of cortical thickness，longitudinal velocity and transverse velocity，which were obtained from the inversion algorithm. For simulation 1， the inversed cortical thickness， longitudinal velocity and transverse velocity are 4.0 mm， 3550 m·s-1and 1950 m·s-1，respectively. In comparison with the true values， the relative errors are 2.6%， 1.4%， and 2.6%， respectively. For simulation 2， the inversed parameters (CTh，VL，VT) are 3.5 mm，3600 m·s-1and 2050 m·s-1and the relative errors are 0%，2.7%，and 2.5%，respectively. Therefore，the proposed multiparameter inversion algorithm is able to estimate the cortical thickness， longitudinal velocity and transverse velocity accurately.

Fig.4.Comparison between extracted dispersion trajectories(black)and theoretical dispersion curves(red continuous lines)using estimated parameters for(a)simulation 1(free plate)and(b)simulation 2(bilayer).

4.2. In vitro results

As an initial example of thein vitroexperiments， figure 5(a)shows the results for the 2.78 mm thick bovine plate.The black experimental trajectories were extracted from the acquired ultrasonic time signals(Fig.4(a))based on the Burg method. Then the corresponding parameters (CTh，VL，VT)were estimated by the inversion procedure as 2.89±0.19 mm，3911±56 m·s-1，and 1817±23 m·s-1，respectively. The relative error for cortical thickness is 3.96%. The theoretical dispersion curves (red lines) based on the estimated parameters (2.89 mm， 3911 m·s-1， 1817 m·s-1) are also shown in Fig.5(a)，and the experimental trajectories are consistent with the theoretical dispersion curves. Figures 5(b)，5(c)，and 5(d)show the results for 3.58 mm，3.83 mm，and 4.00 mm bovine bones，respectively.

The results for estimation of cortical thickness， longitudinal velocity and transverse velocity for these four bovine bones are shown in Table 2. The real cortical thicknesses and the relative errors are also shown in Table 2. The root-meansquare error(RMSE)for cortical thickness(CTh)is 0.11 mm.Both the small error in cortical thickness and the matching of experimental dispersion trajectories with theoretical dispersion curves suggest that the proposed inversion method can accurately estimate the cortical bone parameters.

Table 2. Parameter estimations for the bovine plates and the relative errors for cortical thickness.

Fig. 5. Comparison between experimental dispersion trajectories (black) and theoretical dispersion curves (red lines) using estimated parameters of the invitro bovine cortical bones for different cortical thicknesses: (a)2.78 mm，(b)3.58 mm，(c)3.83 mm，(d)4.00 mm.

4.3. In vivo results

Table 3 shows the estimated parameters for tibiasin vivoby the UGW method. Cortical thickness，transverse velocity and longitudinal velocity are obtained from the multi-parameter inversion algorithm. The extracted dispersion trajectories (black)and theoretical dispersion curves (red continuous lines) based on the inversed parameters for the 10 volunteers are shown in Fig. 6. The results demonstrate that the experimental trajectories basically match well with the theoretical dispersion curves，which suggests the results of Table 3 are reasonable estimates of the three parameters for tibiasin vivo.

Table 3. Parameter estimation for tibias in vivo by the UGW method.

Table 4. DXA detection results for osteoporosis.

To test the practicability of the UGW method， the 10 volunteers are also underwent DXA detection; the results are shown in Table 4. Six volunteers were diagnosed with osteoporosis(T-score≤-2.5)and four volunteers were diagnosed as healthy normal(T-score＞-2.5)

A comparison of diagnosis of osteoporosis with the UGW method and DXA detection is shown in Fig. 7. Red circles denote volunteers with osteoporosis and blue triangles denote healthy normal volunteers， based on DXA detection.Figure 7(a) shows the osteoporosis results tested by inversed transverse velocity. Basically， the transverse velocities of healthy volunteers are higher than those of osteoporotic volunteers except for No. 3， which suggests that the transverse wave velocity is very sensitive to human osteoporosis. Figure 7(b)shows the osteoporosis results based on testing by inversed longitudinal velocity，which is probably not as sensitive as transverse velocity to human osteoporosis because the longitudinal velocities do not obviously distinguish osteoporotic and healthy normal volunteers. Figure 7(c)shows osteoporosis results on the basis of inversed cortical thickness testing.The cortical thicknesses of osteoporotic volunteers are thinner than those of the healthy volunteers except for No. 5， which suggests that when osteoporosis occurs the cortical thickness is reduced. Therefore， transverse velocity and cortical thickness are probably more sensitive to human osteoporosis than longitudinal velocity.

Fig.6.Comparison between extracted dispersion trajectories(black)and theoretical dispersion curves(red continuous lines)using estimated parameters for the tibias of 10 volunteers.

5. Discussion

In the simulation， two models (free plate and bilayer)were used to mimic single cortical bone and coated cortical bone. The estimated cortical thicknesses，longitudinal velocities and transverse velocities were nearly the same as the true values used in the simulations.This suggests that the proposed method is able to evaluate the bone parameters accurately.

In thein vitroexperiments，four bovine plates were used to validate the proposed inversion method. The experimental dispersion trajectories obtained by the Burg method were very clear. The theoretical dispersion curves based on the estimated CTh，VL， andVTmatched well with the experimental ones. Compared with the reference values， the biggest relative error in these four cortical thicknesses was 4.25%，and the RMSE for cortical thickness was 0.11 mm. These suggest that although real cortical bone is anisotropic，the proposed inversion method was able to estimate the bone parameters successfully.

Compared with simulations andin vitrobone plate experiments， ultrasonic signals collected fromin vivotibia suffer much more noise due to the coated soft tissue layer.Therefore，the extracted experimental dispersion trajectories are much noisier for thein vivoexperiments. Nevertheless，the theoretical dispersion curves for the estimated parameters(CTh，VL，andVT)are consistent with the experimental dispersion trajectories，which suggests that the proposed inversion algorithm is also able to evaluate bone parametersin vivo.

Velocities are important parameters for characterizing dispersion of guided waves.[31]When longitudinal and transverse velocities change， the corresponding dispersion curves will change. Therefore， based on the different dispersions of guided waves，the changed longitudinal and transverse velocities could be inversed. When osteoporosis occurs in older people the holes are much larger than in healthy bone， which means that the cortical layer becomes porous(filled with marrow) and less dense. As is well known， longitudinal waves can propagate in solids，liquids and air; however，the longitudinal velocity in a solid is higher than the other two media.On the contrary，transverse waves can only propagate in solid media. Therefore， the longitudinal and transverse velocities for osteoporosis patients should be lower than in health， and the change in transverse velocity should be larger than that of longitudinal velocity because transverse waves cannot propagate through the liquid or air which fills the much larger holes when osteoporosis occurs. Although it is hard to obtain someone’s transverse velocities before and after the occurrence of osteoporosis simultaneously， this was partially demonstrated in our study. The transverse velocities of volunteers with osteoporosis were lower than those of healthy volunteers except for No.3(Fig.7(a));this means that even for different volunteers，the transverse velocity changes a lot when osteoporosis occurs， so that a clear transverse velocity gap could be observed between osteoporotic and healthy volunteers. On the contrary， longitudinal waves propagate not only in solids but also in liquids and air. When osteoporosis develops， the longitudinal velocity will change， but not as much as the transverse velocity. Thus， for different volunteers with different bone masses， the longitudinal velocity is not as good a good biomarker for distinguishing between osteoporotic and healthy bones(Fig.7(b)).

Fig.7. Osteoporosis: between DXA and UGW methods: (a)transverse velocity results， (b) longitudinal velocity results， (c) cortical thickness result.Red circles denote volunteers with osteoporosis and blue triangles denote healthy volunteers.

Many researchers have listed cortical thickness as being a key point in their studies[8，21，35]because it will change when osteoporosis develops in elderly people，and this is also shown in our study(Fig.7(c)). The cortical thicknesses of all six volunteers with osteoporosis was less than 4.0 mm， which suggests that when osteoporosis occurs the cortical thickness becomes thinner than before. Therefore，transverse velocity and cortical thickness could be good biomarkers for differentiating osteoporosis; combining them to diagnose osteoporosis may be a good choice in clinical applications.

There are some important limitations in our study. Only 10 volunteers were recruited. More subjects are needed to test the relationship between osteoporosis and parameters that are inversed by the UGWs. However，there is a clear trend that the inversed transverse velocity and cortical thickness are highly correlated with human osteoporosis. The multi-parameter inversion algorithm could be computationally expensive because it compares experimental dispersion trajectories with all possible theoretical dispersion curves. In our study， it will run for about 15 min to solve the bone inversion procedure for one subject on a computer with an Intel Core i7 4.00 GHz processor. This is sufficient for an offline study; however， it is not acceptable for use in real-time bone diagnosis. Thus，in future work， we will improve the inversion procedure to lower the computational cost. Apart from that， real bone is an anisotropic irregular hollow shell， a middle layer between tissue and marrow. Although an isotropic free plate model is suggested for acquiring cortical bone properties evenin vivo，a proper bone model is worth studying in future work.

6. Conclusion

In this research， a multi-parameter inversion algorithm based on UGWs is introduced to estimate the cortical thickness， longitudinal velocity and transverse velocity of long bones. The feasibility of the proposed algorithm is verified using simulations，in vitrobovine plates andin vivotibias. The results demonstrate that the multi-parameter inversion algorithm was able to evaluate cortical thickness and bulk velocities (VTandVL) accurately. Besides， compared with DXA detection， transverse velocity and cortical thickness are more sensitive to human osteoporosis than longitudinal velocity，and these two parameters could be good biomarkers for divide osteoporotic and healthy bone. Therefore，these two parameters(CTh andVT)of long bones can potentially be used to assess bone status in clinical applications.

Acknowledgments

Project supported by the National Natural Science Foundation of China (Grant No. 12034005)， in part by the Program of Shanghai Academic Research Leader (Grant No.19XD1400500)，and in part by the China Postdoctoral Science Foundation(Grant No.2019M661334).