1 INTRODUCTION

MRI relies on highly accurate gradient fields for signal preparation and spatial encoding. Deviation from prescribed gradient waveforms causes image artifacts and error in quantitative readouts. Much effort is thus dedicated to perfecting gradient hardware and correcting remaining error at the levels of input waveforms, image reconstruction, and image processing.

One chief cause of perturbation is eddy current driven by gradient switching. Eddy currents outside the gradient tube, especially in the magnet and cryostat, can be largely suppressed by active shielding.1-3 In contrast, eddy currents in the gradient coils themselves generally remain and are a key consideration in design optimization.3 Relevant eddy currents can also occur in RF equipment, particularly in RF screens.4-6 Remaining eddy currents are commonly addressed by pre-compensation filters (also termed pre-emphasis) applied to input waveforms, exploiting the linear dependence of eddy currents on the underlying field dynamics.7-13

Besides eddy currents, gradient switching also drives mechanical oscillations, or vibrations, of gradient tubes,14-18 which are particularly strong at mechanical resonance frequencies 19-22. The mechanical behavior of gradient coils has been extensively analyzed, modeled, and accounted for in gradient design.15, 23, 24 As part of these efforts, the acoustic response of gradient systems to input waveforms has been described by linear models, exploiting approximate linearity of electromechanical coupling.15, 21, 25 Naturally, besides acoustics, the mechanical behavior of gradient systems is equally central to considerations of structural integrity. However, it can also contribute to field perturbation via reverse coupling from mechanics to electromagnetics. One prominent effect of this kind is magnetic field oscillation along with mechanical resonance,19 which also alters of coil impedance.26

To capture direct electromagnetic as well as mechanical pathways of field perturbation, gradients have been modeled as general linear time-invariant (LTI) systems.27-30 Such models have proven useful as a basis of pre-emphasis,31 for image reconstruction,28, 32 and for studying field perturbation by devices placed within the gradient range.30

LTI modeling covers diverse gradient behavior and achieves high degrees of accuracy. But the underlying assumption of time invariance cannot be expected to hold strictly. It has been reported that LTI models can retain their utility over long periods of time up to three years.32 However, gradient coils do undergo physical changes in the short term as they heat up during operation. Gradient heating is predominantly due to ohmic losses in the coils, from both bulk and eddy currents, and forms another key consideration in gradient system design.33-36

Regarding electromagnetic behavior, the chief effect of gradient heating is increase in resistivity. Change in overall resistance hardly alters the net current flowing through the coil terminals when using gradient amplifiers with high output impedance17 or feedback control based on sensing of output current.37 However, change in resistivity due to change in temperature does alter eddy currents within coil conductors and in other system parts subject to heating. The mechanical properties of gradient tubes also change with temperature, forming another potential cause of LTI violation.29 A mechanism of particular interest is thermal shift of mechanical resonances, which has been suggested as a cause of time-dependent artifact in EPI time series.38

Against this background, the goal of the present work is to study the effects of gradient heating at the level of system response and to expand the LTI model to include thermal variation. To capture changes in response, measurement of gradient transfer functions is performed rapidly, in 1 s per frame, and at high spectral resolution. Thermal modeling is performed by linear expansion of transfer functions in terms of temperature observables. It is found that this approach permits prediction of thermal response changes due to both change in the lifetime of eddy currents and shifts of mechanical resonances. A global linear model is found to perform virtually as well as tiling the temperature domain with local models. Parts of this work have previously been presented at ISMRM conferences39-41 as referenced in the following.

2 METHODS 2.1 Hardware

All experiments were carried out on a 3T whole-body MRI system (Achieva, Philips Healthcare, Best, The Netherlands), equipped with actively shielded gradients of maximum strength of 40 mT/m and slew rate of 200 mT/(mms). In the gradient tube, the radial order of coils is x, y, and z from inside to outside and the z coil and shield are water cooled. Throughout this study, the system’s built-in eddy current compensation was enabled.

Magnetic field responses were recorded with a 16-probe field camera (Skope Magnetic Resonance Technologies, Zurich, Switzerland) that yields magnetic field evolution in terms of spherical harmonics up to third order.42

