Ultra-wide band radar for prospective respiratory motion correction in the liver

Objective. T1 mapping of the liver is time consuming and can be challenging due to respiratory motion. Here we present a prospective slice tracking approach, which utilizes an external ultra-wide band radar signal and allows for efficient T1 mapping during free-breathing. Approach. The fast radar signal is calibrated to an MR-based motion signal to create a motion model. This motion model provides motion estimates, which are used to carry out slice tracking for any subsequent clinical scan. This approach was evaluated in simulations, phantom experiments and in vivo scans. Main results. Radar-based slice tracking was implemented on an MR system with a total latency of 77 ms. Moving phantom experiments showed accurate motion prediction with an error of 0.12 mm in anterior-posterior and 0.81 mm in head-feet direction. The model error remained stable for up to two hours. In vivo experiments showed visible image improvement with a motion model error three times smaller than with a respiratory bellow. For T1 mapping during free-breathing the proposed approach provided similar results compared to reference T1 mapping during a breathhold. Significance. The proposed radar-based approach achieves accurate slice tracking and enables efficient T1 mapping of the liver during free-breathing. This motion correction approach is independent from scanning parameters and could also be used for applications like MR guided radiotherapy or MR Elastography.


Introduction
Liver diseases were the main cause for about 2 million deaths in 2019 alone, adding up to 3.5% of all recorded deaths worldwide (World Health Organization 2020). The most prominent case being liver cirrhosis which develops from fibrotic tissue. The progression of liver fibrosis may be reversed if detected in early stages (Farci et al 2004) therefore reducing the risk of liver failure and primary liver cancer (Davis and Roberts 2010). Quantitative MR imaging such as T1 mapping or MR elastography has proven to be an effective tool for detecting tissue abnormalities in the liver, providing biophysical parameters that are comparable over time and between patients (Kelekis et al 1996, Heye et al 2012, Kim et al 2012, Manduca et al 2021, Babu et al 2016. MR data acquisitions, especially for quantitative imaging, are time consuming and hence prone to motion artefacts (McClelland et al 2013). This is particularly true for the abdomen where the liver can move up to several centimeters due to respiratory motion (Suramo et al 1984). To avoid breathholds during scans, which restricts the time of a single acquisition and can vary in absolute positions in between measurements (Takamatsu et al 2013), motion tracking becomes necessary. These motion estimates can be used either for gating (Ehman et al 1984), where data is measured only within a certain motion state or the motion estimate is used to adapt the measurement in such a way that displacements are compensated for (slice tracking) (Danias et al 1997). The latter having the advantage of acquiring data continuously and therefore measuring more data in the same Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. amount of time compared to gating. However, such approaches require motion models that uses surrogate motion measurements as input and produces a motion estimate as output (McClelland et al 2013).
Different motion surrogates for slice tracking have been proposed. For example, pencil beam navigators (Ehman and Felmlee 1989) are widely used to track the position of the diaphragm with fast 1D MR acquisitions. However, this approach interrupts the MR image acquisition and affects the magnetization of the liver. This can be overcome by using other surrogate signals such as the pilot tone (Speier et al 2015) or a noise navigator (Navest et al 2018) which track modulations in a generated constant RF signal and the statistical noise in the receiver channels respectively. These methods can generate surrogate signals only when MR data is received beforehand. This can restrict sequence designs if the sequence needs to be modified in order to extract the surrogate for prospective correction. This problem can be overcome by a system which is decoupled from the MR system and hence provides a motion surrogate independent from any MR data acquisition or MR hardware.
Previous work has shown that a radar surrogate signal can provide information from internal organ movements with very high temporal resolution (Schwarz et al 2010). The advantage of the radar approach over commonly used surrogates like respiratory bellows (Ehman et al 1984) lies in its ability to provide signals not just from the body surface, but also from regions below. Additionally, the radar signals are obtained contactless without the need to attach any devices to the patient directly.
In this paper we present radar-based slice tracking (RST), an implementation for accurate prospective motion correction. RST consists of an external program that controls an ultra-wide band radar system, calculates a motion estimate and sends a correction to the running MRI sequence in addition to creating a motion model based on a radar signal calibration. We evaluate the method based on robustness to different motion types, accuracy and long-term stability via simulations with ground truth data and measurements on a moving phantom with controlled parameters. We compare the performance of RST to a respiratory bellow correction as reference in in vivo measurements.

Methods
The proposed prospective RST approach consists of two main parts, as depicted in figure 1. First, a motion model which describes the movement of a region-of-interest (ROI) based on the current radar signal is estimated prior to the diagnostic scan. The radar signal is recorded simultaneously to a short dynamic MRI scan. After selecting a ROI in the MR images a linear motion model λ is formed describing the movement of the ROI as a function of the radar signal. In a second step, λ is then utilized together with a real-time radar signal during a diagnostic scan to update the current slice position. These two parts are explained in more detail in the following section. In the calibration step the radar trajectory δ gets extracted using principal component analysis (PCA) and combined into a linear motion model with the displacement j estimated from an ROI (red box) selected in motion-resolved images. This motion model is used in the following scans to predict displacement of the ROI from a radar response. Based on this displacement the slice position is adapted in real-time during data acquisition to correct for the respiratory motion.

Radar motion surrogate
An array of transmitting and receiving ultra-wide band antennas (1-5 GHz) is positioned above a subject's torso inside the scanner bore. With this array, radar pulses are sent towards the subject. At every interface where volumes of different dielectric properties meet (e.g. air, skin, fat, muscle, liver, lungs), part of the radar signal is reflected and another part continues through the tissue. The sum of all reflections gives the impulse response function ρ(t) for one radar pulse with a duration of 57 ns. It captures information about the subject's internal morphology where contributions which occur later in ρ(t) are expected to originate further away from the antennas including deeper tissues (Thiel et al 2009). By measuring these impulse response functions over time (i.e. with a sampling rate of 84 Hz), changes in ρ(t) due to changes in the tissue interfaces can be identified. It is assumed that the largest amount of change in the time series of these impulse response functions is caused by respiratory motion. It is possible to identify prominent modulations in such a time series using principal component analysis (PCA) (Pearson 1901). The largest components in the series can therefore be identified as respiratory motion. The coefficients of these principal components are used as a motion surrogate δ. A combination of multiple components has been evaluated as part of the experiments.

Motion estimation j
A set of dynamic sagittal MR images is acquired simultaneously with the radar signal. In the reconstructed images an ROI is chosen, (e.g. a tumor/cirrhosis or liver vessels in healthy volunteers). This region defines the area which will be locally corrected by the model λ. The image registration is done by maximizing the crosscorrelation to a reference image. Subpixel translation is estimated by further upsampling the Fourier transformed image in a small neighbourhood around the initial maximum (Guizar-Sicairos et al 2008, Van der Walt et al 2014). For the two main directions of respiratory motion displacement (anterior-posterior and headfeet) a separate linear motion model is built with the simultaneously acquired radar response.

Calibration
After extracting the motion surrogate δ and the registered displacement j in the reconstructed images a motion model λ can be constructed For motion correction a reference respiratory motion state is selected by splitting the respiratory motion into eight states based on the surrogate signal j and choosing the motion state with the highest number of data points. This reference point defines the zero position of the motion model all other positions will be corrected towards. Due to variances in logging times, a temporal delay between δ and j is introduced, which can be corrected for by introducing a shifting term that maximizes the cross correlation between both signals. The radar data is then averaged over the same time intervals used for the MR image reconstruction and the two synchronized data arrays are combined to form a motion model λ. This model is linearly extrapolated and uncapped on any side.

Slice tracking
As depicted in figure 2 in a regular setup the MR host controls the subsequent sequence planner which in return controls the gradients during the scan. In RST we use a standalone program on the MR host which is connected to the radar system and controls the radar measurements. Once the system is activated a stream of radar response functions ρ is recorded and combined with the coefficients of the PCA and used in the motion model λ to predict the current respiratory motion. Both were established in the previous calibration step. The resulting correction parameters are sent to the running sequence via vendor specific socket connections. For prospective motion correction the current slice position is adapted during sequence execution in real-time.

Experiments
To evaluate the proposed motion correction scheme, numerical simulations, experiments with a motion phantom and in vivo scans were carried out.

Radar system
An array of one transmitting and two receiving double-ridged horn antennas (Schwarz et al 2010) are positioned in two locations above the subject's torso inside the scanner. After emitting a radar signal generated by an M-sequence (Sachs 2003) baseband module (MEODAT GmbH, clock signal = 8.95 GHz) the reflected signals are recorded. Using a cross correlator within the hardware, the impulse response function ρ can be calculated. This function describes the wave propagation in a given morphology over a total timespan of 57 ns. It consists of 511 datapoints, each with a dwell time of 112 ps. Impulse response functions are generated by correlating the original pseudo-random sequence with the reflected signal at 512 different offsets. These functions are then averaged and recorded with a sampling rate of 84 Hz allowing for an accurate tracking of respiratory motion. Typical signal to noise ratios for the acquired radargrams are around 100.

MR imaging
All experiments were done on a 3 T system (MAGNETOM Verio, Siemens Healthineers, Erlangen, Germany) using inhouse-developed sequences and reconstruction methods. A 4-channel Siemens Large Flex coil and an 8-channel Siemens Spine coil were used for signal reception. The dynamic MR scan for calibration was carried out via a spoiled gradient echo sequence with a 2D golden angle radial trajectory (Winkelmann et al 2006) with 7000 spokes in sagittal view (TR/TE/FA 3.9 ms/1.73 ms/12°, voxel size = 1.2 × 1.2 × 4 mm 3 , FOV = 300 mm 2 ) with a total scan time of 27 s. Dynamic images were reconstructed using a non-Cartesian iterative SENSE method (Pruessmann et al 1999(Pruessmann et al , 2001 while the time resolution for each image was optimized as shown in the experiments. Measurements for T1 maps used the same trajectory with 3360 spokes for 16 s with an inversion pulse every 2276 ms (Becker et al 2017). (TR/TE/FA 4.6 ms/2.03 ms/5°, voxel size = 1.2 × 1.2 × 6 mm 3 , FOV = 300 mm 2 ). T1 maps were created using a model-based acceleration of parameter mapping

Numerical simulation
To simulate the interaction between the radar and the body tissue a commercial finite-difference time-domain (FDTD) environment was used (RemCom XFDTD 6.4). For simplification the radar transmission was approximated by a Gaussian pulse applied to a Vivaldi antenna (Gibson 1979), which has similar properties to the used double ridged horn antenna but is easier to implement into the simulation framework. The penetration depth was analyzed on a simple multi-layer tissue phantom of varying thickness. The signal loss was simulated for skin, fat and liver tissue. XCAT phantoms were used for motion-resolved simulations for different breathing types. The 4D extended cardiac torso software package (Segars et al 2010) allowed to generate a virtual human phantom in 13 different respiratory states. These were formatted into FDTD compatible meshed objects using empirical dielectric tissue properties (Gabriel 1996). The robustness of motion estimates was tested by evaluating the motion model error for different breathing types. Therefore, we created virtual body phantoms describing both extreme cases of possible motion, one being motion in only head-feet direction the other only in anterior-posterior direction. The effect of antenna positioning was investigated with two receiving antennas at different locations (see figure 2).

Phantom experiments
An agarose phantom was mounted on a small wagon with a fixed track. The phantom moved back and forth along the scanner bore direction, driven by a pushrod with some mechanical clearance. The phantom trajectory is similar to a sinewave. To simulate respiratory motion, the amplitude of the motion of the phantom was set to 25 mm with a cycle length of 10 s. The phantom ensured a reproducible motion pattern over a long time.
The goal of the phantom measurements was to evaluate the accuracy of model predictions, finding the optimum number of principal component weights to combine into a motion surrogate and investigate the longterm stability of a motion model. Five dynamic MR scans were carried out over a total time span of 1.8 h. In all measurements λ was created based on the motion information obtained from dynamic MR images and the simultaneously acquired radar signal. Accuracy was accessed by comparing prediction of λ to the motion information from the dynamic MR images for each time point individually. Long-term accuracy was evaluated by creating a motion model once and testing its performance in all subsequent scans to identify possible drift in the setup caused by heating of components or signal attenuation.
Using combinations of multiple PCA coefficients different motion surrogates were created. The weightings of each coefficient in this sum were optimized to minimize the error to the known trajectory of liver displacement. Multiple motion models were created for one timepoint and their errors were evaluated at all following timepoints. This gives an estimate on which coefficients provide additional motion information to further refine the motion model and which just give unrelated information, which lead to an overfitted model that does not generalize. The RMSE was used as a metric to compare motion models with different surrogates.
For an optimal λ each MR image should be reconstructed from radial k-space data collected over a time frame as short as possible to get a high temporal resolution. At the same time each image should include enough data to ensure high image quality for accurate motion estimation, thus finding the minimum number of radial lines necessary to still perform image registration was a goal of the experiments. Images with varying numbers of radial k-space lines in each dynamic frame were reconstructed and motion models were estimated. The accuracy of the motion models was evaluated to find an optimal number of radial lines for image reconstruction.
The total time delay between respiratory motion and application of the corresponding slice correction was measured by having the phantom move along one image dimension and perform correction along the orthogonal non-moving image dimension. In image space both trajectories were extracted separately by image registration and cross correlated to determine their temporal delay. All comparisons between motion models were based on the root-mean-square-error (RMSE) between the motion registered in the dynamic images and the corresponding prediction based on the motion model.

In vivo scans
Data was measured from two healthy volunteers aged 30.5 ± 3.5 years. The participants gave written informed consent in advance according to regulations of our institutions ethical committee.
This data acquisition was carried out during free-breathing with and without the proposed motion correction approach and during a breathhold for reference.
As a comparing signal for in vivo measurements, we used a Biopac TSD221-MRI pneumatic respiration transducer on a Biopac MP150 acquisition system. The data was recorded with a sampling rate of 0.5 ms. Figure 3 shows the simulated propagation of a radar pulse through a tissue layer. Reflections occur at both the air-tissue-and the tissue-air interface. The radar signal loses intensity the further it penetrates the dielectric material. The simulations showed that the received signal intensity was halved every 10.19 cm it travels through fatty tissue. For skin and liver tissue the signal intensity was halved already every 1.5 cm (see figure 3). Figure 4 shows the performance of the motion model for different numbers of PC coefficients used for the motion surrogate δ. Motion models were built using combinations of PC coefficients ranging from one up to twenty different coefficients used. These models were calibrated at timepoint t = 0 and evaluated at five different later times during the phantom measurements. The mean of RMSE score over all time points is depicted in figure 5. The best model accuracy was achieved by using coefficients from just the first PC. All further experiments therefore used the first PC coefficient to build the motion models.

Radar response
In figures 5(a), (b) radar responses ρ(t) were recorded from two antennas, one upper antenna located closer to the subject's head and one sideward antenna closer to the side of the body (see figure 2). Both signals depict a reflection profile over 57 ns consisting of 511 datapoints, each of these points is treated as a separate channel which changes over the course of different respiratory states. It shows stronger contributions towards the first principal component early on in the signal but still continues to find information in channels throughout the entire response as indicated in red in figures 5(a), (b). Even though the signals from both antennas differ in shape the overall changes in these signals describe the respiratory motion equally well. In figures 5(c), (d) the extracted radar motion surrogate trajectories can be seen for the two different types of respiratory motion -one being motion in only anterior-posterior direction, the other only in head-feet direction. The recorded surrogate trajectories follow the ground truth liver displacements in both respiratory motion types. The motion model in the case of only head-feet motion had a slightly smaller error, but both motion models produce an R 2 close to 1.   Figure 6 shows the model error for motion models calibrated with reconstructed MR images of different temporal resolution. High temporal resolution (b) leads to reconstructed images with large undersampling  artefacts, while low temporal resolution (d) is unable to resolve the phantom movement. The model error drops below 1 mm when the number of lines per image is chosen between 8 and 47. The phantom data shows a local minimum at 12 lines and a steep increase in error above 47 lines (see figure 6). Figures 7(a), (b) show the difference between measured and predicted tissue displacement in the moving phantom and in an in vivo experiment. This experiment was done at multiple times, the resulting average RMSE was for the phantom 0.12 ± 0.01 mm for the anterior-posterior model and 0.81 ± 0.07 mm for the head-feet model. The in vivo experiments resulted in 0.52 ± 0.22 mm for the anterior-posterior model and 0.88 ± 0.31 mm for the head-feet model.

Accuracy and long-term stability
In figures 7(c), (d) measurements on long-term stability of a motion model calibrated at timepoint t = 0 are shown. Evaluations at different time points in the phantom show an error centered around 1 mm in the headfeet model and around 0.15 mm in the anterior-posterior model. This error seems to be consistent for all times up to almost two hours without any obvious drift -Also in the in vivo experiment no drift or an increase in error was seen.

Total delay
In order to determine the delay between the motion and the slice tracking, the motion correction was applied orthogonally to the movement. Figure 8 depicts the trajectories of a phantom physically moving along one axis and the corresponding motion correction applied to a perpendicular axis. Correlating the separated motion trajectory with the correction trajectory resulted in a correction latency of about 77 ms.   Figure 9 compares trajectories from a radar motion surrogate and a respiratory belt surrogate to the registered liver displacement in an in vivo experiment. The RMSE of RST with 0.84 mm is comparable to the error in phantom experiments. A motion model based on a respiratory belt signal resulted in an error almost 3 times higher with a trajectory visually different from the registered liver motion. The reference trajectory was heavily over-/undershot at certain points in time. Figure 10 shows in vivo correction measurements of free-breathing data acquisitions, without and with RST, and a breathhold measurement. In the free-breathing scan displacements in both image directions were visually reduced and the T1 map looks similar to the breathhold acquisition. However, breathhold and corrected acquisition should not be compared directly since they depict different liver positions and the breathhold shows a small drift in addition. Figure 11 depicts the quantitative T1 maps from an uncorrected free-breathing acquisition, a free-breathing acquisition with RST and a breathhold scan. The RST data shows sharper edges for free-breathing, especially in the liver vessels that were used for calibration but also on the upper liver dome and in the kidneys as indicated with white arrows. Image details from RST are comparable to the breathhold data.

Discussion
The proposed RST approach is able to predict displacements in the liver due to respiration and correct for them in a prospective way. The experiments in the simulation environment have shown, that the radar response originates from body layers of different depth and the influence of adipose tissue to signal intensity is only minor. The proposed approach showed an accurate motion prediction in phantom scans with an MAE of 0.6 mm over 2 h which was lower than for a previously proposed pilot tone-setup (Ludwig et al 2021) where errors of 2.5 mm were reported for a scan time of less than 9 min. In vivo experiments suggest that a motion model based on a  Figure 10. In vivo measurements. The trajectory of a blood vessel is depicted in a projection over time in both image directions. In the uncorrected case displacements due to respiratory motion is clearly visible. The corrected case shows a measurement during which radar correction was prospectively applied. The breathhold is an uncorrected measurement as reference.
radar response function results in better motion prediction than a model based on a respiratory belt, especially during times of under-/overshooting the actual displacement where the respiratory belt model delivers untrustworthy predictions. Measurements with prospective RST shows a successful compensation of liver displacements and visibly improved quality in acquired T1 maps. The in vivo data show that the RST approach can successfully be implemented on a running MR system, but the results are limited by the fact that it was tested on just two volunteers and no ground truth motion can be accessed in in vivo scans. Quantitative assessment of resulting T1 maps was not done here and goes beyond the purpose of this work.
The radar signal is attenuated by the tissue. Nevertheless, electrical permittivity of fatty tissue is low, thus less energy is taken from the radar waves to move charges within the tissue. Figure 3 shows that even after 10 cm of fat tissue, still 50% of the signal amplitude is available, making RST usable for patients of various body mass indices (BMIs). This is also supported by figure 5 which shows that the PCA detects signals related to the respiratory motion not just at the beginning of the response function (i.e. first reflection on the body surface) but also contributions from later time points. These contributions can be interpreted as reflections from dielectric interfaces inside the body. The strong signal attenuation in liver tissue however indicates that volumes (like small tumors) deep inside the liver are likely to not be measured directly. Since they do not move autonomously inside the liver it is sufficient to detect changes in the liver boundaries and interpolate movements with the linear calibration.
Antennas positioned at different locations record varying radar responses, since they are exposed to other morphologies. The predominant changes in these signals however stayed similar. Thus, using multiple receiving antennas at once may only lead to a small improvement of the surrogate signal for liver displacement.
Future MR measurements on the heart might cause a risk of signal shielding since cardiac coils are typically denser than the flex coils used in this work.
Regarding the variability of different breathing motions themselves, the simulations have shown that the radar system provides reliable motion correction models in both extreme cases of only anterior-posterior and only head-feet motion. The prediction error was slightly lower in the head-feet model even though the motion change was parallel to the antenna opening. This indicates that the source of these changes originated from deeper tissue interfaces which provide a smaller signal intensity but describe the desired motion more accurately. On the other hand, displacements in head-feet direction have a higher amplitude and would therefore lead to a more distinctive change in the radar response. The tissue interfaces are not perfectly planar either, so a displacement in only head-feet direction will result in a small change in anterior-posterior direction simply due to the difference in elevation. Increasing the number of principal components used as a motion surrogate generally increased the model error. When combining components their linear factors were optimized to match the registered displacement in the calibration best. These factors do not seem to generalize for later scans which suggests that information useful to predict motion in the presented way are restricted to only the first or perhaps up to the first three components. In vivo studies may show a more complex interaction of different motions in the first principal components. To isolate the respiratory part of this motion a frequency bandpass filter can be useful, as it has been proposed in prior work (Thielemans et al 2011).
For the reconstructed images to capture respiratory motion a sampling frequency of at least twice the respiratory frequency is necessary (Nyquist-Shannon sampling theorem). Ideally the sampling frequency is chosen as high as possible to resolve the motion trajectory with most accuracy. Phantom experiments show that images reconstructed with less than 12 radial lines increase the model error due to dominant undersampling artefacts. Since features in in vivo experiments are less defined than in a phantom, 23 lines per image were chosen which still allowed good registration while providing a temporal resolution of 89.4 ms.
Experiments on a moving phantom have shown, that the RMSE of a motion prediction lies in the range of about 5% of the total movement amplitude. The error caused by the delay of the entire motion prediction pipeline (80 ms) could be further reduced by incorporating a forward prediction, e.g. in the form of a Kalman filter (Spincemaille et al 2008, Bacher et al 2018.
A drawback of the having the radar antennas in the scanner lies in the fact that its current implementation restricts the 70 cm bore by an additional 10 cm from the top and 4 cm on each side due to the way they are mounted. Tissue heating for the radar is very unlikely with an RF output power of 1 mW which is magnitudes lower than Wi-Fi e.g. from mobile phones not to mention way lower than the MR system itself. The antennas operate outside the MR frequency of 123 MHz. Cables connect them to the signal generators which are positioned outside the MR scanner room to prevent any RF noise degrading image quality.
Here RST was used to enable T1 mapping during free-breathing but it could be applied to many other approaches such as MR elastography (Muthupillai et al 1995). Even radiotherapy (Raaymakers et al 2008) is possible because RST is independent from the scanner and does not need an MR acquisition running to predict tissue motion.

Conclusion
This work presents a prospective respiratory motion correction method utilizing an ultra-wide band radar signal as a surrogate signal. A short calibration scan is used to obtain a fast quantitative motion signal describing displacement in a selected region of interest in the body. The model corrections have shown to be accurate, suitable for different motion types, as well as stable over time. First in vivo applications are promising as they show an improvement in image quality for free breathing T1 mapping acquisitions.

Ethical statement
The research related to human use complies with all the relevant national regulations, institutional policies and was performed in accordance with the tenets of the Helsinki Declaration. The examination was approved by the local ethics committee ('Ethikkommission der PTB') and is in accordance with the relevant guidelines and regulations. Written informed consent for examination and publication was received from all patients prior to the examination. No approval ID number was specified.

Conflicts of interests
Peter Speier is an employee of Siemens Healthcare, Erlangen (Germany).

List of patents
UWB antenna system (DE102008019862B4).