Thermal imaging of the gradients in operation was performed with an infrared camera (Compact, Seek Thermal, California, USA) mounted on a mobile phone. The camera was placed 2.5 m from the patient end of the MRI system at an angle of to capture at least a quarter of the bore wall. For direct view of the gradient tube, the bore liner and the RF body coil were removed. Temperature measurements were performed with a total of 19 fiber-optic sensors (Neoptix, Quebec City, Quebec, Canada) with temporal resolution of 1 Hz and relative accuracy specified as C.

2.2 Heating and temperature measurement

The spatial patterns of gradient heating are frequency dependent. To capture frequency differences, gradient heating was induced separately by DC and AC operation of one gradient channel at a time, separated by at least 15 min for approximate cool-down. DC heating was performed over 5 min, using 1.5 s intervals of constant amplitude at 25 mT/m separated by short breaks due to software limitations. AC heating was achieved with 4-minute EPI sequences of 31 mT/m peak gradient strength and variable fundamental frequency of the readout gradient. Low, medium, and high frequencies were targeted by variation of scan parameters at the limit of the gradient specifications, avoiding mechanical resonance frequencies. The resulting frequencies were 355 Hz, 383 Hz, and 390 Hz for the x, y, and z gradient, respectively (low range), 761, 880, and 751 Hz (medium range), and 1029, 1143, and 999 Hz (high range). The DC and medium-frequency protocols were used for training response models and will be referred to as “training sequences.” The low- and high-frequency protocols were used only for validation across frequency bands and are referred to as “cross-validation sequences.”

Placement of the temperature sensors was based on initial infrared imaging during the 12 heating scenarios,41 each followed by subsequent cooling down. Heating via the z gradient led to only minor temperature change due to the active cooling of this coil. To cover a somewhat greater temperature range, the z heating sequences were therefore performed twice in immediate succession.

The resulting infrared image series, taken at 20 s intervals, were analyzed to determine sensor positions that capture the chief temperature dynamics. In a first step, candidate sets of 13 positions were chosen randomly and the covariance matrix of the 13 temperature time courses was calculated for each set. As a metric of temperature variation captured, the determinant of the covariance matrix was computed. From 100,000 random sets of positions, that with the largest determinant was selected. In a second step, this set was refined by maximizing the determinant by local gradient ascent. At the resulting positions on the inner surface of the gradient tube, 13 temperature sensors were mounted. Another five sensors were mounted on the water cooling circuits (one on the inflow and four on different outflow channels) and one was used to record room temperature.

2.3 Measurement of transfer functions An LTI system is fully determined by its impulse response h(t), which relates the system input, i(t), to its output, o(t), by a convolution: (1)In the frequency domain, convolution simplifies to multiplication and deconvolution to division: (2)where denotes the angular frequency, and the Fourier transforms of the input and output, and the system’s transfer function. Viewing a gradient chain in these terms, its output is the amplitude of the gradient field that it generates while its input may be considered at different levels such as the coil terminals, the analog input to the power amplifier, or the preceding digital waveform definition. The latter perspective is chosen in the present work and inputs are defined via the system console. Gradient transfer functions were measured as described in Ref. [29] using a frequency-swept pulse.43 The sweep pulse c(t) was of the form (3)with the linearly increasing frequency (4)with bandwidth BW = 30 kHz, pulse length = 40 ms, and amplitude A(t) starting at 31 mT/m and gradually decreasing with increasing frequency to stay within slew rate constraints. The pulses were played out via the MR system console, which was also programmed to trigger acquisition with the field camera.

The duration of the field readout determines the frequency resolution of resulting transfer functions. To extend the effective readout duration beyond the lifetime of the field probes, data from multiple successive acquisitions were concatenated. The sweep pulse was played out four times in a row with a TR of 250 ms, triggering field acquisition at different times as shown in Figure 1. After each excitation, the field probe signals decay due to relaxation, causing corresponding decline in SNR and noise increase in field estimates (Figure 1, right). For smooth concatenation, adjacent measurements were taken with 10 ms of overlap. With this procedure, an effective acquisition duration of 200 ms was achieved, yielding frequency resolution of 5 Hz based on data collected in just under 1 s. The gradient transfer function was obtained by frequency-domain division of the gradient field output and sweep input (Equation 2).

Input and field measurement during and after a 40-ms long sweep pulse. The SNR of the field probe readouts drops due to decay of the probe signal, causing increase in field noise (right). Therefore, four measurements are concatenated at , 70, and 130 ms, yielding the transfer function at a frequency resolution of 5 Hz in less than 1 s

Time series of transfer functions were collected after each of the aforementioned heating protocols, for 8 min and once every 10 s, which was found to suffice based on the dynamics observed. Temperature was monitored throughout. Temperature and transfer function recordings were synchronized based on the console trigger, using linear interpolation to align the temperature data.

2.4 Thermal model To incorporate thermal behavior, the LTI system model is expanded by temperature dependence of the impulse response and, equivalently, of the transfer function. In doing so, it is assumed that the system is time invariant in the short term, at the scale of the duration of its impulse response, but changes its characteristics at the scale of seconds and beyond. Temperature change is assumed to amount to a small perturbation, permitting linear expansion of the temperature-dependent transfer function: (5)where the index i counts the temperature sensors, denotes the output of the i-th sensor, its output in the cold state, and the transfer function in the cold state. For the discrete experimental data, Equation (5) assumes the matrix-vector form: (6)in which the transfer functions and are vectors with rows corresponding to discretized frequency and the temperature dependence is described by the matrix , with the number of sensors. is the vector of the temperature excursions: Given and , this model permits calculation of the temperature-dependent transfer function by Equation (6). Calibration of this model amounts to determining and . This was done by fitting to the data obtained with the training sequences described above. To this end, Equation (6) is extended by horizontal stacking of the transfer function and temperature vectors and least-squares solving for and : (7)where the subscripts of and count the temperature states probed for training. Such modeling is performed globally as well as locally in the temperature domain. Each local linear model is centered about a given , for which a prediction is desired, and its training is restricted to temperature states within a sphere of radius R, which will be referred to as the temperature radius. More precisely, a local model of radius R relies on data taken in temperature states that fulfill (8)A global model uses all training data and establishes a globally linear temperature-to-transfer-function relationship equivalent for all . 2.5 Model performance Model performance was assessed in terms of five metrics. First, self-consistency was characterized by the model’s root-mean-square error (RMSE) of reproducing the training data, which is given by the normalized Frobenius norm of the residual: (9)With insufficient diversity of training data, Equation (7) results in overfitting. In this case, the model will exhibit a high degree of self-consistency but perform badly at predicting scenarios not included in the training. As an indicator of this regime, the condition number of the temperature training matrix was calculated, which is large in the case of overfitting and settles upon robust over-determination of the model. Model prediction performance was assessed in terms of the RMSE of predicted transfer functions, calculated in three ways. First, the frequency-resolved RMSE was determined by taking the mean only over time, that is, over temperature states: (10)Second, the error at each time (i.e., in each temperature state) was determined as the RMSE taken over frequency from 0 to 2 kHz: (11)Third, the total RMSE over time and frequency (0–2 kHz) was calculated, resulting in a single number for each set of predictions: (12) 2.6 Model validation

The model was trained by repeated heating of the system with the training sequences: on three different days within 1 week for the y and z gradients and, for a more extensive study, on seven different days within 3 weeks for the x gradient. Self-consistency and conditioning of local models were determined as functions of the temperature radius.

After training, on an additional day, transfer functions measured after heating with the same sequences were predicted based on the temperatures recorded at the same time. Global modeling and local modeling with variable temperature radius were compared in terms of .

To study the relative importance of the different temperature sensor positions, predictions based on single sensors, subsets of sensors, and all sensors were compared.

The sufficiency of the training sequences was assessed in two ways. First, one training sequence at a time was excluded from training and used to test predictions instead. Second, the model was determined using all training sequences and tested by predicting transfer functions obtained after heating with the cross-validation sequences.

2.7 Simulated imaging

To illustrate heating effects and their correction at the image level, EPI and spiral imaging were simulated based on measured and modeled transfer functions. Assuming a FOV of 23 cm, single-shot EPI was simulated with readouts in the x direction and two readout frequencies, one matching a mechanical resonance just below 1 kHz (EPI on-res, resolution = 2.1 mm) and one off any mechanical resonance (EPI off-res, resolution = 1.5 mm). Single-shot center-out spiral imaging was simulated with resolution = 1.5 mm.

Based on the nominal k-space trajectories of these readouts, actual trajectories were obtained by distortion according to Equations (1) and (2), using three kinds of transfer functions: (i) those measured in a warm state (with temperature readings up to C) and a hot state (up to C), (ii) model estimates for these states, and (iii) the transfer functions measured in the cold state. Raw image data from a digital 2D brain phantom were generated by sampling its Fourier transform along the trajectories measured in the heated system (i). Image reconstruction was then performed with the model-based (ii) and cold-state trajectories (iii).

3 RESULTS 3.1 Temperature patterns and sensor positions

The temperature patterns obtained by infrared imaging are displayed in Figure 2, showing the respective hot coil at a moment when a clear pattern is seen. As is to be expected, the thermal pattern depends on the gradient direction as well as the frequency of the heating sequence. The heating patterns of the x and y gradients are similar up to rotation. The temperature pattern of the x coil is most pronounced because it is the innermost coil. The weak heating pattern of the z gradient reflects active cooling and the largest distance from the surface.

Infrared images of the inner surface of the gradient tube, showing temperature patterns generated by DC and EPI heating. Based on time series of such images, the positions of 13 temperature sensors were optimized for capturing temperature dynamics (green circles)

Based on the infrared recordings, 13 sensor positions (green circles in Figure 2) were found by the optimization procedure described in the Methods section. Temperature monitoring with these sensors (1 to 13) and the six remaining sensors recording cooling-system (14–18) and room temperature (19) is shown at the bottom of Figure 6.

3.2 Effects of heating on transfer functions

Figure 3 shows an example of thermal change in the transfer functions of the x, y, and z gradients, caused by heating with medium-frequency EPI sequences with corresponding readout direction. Two chief effects are observed. First, as the gradients heat up, the magnitudes of the transfer functions shift upwards, along with shifts in phase. Second, the frequencies of mechanical resonances shift down by amounts in the same range as the linewidths.

Examples of transfer functions observed in a hot state and the cold state. The chief effects of heating are a general increase and phase shift of the field response and downshift of mechanical resonances

3.3 Model self-consistency and conditioning

Figure 4 shows the conditioning, self-consistency, and prediction error of the thermal model as a function of temperature radius. With training as used in this study, the model is seen to be ill-conditioned for radii smaller than C. In this case, few of the temperature states visited during training lie within the temperature radius so that the matrix contains critically few columns, underdetermining the model. The conditioning settles beyond C, indicating robust over-determination. Even in the over-determined regime, the condition number per se is still large. This is to be expected because temperature is naturally correlated between sensor positions, resulting in large first singular values of the training matrix. Under-determination at low temperature radii is also reflected by exceeding self-consistency (low ) and large prediction error. settles at 1.5–2.0 per mil for all gradient axes and the prediction errors settle at only somewhat larger values. Stable self-consistency and prediction even at the largest radii indicate validity of the global linear model.

Model conditioning, self-consistency () and of prediction as functions of the temperature radius of the model. At small radius, ill-conditioning and high self-consistency indicate over-fitting. Stable criteria beyond temperature radius of C reflect robust modeling at prediction errors of 0.2–0.3 per mil. Stable criteria also at the largest temperature radii indicate validity of the global linear model 3.4 Global model

Figure 5 shows the performance of the global model in terms of prediction error as a function of frequency (), compared with the deviation from the transfer function in the cold state. The model accuracy reaches almost the noise level with only small residual error predicting the behavior of mechanical resonances.

The prediction performance of the global model is virtually independent of temperature as shown by the same criteria as functions of time in Figure 6. Error incurred by a reference taken in the cold state closely relates to the underlying temperature dynamics shown at the bottom of the figure. It converges back toward zero as the system returns to the cold state. The cool-down and thus the convergence are slow, however, with multiple time constants up to the order of 10 min.

The figure also illustrates that heating changes the behavior not only of the gradient coil used for heating but of all coils, as may be expected due to thermal and mechanical coupling in the same structure. Nonetheless, heating by AC operation with the EPI sequences does affect the actually driven coil the most, unlike with DC heating. Notably, both DC and EPI heating via the actively cooled z coil alter the transfer functions substantially although the observed temperature changes remain moderate.

Performance of the globally linear model in terms of prediction . Prediction accuracy reaches almost noise level. Subtle residual error remains only at the sites of mechanical resonances Prediction error as a function of time during temperature change. Top: Prediction of the global model and the cold reference. Bottom: corresponding sensor temperatures 3.5 Local models

Figure 7 illustrates the transition from the globally linear model to locally linear models with temperature radii ranging from C down to C. As the plots show, the prediction performances of global and local modeling are broadly comparable. The local approach is somewhat superior at modeling thermal change of the mechanical resonances, suggesting that their largest shifts may breach the limits of global linear expansion. At C radius, the local model starts to fail due to underdetermination as observed in Figure 4. In this data, the best radius is at about C. However, in view of the moderate benefits of local modeling, only global modeling was used in the remainder of this study.

Comparison of global and local modeling in terms of prediction . The two approaches perform broadly similarly. However, local modeling at the best temperature radius of C is somewhat better at capturing changes in mechanical resonance effects 3.6 Choice of temperature sensors

Figure 8 reports the dependence of model performance on the underlying choice of temperature sensors. Use of all 19 sensors deployed in this study achieves the best predictions. However, notably, even a single well-placed sensor such as the most predictive in-bore sensor (#7) or one of the outflow sensors on the cooling circuit yield useful predictions, reducing error to the range of 20%–30% relative to the cold reference. Use of more sensors brings down the error significantly further and is helpful particularly for modeling mechanical resonances. Room temperature monitoring (not shown separately) added only negligible predictive value.

Dependence of prediction error () on the placement and selection of temperature sensors, using the global model. The best results are achieved with all 19 sensors. However, even single-sensor monitoring in-bore or of outflowing cooling water reduces prediction error significantly relative to the cold reference (Figure 5) 3.7 Prediction of untrained scenarios

Figure 9 illustrates the necessity of training with different heating scenarios. The model performance deteriorates when training the model with only five of the six training sequences and using it to predict transfer functions observed after heating with the skipped sequence (Figure 9, top). However, upon use of all training sequences, that is, with both AC and DC heating via each of the gradient chains, thermal change was predicted rather well for the validation sequences with AC heating at lower or higher frequency than during training (Figure 9, bottom).

Top: Prediction error () of x gradient transfer functions with incomplete training: One of the training sequences was excluded from model fitting and transfer functions after heating with the excluded sequence were predicted. For z gradient heating, the prediction still worked relatively well, but it failed for the others. Bottom: When relying on all training sequences, the heating effects of the validation sequences with untrained EPI readout frequencies were well predicted Simulated imaging in a warm state (C, left) and a hot state (C, middle), with three single-shot readouts (from the top): EPI off any mechanical resonance, EPI on a resonance just below 1 kHz, and spiral. The shown images were obtained by reconstruction based on cold-state transfer functions (“cold H”) and based on the thermal model (“model H”). Right: Effective trajectories in the cold state (blue), in the hot state (black), and according to the model (red), and their differences hot-cold (“error cold,” blue) and hot-model (“error model,” red). Bottom: spectra of the readout waveforms (x gradient), superimposed with error in the transfer function obtained with the model (black) vs the cold reference (grey) 3.8 Simulated imaging

Figure 10 shows the results of the imaging simulations along with the underlying trajectories. In the EPI trajectories, the strongest heating effect is oscillatory deviation in the readout direction, which is caused by change in both the transfer amplitude and the delay at the readout frequency. On-resonance, large change in delay due to shifting of the resonance caused greater net trajectory differences, up to 0.2 times the sampling interval, although the maximum excursion of the trajectory deviates less from nominal than off-resonance. Reconstructions based on cold-state trajectories exhibit blurring and ghosting. The ghosting is stronger in the on-re