Miriam Löcke

Jeremías Garay

Pamela Franco

Sergio Uribe

Cristóbal Bertoglio

A comparison of phase unwrapping methods in velocity-encoded MRI for aortic flows

Purpose: The phase of a MRI signal is used to encode the velocity of blood flow. Unfortunately, phase unwrapping artefacts may appear when aiming to improve the velocity-to-noise ratio (VNR) of the measured velocity field. This study aims to compare various unwrapping algorithms at hand of ground-truth synthetic data generated using CFD simulations. Methods: We compare four different phase unwrapping algorithms on two different synthetic datasets of 4D MRI flow as well as 26 datasets of 2D PC-MRI acquisitions including the ascending and descending aorta. The synthetic datasets are constructed using CFD simulations of an aorta with a coarctation, with different levels of spatiotemporal resolutions and noise. The error of the unwrapped images was assessed by comparison against the ground truth velocity field in the synthetic data and dual-VENC reconstructions in the in-vivo data. Results: Using the unwrapping algorithms, we were able to remove aliased voxels in the data almost entirely, reducing the L2-error of the data compared to the ground truth by 50%-80%. Results indicated that the best choice of algorithm depend on the spatio-temporal resolution and noise level of the dataset. Temporal unwrapping is most successful on datasets with a high temporal and low spatial resolution (δt = 30ms, h = 2.5mm), reducing the L2-error by 70%-85%, while Laplacian unwrapping performs better with a lower temporal resolution or better spatial resolution (δt = 60ms, h = 1.5mm), especially for SNR 12 as opposed to SNR 15, with an error reduction of 55%-85% compared to the 50%-75% achieved by the Temporal method. The differences in performance between the different methods were found to be statistically significant.

Introduction

Phase contrast MRI is an important technique for evaluating vascular hemodynamics [START_REF] Markl | 4D flow MRI[END_REF][START_REF] Markl | Time-resolved three-dimensional phase-contrast MRI[END_REF]. This technique encodes velocity information in the data phase, so all values are restrained to the range (-π, π), where the maximum encoded velocity V enc is mapped to π. If V enc is set lower than the actual velocity, an image error known as "phase wrapping" occurs as the phase value is wrapped back into the range (-π, π).

Conversely, choosing a too high V enc decreases the velocity-to-noise ratio (VNR), thus lowering the quality of the data [START_REF] Conturo | Signal-to-noise in phase angle reconstruction: Dynamic range extension using phase reference offsets[END_REF]. Since the range of measured velocities in MRI can be very large, the selection of the proper V enc may involve repeated scans. This makes it attractive to post-process the acquired velocity data to correct the wraps introduced by using a low V enc by phase unwrapping techinques. This can be performed in a number of different ways by assuming regularity in the spatial [START_REF] Loecher | Phase unwrapping in 4D MR flow with a 4D single-step Laplacian algorithm[END_REF][START_REF] Bhalerao | Unwrapping phase in 3D MR phase contrast angiograms[END_REF][START_REF] Dong | Phase unwrapping with graph cuts optimization and dual decomposition acceleration for 3D high-resolution MRI data[END_REF] or the temporal dimension [START_REF] Untenberger | Spatiotemporal phase unwrapping for real-time phase-contrast flow MRI[END_REF][START_REF] Itoh | Analysis of the phase unwrapping problem[END_REF][START_REF] Xiang | Temporal phase unwrapping for cine velocity imaging[END_REF] or both [START_REF] Loecher | 4 D Gradient Based Phase Unwrapping for PC-MR Flow Data[END_REF][START_REF] Barnhill | Real-time 4D phase unwrapping applied to magnetic resonance elastography[END_REF].

However, so far, there has been no thorough comparison of these different algorithms. For the vast majority of the algorithms, the code is also not publicly available and hence methods cannot be compared with the same data sets. Moreover, ground truth data in subjects or experimental work would involve for instance multi-venc acquisitions -in order to get high VNR perfectly unwrapped images -which require advanced scan protocols not yet available in clinical practice.

Therefore, the aim of this article is therefore to compare four different phase-unwrapping methods using synthetic data generated using CFD simulations where a ground truth is available. We establish the advantages and drawbacks of the various methods and the cases where they are appropriate. Additionally, the code we used for each of the algorithms will be made publicly available. Moreover, we confirm some of the findings on 2D PC-MRI data sets published previously and acquired with multiple V enc values.

The remainder of the paper is structured as follows. In Section 2 the phase wrapping problem is introduced and we present the different unwrapping methods. In Section 3 we describe the details of the synthetic and subject datasets as well as the implementation details of the algorithms. Section 4 contains the results of the unwrapping algorithms for the different datasets, which are then compared and discussed in Section 5. Finally, Section 6 gives some conclusions.

Theory

Problem Statement

In PC-MRI, the velocity information is encoded into the phase of the transverse component of the magnetization such that there is a linear relationship between phase contrast φ and velocity

v φ = v • π V enc (1)
where the parameter V enc is set by the user of the MRI scanner. Since the phase can only be measured in the interval [-π, π), the V enc parameter decides the range in which the velocity data can be encoded. As a result, when the true velocity exceeds the V enc , velocity aliasing occurs, i.e. the estimated velocity will be u true -2n • V enc instead of u true . Here n ∈ Z depends on the relation between u true and V enc . Analogously, the estimated phase will be φ w = φ true -2n • π where the true phase may now exceed the interval [-π, π).

Out of the available methods, we review the three main algorithms solving the 4D phase unwrapping problem: one relying on temporal smoothness, one relying on spatial smoothness for solving a Laplacian equation over the entire domain, and one using both temporal and spatial smoothness to iteratively unwrap the picture pixel-by-pixel by assigning each pixel a probability of being wrapped.

Phase unwrapping methods

Temporal unwrapping

Temporal phase unwrapping has first been introduced in [START_REF] Itoh | Analysis of the phase unwrapping problem[END_REF]. This method is based on the assumption that at each voxel the phase only varies slowly in time or that the temporal resolution is high enough. Inter-voxel dependency is neglected.

Given a time series of N phase maps φ w (x, t j ), j = 1, • • • , N , the differential phase maps are computed as:

D j (x) = φ w (x, t j ) -φ w (x, t j-1 ), j = 2, ..., N. (2)
According to the aforementioned assumption, these differential maps cannot contain any phase wraps of their own. Therefore any value |D j (x)| greater than π has to be the result of a phase wrap occurring in one of the phase maps. To regain the "correct" differential value, it is sufficient to wrap the differential phase maps into the range (-π, π). Once the wrapping-free differential phase maps have been computed, the unwrapped phase maps can be calculated by integrating over the differential maps, starting at a reference time frame. This requires the existence of a timeframe that does not have any wrapped pixels, which will have to be manually selected. When trying to measure blood velocities, the cardiac cycle contains times with very low velocities, e.g. in diastole, so such a reference timeframe should exist when the temporal resolution is high enough [START_REF] Xiang | Temporal phase unwrapping for cine velocity imaging[END_REF]. The unwrapped phase φ u with a reference timeframe φ w (x, t ref ) at time t ref can then be described as

φ u (x, t j ) =      φ w (x, t ref ) + j i=ref +1 D i (x) for t j > t ref φ w (x, t ref ) -ref i=j+1 D i (x) for t j < t ref φ w (x, t ref ) for t j = t ref (3)
This algorithm is simple to implement and efficient in terms of computational effort. However, it is requires a sufficiently high temporal resolution as well as a wrap-free reference frame, which has to be manually selected. The method can unwrap multiple wraps (which occur if the V enc is less than 50% of the maximal velocity), as long as they occur in separate timesteps, i.e. a wrapped pixel in one timeframe being wrapped again in a later timeframe. This method is sensitive to noise and may falsely unwrap pixels due to noise. Movement of the aorta due to the contraction of the heart also negatively impacts the results for this method as it breaks the continuity of the pixels between the timeframes.

4D gradient-based unwrapping

In gradient-based unwrapping methods [START_REF] Loecher | 4 D Gradient Based Phase Unwrapping for PC-MR Flow Data[END_REF] each pixel is assigned a probability that it is wrapped based on the gradients with the neighbouring pixels, and all pixels with a probability higher than a certain threshold are unwrapped. This method evaluates the data pixel-bypixel and often requires several consecutive iterations. The method assumes smoothness in both space and time. The probability that a pixel is wrapped is estimated by the gradients of the neighboring pixels to the evaluated pixel via

P (x, t j ) = i∈N (x) w s (φ w (i, t j ) -φ w (x, t j )) + t i ∈N (t j ) w t (φ w (x, t i ) -φ w (x, t j ) 2π • (#N (x)w s + #N (t j )w t ) (4)
where N (x) is the set of the spatial neighbours of the pixel x and N (t j ) is the set of the temporal neighbours of the time position t j . w s , w t are weighting factors. It is recommended to set w s to 1 for the spatial dimension and w t to 2.5 for the temporal dimension. Each step is detailed in Algorithm 1.

It is suggested to use a dual threshold approach, where the algorithm first completes a number of iterations with a low threshold in order to break up larger regions of wrapped, and then uses a higher threshold for further iterations to accurately unwrap the remaining pixels. The suggestion for the optimal thresholds are r low = 0.32 and r high = 0.75, for 50 iterations each. This algorithm is very simple to implement and requires little computational effort per iteration. However, depending on the number of iterations, the time required may increase above that of other methods. This method can also only unwrap single wraps and therefore requires that V enc > 50% of the maximal velocity to unwrap successfully.

Algorithm 1: 4D gradient based unwrapping algorithm

Given a wrapped phase image φ w , a set of thresholds r low and r high , and iteration numbers n low , n high :

Set φ = φ w Repeat n low times: foreach (x, j) ∈ Ω × {1, • • • , N } do compute P (x, t j ) if P (x, t j ) > r low then φ(x, t j ) = φ(x, t j ) + 2π if P (x, t j ) < -r low then φ(x, t j ) = φ(x, t j ) -2π Repeat again n high times with r high instead of r low Return φ

Laplacian unwrapping

Laplacian unwrapping is a popular single-step technique that can be used with 2D, 3D, and 4D Flow MRI data. Starting from the relationship

φ(x, t) = φ w (x, t) + 2πn(x, t) (5)
the method assumes the continuity and differentiability of φ and φ w and aims to first find the Laplacian of the true phase. Specifically, the method is obtained by first showing the relation: ∆e iφ(x,t) = -e iφ(x,t) ∇φ(x, t) 2 2 + ie iφ(x,t) ∆φ(x, t), with ∆ and ∇ the spatio(-temporal) Laplacian and gradient operators, respectively. By combining it with Equation ( 5), the following equation is obtained:

∆φ(x, t) = Im(e -iφw(x,t) ∆e iφw(x,t) ). ( 6
)
Simplifying the right hand side leads to

∆φ(x, t) = cosφ w (x, t)∆sinφ w (x, t) -sinφ w (x, t)∆cosφ w (x, t) (7)
Therefore, the unwrapped phase is computed by solving the problem:

φ u (x, t) = ∆-1 (cosφ w (x, t)∆sinφ w (x, t) -sinφ w (x, t)∆cosφ w (x, t))
where ∆ denotes the Laplacian operator but applying some set of boundary conditions. In the case of a spatio-temporal Laplacian operator for 4D datasets, it's necessary to introduce a scaling constant for the temporal dimension. In practice, this scaling is quite robust and was simply set to 1s 2 /cm 2 , as it was in [START_REF] Loecher | Phase unwrapping in 4D MR flow with a 4D single-step Laplacian algorithm[END_REF]. Note that, in discrete data, the Laplacian operation on the right-hand-side's terms is approximated by the sum of the second derivatives along each dimension.

The results differ depending on the chosen boundary conditions -Periodic, Dirichlet or Neumann. In this work, periodic boundary conditions were used on all boundaries, as the implementation is provided by [START_REF] Loecher | Phase unwrapping in 4D MR flow with a 4D single-step Laplacian algorithm[END_REF].

The Laplacian method can be computationally expensive, depending on the method used to apply the inverse Laplacian and the size of the input data, but it can unwrap nested wraps and it is reported to handle noise well since it provides a simultaneously smoothing field.

Methods

Synthetic data

A computational fluid dynamics simulation in physiological regimes was generad as detailed in Appendix A. Let u true be the velocity field generated by this simulation. Then the MRI magnetization measurements were modeled in the following way:

M u meas = C exp(i(φ 0 + πu true /V enc )) + u (8)
where φ 0 is the background phase, C corresponds to the magnitude, and is a Gaussian measurement noise with a mean of zero. PC-MRI requires an additional measurement to recover an estimate of u true since the background phase is unknown. This can by done by turning off the motion encoding, resulting in the magnetization

M 0 meas = C exp(i(φ 0 )) + u . (9)
Now the velocity can be estimated from the two magnetizations via the equation

u meas = ∠ exp(i(∠M u meas -∠M 0 meas )) π V enc . (10)
These synthetic measurements are interpolated into a tetrahedral box mesh placed around the aorta geometry. We consider two different resolutions for the box mesh: 1.5×1.5×1.5mm 3 and 2.5 × 2.5 × 2.5mm 3 . The reference solution is also undersampled in time with τ 1 = 0.03s or τ 2 = 0.06s, leading to a total of 28 respectively 14 time steps. Two different noise levels of 15dB or 12dB in the complex magnetization was applied, and two V enc s were set to 120% and 60% of the maximum velocity, resulting in 192cm/s and 96cm/s, respectively. Only the velocity in the foot-head direction was acquired.

Next to the 4D Flow aorta measurements, we also generated synthetic 2D "slice" measurements by selecting a single 2D slice perpendicular to the foot-head direction and interpolating the velocity into a rectangular hexahedral mesh with a spatial resolution of 2.5 × 2.5mm 2 . The sub-sampling in time as well as the noise levels remained the same. The chosen V enc s were 153cm/s and 77cm/s corresponding to 120% and 60% of the maximum velocity respectively.

Fifty independent realizations of each noise level per V enc and spatial resolution were made. The results in Section 4 therefore show the averages and standard deviations over these fifty realizations.

Subject's data

We used a database of 2D PC-MRI of 26 subjects originally reported in [START_REF] Carrillo | Optimal Dual-VENC (ODV) Unwrapping in Phase-Contrast MRI[END_REF]. Measurements were acquired in a 1.5T Achieva scanner with a 5 elements cardiac coil. The following scan parameters were used: FA = 15 • , TR = 5.5ms, TE = 3.7ms, matrix size = (320, 232). The acquired data are slices perpendicular to the ascending aorta just above the valsalva sinus, in 2D plus time. The spatial resolution was 1mm × 1mm and the temporal resolution ranged from 35ms to 48ms depending on the patient's heart rate. Two different V enc s were used: 150cm/s and 75cm/s.

The MRI images were segmented using MATLAB so that non-zero entries are only present in the ascending and descending aorta.

To investigate the noise sensitivity, the experiments were repeated with a Gaussian noise with a zero mean and a standard deviation of 0.15π added to the phase of the segmented measurements.

Implementation

All algorithms were implemented in both Python and Matlab. The codes are included as supplementary material, including the input CFD data and scripts generating the results.

For the Temporal method, the reference timeframe was selected at diastole, where the flow velocity is low enough to ensure to be aliasing-free. For the synthetic data, this frame was the first time frame. For the volunteer data, the last time frame was used instead. The selection of the reference timeframe does not change the results if it is aliasing-free.

For the Probability method, we used the threshold values suggested in [START_REF] Loecher | 4 D Gradient Based Phase Unwrapping for PC-MR Flow Data[END_REF], i.e. t low = 0.32 and t high = 0.75. Their work suggests using 50 iterations of either threshold, but we adjusted this value empirically to 10 iterations of either threshold.

The implementation of the Laplacian method is a translation of the code provided with [START_REF] Loecher | Phase unwrapping in 4D MR flow with a 4D single-step Laplacian algorithm[END_REF] from MATLAB to Python. This provides both an implementation in 3D as well as 4D.

Error assessment

We compare the unwrapping methods on the error norm

e = u unwrapped -u ref u ref (11)
where u ref is a reference solution. For the synthetic data, u ref is the ground truth, i.e. the result of the CFD simulation without noise or aliasing. For the in-vivo data, the reference solution is gained by unwrapping the data with the multi-V enc algorithm OMME [START_REF] Herthum | Multiple motion encoding in Phase-Contrast MRI: A general theory and application to elastography imaging[END_REF], which is an extension of the dual-V enc unwrapping algorithm presented in [START_REF] Carrillo | Optimal Dual-VENC (ODV) Unwrapping in Phase-Contrast MRI[END_REF], with V enc1 = 150cm/s and V enc2 = 75cm/s.

To assess the statistical significance of the difference between the unwrapping results, a Friedman test with a post-hoc Wilcoxon Signed-Rank test was used, with p < 0.05 being considered statistically significant after a Bonferroni correction to account for multiple comparisons. The statistical analysis was performed in Python.

Results

Computing times

The processing times of the studied algorithms, running on a notebook with 16GB RAM and an Intel Core i5 processor, are shown in Table 1. and4 show the errors and standard deviation for each unwrapping method for all considered spatial and temporal resolutions for the Aorta mesh and the slice mesh respectively. Also included are the errors for the wrapped data and unwrapped with OMME [START_REF] Herthum | Multiple motion encoding in Phase-Contrast MRI: A general theory and application to elastography imaging[END_REF] using V enc1 = 120% and V enc2 = 60% of V max . The reference solution used to compute the error is the ground truth of the fluid simulation without noise.

Algorithm Average time

In Figure 2(a), it can be seen that all unwrapping methods lead to an improvement in the error over the original data. The Temporal method is very close to the baseline OMME solution, whereas the other methods perform worse. The 4D gradient method has a significantly higher error with a larger standard deviation than the other methods. The pattern persists with the increase in noise for SNR 12.

Using the finer grid in 2(b), all methods except 4D gradient perform significantly better, with the Spatiotemporal Laplacian now achieveing similar errors to the Temporal method for both noise levels. Additionally, the standard deviation of the errors has decreased.

For SNR12, the differences between the Temporal results and the OMME results were not statistically significant (p = 0.5 for h 1 = 2.5mm and p = 0.11 for h 3 = 2.5mm), indicating that the Temporal method performs at the same level as OMME with regards to both means and variance in these cases. All other methods were pairwise statistically significantly different (p < 0.05).

For the lower temporal resolution of δt = 60ms in 2(c), the Temporal method has a significantly higher error and larger standard deviation as compared to δt = 30ms. Also, the method's sensitivity to noise appears to have increased. The other methods are not impacted significantly by the change in temporal resolution. As a result, the 4D Laplacian now has similar errors as the Temporal method and lower ones for SNR 12. This pattern persists in 2(d) with the higher spatial resolution. This improves the error values for all methods but the Temporal method, leading to the 4D Laplacian having a significantly lower error and standard deviation. For SNR 12, the 3D Laplacian outperforms the Temporal method as well.

For the lower temporal resolution, the differences in the performance of the methods were in each case found to be statistically significant for each pairwise combination of methods.

Figure 3 shows the results of the unwrapping methods on the Aorta mesh on an example 2D slice including the ascending and descending aorta. It can be seen that the Temporal method and both Laplacian methods unwrap perfectly or nearly perfectly, whereas the 4D Gradient only leads to a very partial unwrap.

2D data

The results on the 2D slices with hexagonal elements are depicted in Figure 4. Here, all methods except for the 4D gradient method perfom very similarly for δt = 30ms. The Spatiotemporal Laplacian has a much lower standard deviation than the other methods, however, especially for SNR 12. For the coarser temporal resolution, the Temporal method does significantly worse than the Laplacian methods, with an error twice as large as that of the Spatiotemporal Laplacian and a much larger standard deviation. The 4D gradient method does not lead to a significant reduction in the error. The only methods without statistically significant differences were OMME and Temporal for δt = 30ms (p = 1 and p = 0.17, respectively) as well as OMME and 4D Laplacian for SNR15, δt = 30ms (p = 0.28), Temporal and 4D Laplacian also for SNR15, δt = 30ms (p = 1), and Temporal and 3D Laplacian for SNR 15, δt = 60ms (p = 0.13).

Examples for the unwrapping results on the 2D slice data can be seen in Figure 5, with a temporal resolution of δt = 30ms and a noise level of SNR 12. It can again be observed that the Temporal method unwraps most wrapped voxels, but also generates a falsely wrapped voxel with a very high velocity. The Spatial Laplacian leads to a good unwrap with only a few wrapped pixels remaining, while the Spatiotemporal Laplacian unwraps perfectly. However, the 4D gradient method only leads to a partial unwrap, with many aliased voxels remaining.

Subject's data

The errors of the unwrapping methods on each of the volunteer MRI acquisitions are depicted in Figure 6. The reference solution used here is acquired by applying OMME to the MRI It can be seen that the reduction in error is lesser than for the synthetic data Still the error is on average reduced by unwrapping by more than 50%. The Temporal method and Spatial Laplacian method perform similarly on this dataset, while the Spatiotemporal Laplacian has the lowest error and standard deviation of all methods. The 4D gradient method, as with the synthetic data, only marginally reduces the error. For the errors on the noised data in Figure 6 show that the errors for all methods increased slightly, but the comparative performance of the methods remained the same.

(a) h 1 = 2.5mm, δt = 30ms (b) h 2 = 1.5mm, δt = 30ms (c) h 1 = 2.5mm, δt = 60ms (d) h 2 = 1.5mm, δt = 60ms
For this type of data, the differences between Temporal and 3D Laplacian were not statistically significant for both the original and noised data (p = 0.79 and p = 0.35, respectively). Additionally, for the noised data, there were no statistically significant differences between the input data and the results of the 4D Gradient method (p = 0.88), confirming that no effective unwrapping was done. All other methods were statistically significantly different from both the input data and each other.

Figure 7 shows example pictures of the unwrapping results from volunteer 7. The Temporal method unwraps very well except for several voxels on the edges of the ascending and descending aorta, where the movement of the aorta disrupts the temporal continuity of the voxels between one timestep and the next. The spatial Laplacian method fails to unwrap a small area on the edge of the ascending aorta, but unwraps everything else successfully .The spatiotemporal Laplacian methods removes all aliased voxels. In contrast, the 4D gradient method removes only a subset of aliased voxels, leading to a partial, fractured unwrapping.

The results of the data with additional noise are shown in Figure 8. It can be seen that all the methods except 4D Gradient deliver similar results to the original data except for a few additional falsely wrapped voxels, while the unwrapping results of 4D Gradient are significantly worse, with no effective unwrapping done.

Discussion

The results on the different datasets are consistent with each other in regards to the relative performance of the different unwrapping methods. This shows that the performance of the uwrapping algorithms on the synthetic data used here is likely a good indicator of their performance on in-vivo data.

On each of the datasets, the Temporal method and the Spatiotemporal method generally outperform the Spatial Laplacian and the 4D gradient method. The 4D gradient method fails to perfom adequately, as the resulting images are still heavily aliased, and is additionally more computationally expensive than the other methods due to the number of iterations. The Spatial Laplacian offers no advantage over the Spatiotemporal Laplacian and performs worse than the Temporal method in almost all cases as well, with the exception of the PC-MRI acquisitions, where there is no statistically significant difference between the Spatial Laplacian and the Temporal method. It is only slightly less computationally expensive than using the higher-dimensional Laplacian and significantly more expensive than the Temporal method.

With a small timestep, the Temporal method achieves equal or better results than the Spatiotemporal method, especially on a coarser mesh. It is also significantly cheaper to compute than the 4D Laplacian. However, it is more sensitive to noise, and has a much larger standard deviation, especially on the 2D datasets. As a result, the Spatiotemporal Laplacian outperforms the Temporal method on datasets with a finer grid, a larger timestep, or more noise. It is overall the most reliable method and often generates completely unaliased images.

An interesting aspect is that the Spatiotemporal method does not drastically outperform the Temporal method on the full 4D Aorta datasets, despite having the advantage of using an additional spatial dimension as compared to the PC-MRI slice datasets.

A limitation of this study is that we have not covered here the performance of the algorithms of data with a lower V enc (< 50% of the maximum velocity). That scenario is however, relatively unrealistic in clinical practice.Moreover, the algorithms were also not tested on full 4D Flow in-vivo MRI acquisitions, nevertheless the synthetic data presented here is realistic in terms of spatio-temporal resolution and noise.

Conclusion

In this work we compared four different algorithms to correct the aliased data and retrieve correct velocity values. Using two different synthetic datasets and in-vivo data, we have shown the advantages and disadvantages of the algorithms and have come to the conclusion that the Spatiotemporal Laplacian described in [START_REF] Loecher | Phase unwrapping in 4D MR flow with a 4D single-step Laplacian algorithm[END_REF] outperforms the others in terms of performance and reliability with acceptable computational cost. The Temporal method is computationally cheaper, but more sensitive to noise and low temporal resolution, and is therefore only preferable if the temporal resolution is high compared to the spatial resolution or computational capacities are highly limited. The other two algorithms (Spatial Laplacian and 4D gradient) offer no advantages compared to these two.

for the velocity u(x, t) and the pressure p(x, t) and with the density ρ = 1.2gr • cm 2 , the dynamic viscosity µ = 0.035P , and with P (t) given by the three-element Windkessel model :

             P = R p, Q + π Q = Γ u • n C d, dπ dt + π R d, = Q . ( 14
)
Here R p, and R d, represent the resistance of the vasculature proximal and distal to Γ , respectively, while C d, models the compliance of the distal vessels. n is the exterior normal vector of ∂Γ, and the initial conditions u(x, 0) = 0, π 1 (0), . . . , π 4 (0) = 85mmHg are given. The discretized model is solved using a modified semi-implicit Chorin-Temam scheme with a backwards Euler scheme for the time derivative. The boundary condition at the inlet is a pulsatile flow described by

u inlet = -U f (t) n, ( 15
)
where U = 75cm/s is a constant amplitude and f (t) is defined as Table 2: Numerical values of the three-element Windkessel parameters for every outlet

f (t) = sin( πt T if t ≤ T Γ 1 Γ 2 Γ 3 Γ 4 R p (dyn • s • cm -5 )
The reference solution for the velocity was established by solving the described model on a mesh consisting of 2,752,064 tetrahedrons and 510,755 vertices, resulting in an average element diameter of 1.1mm. The time step is set as τ = 0.001.

  (a) CFD velocity solution at t = 0.24s. (b) 3D aorta and slice mesh used for the measurements

Figure 1 :

 1 Figure 1: Reference solution of the blood flow from the CFD simulation and meshes used for simulating measurements

Figure 2 :

 2 Figure 2: Error by each method for different mesh sizes and δt on the Aorta mesh

Figure 3 :

 3 Figure 3: Examples of the unwrapped velocity on the Aorta dataset, sliced at z = 5.7cm, with SNR 12, δt = 60ms, at t = 0.18, on the coarse mesh with spatial resolution h 1 = 2.5mm.

  (a) δt = 30ms (b) δt = 60ms

Figure 4 :

 4 Figure 4: Error by each method for different δt on a slice taken at z = 5.7cm.

Figure 5 :

 5 Figure 5: Examples of the unwrapped velocity on a slice taken at z = 5.7cm, with SNR 12, δt = 30ms, at t = 0.18

Figure 6 :

 6 Figure 6: Average error and standard deviation by each method on the 2D PC-MRI data from 26 different subjects.

Figure 7 :

 7 Figure 7: Examples of the unwrapped velocity on MRI data acquisition from volunteer 7 at t = 4

Figure 8 :

 8 Figure 8: Examples of the unwrapped velocity on MRI data acquisition with added noise from volunteer 7 at t = 4

Figure 9 :

 9 Figure 9: 3D aortic model geometry with boundaries

  dyn • s • cm -5 ) 7200 11520 11520 4800 C (dyn -1 • cm 5 ) 4 • 10 -4 3 • 10 -4 3 • 10 -4 4 • 10 -4

Table 1 :

 1 Processing times of the algorithms for one realization

		Aorta (h = 2.5mm) Aorta (h = 1.5mm) 2D slice PC-MRI
	Temporal	0.131s	0.247s	0.019s	0.011s
	Spatial Laplacian	0.605s	1.344s	0.076s	0.123s
	Spatiotemporal Laplacian	0.933s	1.966s	0.074s	0.247s
	4D Gradient	81.915s	163.001s	8.081s	0.307s
	4.2 Results
	4.3 Synthetic data
	4.3.1 3D data
	Figures 2

t T (t -T )e -k(t-T ) if T c > t > T (16)with T = 0.36s being the opening time of the valve, T c = 0.80s the duration of the cardiac cycle, and 1/k = 1/70s the typical time for the closing of the aortic valve. The parameters of the Windkessel boundary conditions are given in Table2.

Acknowledgements C.B. and M.L. acknowledge the funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 852544 -CardioZoom).

A Reference solution setup

The synthetic flow data is generated on the domain Ω, which is the lumen of the aortic geometry shown in Figure 9. The boundary ∂Ω of the geometry is partitioned as

where Γ in is the inflow boundary, Γ w the aortic wall, and Γ the outlet boundaries. On this domain we consider the incompressible Navier-Stokes equations


Y Wang

C Toumoulin

H Z Shu

Z D Zhou

J L Coatrieux

Vessel extraction in coronary X-ray Angiography

Keywords: coronary angiography, centerlines and contours extraction, 2D geometric moments, structure tracking

This paper describes a method to extract the vascular centerlines and contours in coronary angiography. The proposed approach associates geometric moments for the estimation of a "cylinder-like model" and relies on a tracking process. The orientation of the cylinder axis and its local diameter are computed from the analytical expressions of the geometric moments of up to order 2. Experimental results are presented on several images of two sequences that show the efficiency of the method.

I. INTRODUCTION

Coronary network reconstruction from X-ray biplane angiography, has motivated a lot of research both for diagnosis and therapeutic purposes [START_REF] Chen | Kinematic and Deformation Analysis of 4D Coronary Arterial Trees Reconstructed from Cine Angiograms[END_REF][2][3] [START_REF] Sarwal | Three-dimensional reconstruction of coronary arteries from two views[END_REF]. The objective is to reach a quantitative description of coronary geometry and motion, which are required to evaluate the degree of severity of pathology or in the frame of a catheter-based interventional procedure, to help the guidance of the catheter and follow the revisualization process. Because structures are displayed in a 2D format, geometric features and curvature dynamic variations of these structures are difficult to handle. Therefore the 3D coronary reconstruction, which is performed from two views, has to rely on preliminary segmentation, branch formation and matching stages.

This paper is focused on the edge lines and centerlines extraction of the coronary tree. Many methods have been developed to extract these features [START_REF] Toumoulin | Vascular network segmentation in subtraction angiograms," a comparative study[END_REF] [6] [START_REF] Coatrieux | Computational vision and structural modeling in cardiac vascular network reconstruction[END_REF]. But the complex appearances of these structures in the image frame (high curvatures, loops, superpositions, crossings of vessels, and so on) make this detection difficult. This explains why several works have used an interactive mode for the selection of sensitive points for active contour control [START_REF] Chen | Kinematic and Deformation Analysis of 4D Coronary Arterial Trees Reconstructed from Cine Angiograms[END_REF] or for a rough piecewise approximation of the segments followed by a dynamic programming path search [START_REF] Mourgues | 3D+T modeling of coronary artery tree from standard non-simultaneous angiograms[END_REF]. The proposed method is fully automatic and relies on the 2-D geometric moments and the vessel modeling by a cylinder to extract its central axis and local diameters (section II). Section III briefly describes the tracking process. Preliminary results are then presented in section IV, on several angiographic images. Section V concludes on the interest of the method and gives some perspectives.

II. LOCAL SHAPE MODELING AND GEOMETRIC

MOMENTS Moments are integral-based features and are therefore robust against noise. Furthermore low order moments have a direct geometric interpretation. Their application provides consequently an efficient way to compute location, orientation and size features of elongated bright structures.

A. Local vessel modeling

The vessel was locally modeled by a cylinder of orientation θ and diameter d (Fig. 1). In the Cartesian coordinate system (E1), θ represented the angle between the axis of the cylinder and the Ox axis. The vessel intensity was denoted I v and the background intensity I b , with I v > I b .

B. Feature extraction and Geometric moments

The (p + q) th order of geometric moments applied on an N x × N y square window, for an image function, f (x, y) is defined as [START_REF] Shling | Multiresolution moment filters: Theory and applications[END_REF]:

m pq = Nx p=1 Ny q=1
x p y q f (x, y) with p, q ∈ Z +

(1) Furthermore, the (p + q) th central geometric moment is defined so as to normalize moment calculation with respect to the image centroid, thus yielding moments invariant to object translation:

µ pq = p q (x -x) p (y -ȳ) q f (x, y) (2)
where (x, ȳ) represents the coordinates of the local centroid inside the square window: Central moments were applied on an isotropic grid using a circular window. Each pixel inside this window was weighted according to its membership of the circle, by considering a sub-pixel decomposition (Fig. 2). The moment computation was therefore reduced to mask convolutions:

x = m 1,0 m 0,0 and ȳ = m 0,1 m 0,0 (3)
µ pq = f (x, y) * w pq (x, y) (4)
where w pq represents the moment filter mask of order (p, q):

w pq (x, y) = (x -x) p (y -ȳ) q A(x, y) (5)
The distance between the window center and its gravity center allows to detect whether the window is centered on a bright structure or not. Therefore the central moments can be expressed in terms of moments m p,q and this gravity center:

µ 2,0 = m 2,0 -m 0,0 x2 (6) µ 0,2 = m 0,2 -m 0,0 ȳ2 µ 1,1 = m 1,1 -m 0,0 xȳ
The local orientation of the structure is given by these three central moments of the second order. In our case, the vessel, being mapped onto an ellipsoid centered at (x 0 , y 0 ), the orientation angle with respect to the x-axis is given by [START_REF] Boldak | [END_REF]:

θ = 1 2 arctan 2µ 1,1 µ 2,0 -µ 0,2 (7)
The local vessel diameter is then expressed in terms of the zero order moment and the vessel and background intensities:

d = 2 m 0,0 -I b πR 2 π(I v -I b ) 1/2 (8)

III. TRACKING PROCESS

The initialization of the tracking was performed from a list of seed points, automatically extracted at a preliminary stage. From one seed point P i , its position was first refined through an iterative process to make it converge towards the central axis. Then the local diameter was estimated using a multi-scale local moment computation to take into account the size variation of the vessel in presence of pathology (stenosis, aneurysm). Afterwards the tracking process was carried out by shifting a circular window towards the next point P i+1 according to the estimated direction at point P i . The incremental displacement between these two points was made adaptive and depended on the vessel size and curvature [START_REF] Boldak | [END_REF] [10].

A. Estimation of the P i+1 Position

The position P i+1 inside the cylinder corresponded to the center of gravity in the circular window. Nevertheless, the orientation estimated at a given position of the window does not guarantee a precise positioning at the next point along the vessel in particular when highly curved vessels are tracked. An iterative centering was therefore performed based on the first order moments to move the point towards the vessel center (Fig. 3). It consisted in:

1. Computing the gravity center inside the window, 2. Moving the center of the window to its gravity center 3. Reiterating step 1 and 2 till the center of the window coincides with its gravity center.

B. Estimation of the vessel and background intensities

The intensities inside (I v ) and outside (I b ) the vessel were made adaptive at each step of the tracking to take into account the heterogeneity of the intravascular and environmental properties. The background intensity I b was set to the mean value computed outside the vessel in its normal direction. In the same way, the vessel intensity I v was equal to the mean intensity computed along the estimated vessel direction.

C. Multi-scale detection strategy

The size of the vessels can widely vary (from 25 to a very few pixels) when considering the overall coronary tree or when diseases are concerned (stenosis, aneurysm). In addition, the analytical diameter formulation as reported in (8) depends directly on the window size as well as Fig. 3. Tracking process and iterative centering of the point P i+1 the precision of the localization of the vessel center P i . Indeed, the vessel corresponds to brighter structures in the image. Therefore optimal parameters are achieved when the spherical window is centered on the vessel axis and its size exactly fits the vessel size. We made this window size adaptive to adjust it, at each step of the tracking, according to the estimated diameter of the vessel. Considering R i-1 , the local estimated radius of the vessel at point P i-1 , the local diameter computation at point P i , was performed at different scales according to the following algorithm:

Algorithm Multi-scale diameter computation R i,n=1 = R i-1 Repeat 1 Estimate P i , θ i , d i , I b , I v 2 If (R i,n < d i /2) R i,n+1 = R i,n + 1 Else if (R i,n + 1.25 > d i /2) R i,n+1 = R i,n -1 Else R i,n+1 = R i,n 3 n = n + 1; while (R i,n = R i,n-1 ) End
The local direction θ i was computed on a window centered at point P i and with radius (R i,n + 1).

D. Re-initialization procedure

The tracking can be stopped when meeting high curvature or superposition of vessels. A square box, twice the size of the computation window, centered on the last extracted point P i , was then open and an intersection of vessels with the perimeter of the window was searched. The previously visited pixels were marked in order to avoid a potential backtrack.

E. End Conditions

The stopping conditions, checked at each step, allow testing the location of the estimated point (inside/outside). This decision relies on the mean intensity and standard deviation values estimated within the current tracking window compared to a reference corresponding to a minimal diameter vessel d min that can be detected. In any case, the vessel intensity must be higher than the background one (I v > I b ).

F. Branch formation

The tracking allowed to extract points located on the central axis of the vessel and to estimate the local diameter at each of these points. The objective was then to establish the continuity between these points (centerlines and edge lines respectively) in order to reconstitute the vascular branches. Centerlines were first reconstituted based on intensity and directional similarity criteria.

A branch was delimited by a bifurcation point or a end point. Then for each centerline, edge points were recovered in the normal direction to the vessel. The spatial location of each point, related to the central axis, was used to construct the left and right contours (Fig. 4). At each bifurcation point, a circle was drawn whose radius was equal to the largest diameter among the estimated radii in the neighborhood of the bifurcation at each starting centerline. Then, edge points were picked up at the level of the intersection of each centerline with the circle. The continuity was then obtained by simply connecting the edge points, located in the same side related to the central axis.

IV. RESULTS

The algorithm was applied on angiographic subtracted images (512 × 512, 8 bits) acquired in standard conditions (Left and Right Anterior Oblique views) with a General Electric Digital Subtraction Angiography system (Fig 5). Each sequence includes 15 images. The images were first convolved with a Gaussian filter of size 9 × 9 and standard deviation equal to 0.5 to smooth the noise. A morphological top hat operator was then applied with a spherical structuring element of size 41 × 41 to enhance the structure of interest. Four horizontal lines were then drawn on each image and points corresponding to maximum intensity points, were extracted on each line. These points were normally located on the crest line of the main vascular branches. They were kept as initial seed point to initialize the tracking process.

The tracking process was applied on each image of each sequence. For illustration, extracted centerlines and contours are shown in Fig. 6 and Fig. 7 for two images of each sequence. Results are superposed on the original images V. CONCLUSION The tracking process allows extracting the centerlines and contours of the main vascular branches within 10 seconds per image on a Pentium-4 PC with 1.5 GHz CPU and 512 Mb RAM. Results show a good efficiency of the method on the set of the images of each sequence. Some refinements are nevertheless necessary to obtain more accurate and smooth contours. The next stages will be devoted to the evaluation of the algorithm to other angiographic sequences: (1) its robustness with regard to the presence of other possible conformations in the images and (2) the accuracy of the detection.

  This work is supported by both the National Natural Science Foundation of China (Grant No.60272045) and the program for New Century Excellent Talents in University.

Fig. 1 .

 1 Fig. 1. Vessel modeling: (a) initial coordinate system (E1); (b) coordinate system after alignment on the y axis (E2)

Fig. 2 .

 2 Fig. 2. 2-D Moment mask: Sub-pixel decomposition and weights A i,j associated to each pixel for a circular window of radius R = 1 (defined in a square window of size N = 3)

Fig. 4 .

 4 Fig. 4. Edge point recovering: for each extracted central axis point (blue point), contour points are got from the local estimated diameter in the normal direction to the vessel (red points).

Fig. 5 .Fig. 6 .Fig. 7 .

 567 Fig. 5. Coronary angiography: Two orthogonal viewpoints acquired from a GE DSA system (Left and Right Anterior Oblique views). Each sequence includes 15 images.


Han Deng

Jingyang Fang
email: jingyang.fang@duke.edu

Yang Qi
email: qiyang@nwpu.edu.cn

Vincent Debusschere
email: vincent.debusschere@grenoble-inp.fr

Yi Tang
email: yitang@ntu.edu.sg

Dynamic Performance Limitation and Enhancement of Grid-Forming Converters

Keywords: Grid-forming converter, inner loop control, coupling effect

Grid-forming converters (GFMCs) are emerging as promising candidates in more-electronics power systems for their capabilities of forming voltage and frequency. Regarding the control of GFMCs, lots of research have been done on the active and reactive power control loop and the voltage control loop of standalone voltage source converters. However, little literature elaborates the design of the voltage and current controllers when GFMCs are connected to the grid. This paper analyzes the stability of standard dual-inner-loop GFMCs. The requirements for the stability are analyzed and the facts that limiting the control bandwidth are identified. Specifically, it is revealed that in the synchronous dq frame, the coupling terms introduced by the grid impedance have a non-negligible impact on system dynamics. More importantly, the conventional dual-loop control scheme inherently limits the control bandwidth of GFMCs to be less than 50 Hz. To improve the speed of tracking the reference, a novel control scheme with a band-pass filter and a band-stop filter is proposed. With the proposed control scheme, the bandwidth of the GFMC is increased without causing instability.

I. INTRODUCTION

The grid-forming converters (GFMCs) are emerging as an attractive solution for grid-support due to their merits of voltage and frequency regulation [START_REF] Lasseter | Grid-forming inverters: A critical asset for the power grid[END_REF]. The basic idea for GFMCs is to embed the characteristics of synchronous machines, e.g. self-synchronization, droop and inertia behaviors, into the power control of voltage source converters (VSCs) [START_REF] Beck | Virtual synchronous machine[END_REF], [START_REF] Qing-Chang | Selfsynchronized synchronverters: Inverters without a dedicated synchronization unit[END_REF]. Compared with the grid feeding converters (GFDCs) whose currents are controlled directly, the GFMCs are independent of phase lock loops (PLLs) and more suitable in weak grids. Furthermore, GFMCs can form grid voltage and frequency even without a reference voltage source [START_REF] Lasseter | Grid-forming inverters: A critical asset for the power grid[END_REF].

The design and control of outer loop power controllers have been widely discussed. [START_REF] Zhong | Synchronverters: Inverters that mimic synchronous generators[END_REF] emulates the synchronous generators behaviors with the VSCs. [START_REF] Wu | Small-signal modeling and parameters design for virtual synchronous generators[END_REF] builds the small signal model of the GFMC and gives the parameter design guideline of the power loop controllers. The equivalence between low-pass-filtered power feedback and virtual inertia in droop controlled GFMC is proved in [START_REF] Suul | Equivalence of virtual synchronous machines and frequency-droops for converter-based microgrids[END_REF]. [START_REF] Liu | Comparison of dynamic characteristics between virtual synchronous generator and droop control in inverterbased distributed generators[END_REF] compares the power loop oscillation damping methods in GFMC controls. A generalized stability requirement for individual GFMCs power controller design regarding the whole interconnection system stability is proposed in [START_REF] Qi | Passivity-based synchronization stability analysis for power-electronic-interfaced distributed generations[END_REF].

Nonetheless, most research works use direct voltage control [START_REF] Zhong | Synchronverters: Inverters that mimic synchronous generators[END_REF] or ignore the inner loop control [START_REF] Wu | Small-signal modeling and parameters design for virtual synchronous generators[END_REF], [START_REF] Qi | Robust power sharing control for parallel three-phase inverters against voltage measurement errors[END_REF], i.e. the voltage control, for the assumption that the bandwidth of voltage control loops is much faster than that of the power loops. In this case, the VSCs are considered as ideal voltage sources with no delay and tracking error in GFMCs [START_REF] Wu | Small-signal modeling and parameters design for virtual synchronous generators[END_REF]. However, the inner loop controllers are important for precise voltage tracking, harmonics filtering and LCL resonance damping [START_REF] Wang | Stability analysis and controller synthesis for single-loop voltage-controlled vsis[END_REF], [START_REF] Qi | Two-dimensional impedanceshaping control with enhanced harmonic power sharing for inverterbased microgrids[END_REF]. Although faster tracking speed is expected from GFMC to omit the voltage dynamics when designing the outer loops, few research works discussed the design and control for GFMCs voltages. The voltage tracking speed of the GFMC, and whether the voltage control loop can be omitted especially in GFMCs with fast dynamics, have not been discussed thoroughly. [START_REF] Wang | Stability analysis and controller synthesis for single-loop voltage-controlled vsis[END_REF] analyzes the stability requirements and gives design guidelines for single loop voltage control of VSCs under isolated operation. Yet, [START_REF] Fang | More-Electronics Power Systems: Power Quality and Stability[END_REF] reports that the LC-filtered isolated model has different characteristics from the grid-tied VSC model. [START_REF] Liao | Evaluation of voltage regulators for dual-loop control of voltage-controlled vscs[END_REF] analyzes the bandwidth of the grid-forming converter controlled by double loop proportional-resonance (PR) controllers with LC-filtered islanded model and proposed some methods of increasing the bandwidth. Nevertheless, it is hard to determine the tracking speed through the bandwidth of the model controlled by PR controllers.

IEEE COPYRIGHT -https://www.doi.org/10.1109/ECCE47101.2021.9595056

v gfdref ω o t v giabc SPWM PWM 1-6 dq abc dq abc v gfabc v gfqref (0) G v (s) G v (s) dq abc dq abc ω o t i giabc G c (s) G c (s) dq abc dq abc ω o t v gid v giq L gi i giabc v gfabc L gg i ggabc v ggabc C gf Grid v gdc + - + - Fig. 1. System schematic diagram.
To further explore the tracking speed of the GFMC and better design its controller, this paper gives the requirements for stabilizing the converter and analyzes the small-signal model of the dual-loop voltage control for grid-forming converter in the synchronous dq frame. This paper also points out that when analyzing the tracking speed for GFMC with an LCL filter in the rotational frame, the coupling effect on the grid side impedance cannot be ignored. It is revealed that under traditional dual-loop proportional-integral (PI) control schemes, the bandwidth is hard to exceed 50 Hz. To increase the tracking speed of the voltage control, a novel control scheme is proposed in this paper.

The remaining of this paper is organized as follows. Section II presents the system model and derives basic stability requirements. Section III analyzes the influence of the coupling effect on the system dynamics. Section IV proposes an enhanced voltage control strategy which can improve the tracking speed. The simulation results are shown in Section V. Section VI concludes the paper.

II. SYSTEM MODEL

Fig. 1 shows the system schematic diagram, where a GFMC connects to the grid through an LCL filter composed of L gi , C gf and L gg , and the ideal voltage sources v ggabc emulate the grid. The grid side inductance L gg includes the filter inductance and the transmission line inductance. The symbols v gfabc denote the voltages over the capacitor C gf , i giabc and i ggabc denote the converter currents and the grid currents, respectively. The voltage and current controllers, i.e. G v (s) and G c (s), mainly regulates the capacitor voltage v gfabc and feedback the converter current i giabc to damp the LCL resonance [START_REF] Tang | Exploring inherent damping characteristic of lcl-filters for three-phase grid-connected voltage source inverters[END_REF]. The voltage and current controllers are implemented in the dq synchronous frame to better get the bandwidth of following a dc reference, where ω o t stands for the phase angle of the dq frame and the subscript dq means the values in the dq frame. The reference values for the voltage controllers, i.e. the phase angle ω o t and the voltage amplitude V gfref or v gfdref , result from the active and reactive power controllers, respectively. As the power controllers are out of the scope of this paper, they are not shown in the schematic diagram.

Fig. 2 illustrates the block diagram of the small-signal voltage control model in the synchronous dq frame. G v (s) and G i (s) represents the voltage and current controller, respectively, which employ the PI controllers and are expressed as

G v (s) = G vpi (s) = K vp + s/K vi , (1)
G i (s) = K cp + s/K ci , (2)
where K vp and K vi , as well as K cp and K ci are the proportional and integral gains for the voltage as well as current controllers, respectively. G d (s) denotes the computation and the pulse width modulation (PWM) delay function

G d (s) = e -1.5sTs , (3)
where T s is the sampling and switching period. Previous research proves that under certain assumptions, the inverter controlled by PI controllers in the dq frame has almost equivalent dynamics with inverters controlled by PR controllers in stationary frames [START_REF] Zou | Stationary frame equivalent model of proportional-integral controller in dq synchronous frame[END_REF], i.e.

G v (s) = G vpr (s) = K vp + K r s/(s 2 + ω 1 2 ), (4)
where K r = 2K vi is the resonance gain, ω 1 = 2πf n is the nominal angular frequency. The open-loop gain of the voltage control model in the stationary frame can be derived by ignoring the coupling terms in Fig. 2 and transforming the signals in the rotational frame into those in the stationary frame, i.e.

G ol v (s) = Gv(s)Gi(s)G d (s)sLgg s 3 LgiC gf Lgg+Gi(s)G d (s)LggC gf s 2 +(Lgi+Lgg)s+Gi(s)G d (s) = Gv(s)Gi(s)G d (s)sLgg LgiC gf Lggs(s 2 +2ςωrs+ωr 2 )+Gi(s)G d (s) , (5)
where

ω r = L gi + L gg L gi C gf L gg , ς = G i (s)G d (s) 2L gi ω r , (6)
whose bode diagram is shown in Fig. 3. Other system parameters and control parameters are shown in Table I and Table II, respectively. As observed in Fig. 3, the stability of this system is affected by the magnitude and the phase at the resonance peak. Here we only discuss one common situation, which requires

|G ol v (jω r )| = | jGv(jωr)ωrLgg -C gf Lggωr 2 +1 | = | -jG v (jω r )ω r L gi | < 1, (7)
as there will be a -180°phase drop around the resonance frequency ω r . A controllable resonance peak requires that the damping coefficient in ( 6) is larger than 0, i.e.

ω 1 C gf ω 1 C gf ω 1 L gi ω 1 L gi G v (s) G i (s) G d (s) 1/sL gi 1/sC gf v gfdref v gfd v gid i gid i gidref 1/sL gg i ggq v ggd G v (s) G i (s) G d (s) 1/
V dc 200 V Re (G i (jω r ) G d (jω r )) > 0.
The requirement can be further translated into

Re(G i (jω r )G d (jω r )) = Re((K cp + K ci /jω r )(cos 1.5T s ω r + j sin 1.5T s ω r )) = K cp cos 1.5T s ω r + (K ci /ω r ) sin 1.5T s ω r = K cp 2 + (K ci /ω r ) 2 cos(1.5T s ω r -atan(K ci /ω r /K cp )) > 0 ⇒ 1.5T s ω r -atan(K ci /ω r /K cp ) < π/2 (8
) which further yields a conservative requirement for the LCL filter parameter design as

f r < f s /6, (9)
for all values of K cp and K ci , where f r = 2πω r , f s = 1/T s .

As there is little difference near the resonance frequency for PI and PR controller when K vi = K r [START_REF] Pan | Capacitor-currentfeedback active damping with reduced computation delay for improving robustness of lcl-type grid-connected inverter[END_REF], the requirement is also remains the same for the PI-controlled converter in the synchronous frame. However, the PR controller is only effective around the nominal frequency. Furthermore, when designing the controller parameters, the bandwidth, or in other words, the tracking speed for the ac reference is hard to achieve from the design process. As a result, to further understand and model them in the GFMC system, the model of voltage control is further built in the synchronous frame as shown in Fig. 2.

III. COUPLING EFFECT ANALYSIS

As seen in Fig. 2, due to the coupling terms, the transfer function of the system is hard to derive. In uninterrupted power supply (UPS) systems or grid-feeding converter systems, it is a common practice to ignore the coupling terms or add decoupling terms when designing the controllers in dq frame [START_REF] Blasko | A novel control to actively damp resonance in input lc filter of a three-phase voltage source converter[END_REF]. However, due to the coupling effect on the additional grid impedance, the coupling terms cannot be ignored anymore. As the transfer function and the open-loop gain of Fig. 2 are hard to derive, we build the state-space model for Fig. 2, whose expression is given in the appendix. The delay function is linearized by the Pade approximation, which yields

G d (s) = e -1.5sTs = (s -0.75T s )/ (s + 0.75T s ) . (10)
The comparison between the open-loop gains with and without the coupling terms is shown in Fig. 4. As observed from Fig. 4, the bode diagrams with and without the coupling terms are different, especially in the low-frequency region. The PI controller is employed to achieve zero steady-state error. However, when ignoring the coupling terms, the open-loop gain model does not show infinite gain when the frequency converges towards 0. Furthermore, near the resonance frequency, the model without the coupling terms only has one resonance peak. In comparison, the model with the coupling terms has two similar peaks around f r ± 50 Hz and one valley around the resonance frequency f r . ( 7) gives an approximate design criterion for the voltage controller as

K cp = 0.1 K cp = 1 K cp = 10 K cp = 0.1 K cp = 1 K cp = 10 Magnitude
|G ol v (jω r )| = | -jG v (jω r )ω r L gi | = |(K vp + K vi /jω r )ω r L gi | = L gi (K vp ω r ) 2 + K vi 2 < 1 . ( 11
)
Although the bandwidth or the tracking speed is highly related to K vp , [START_REF] Qi | Two-dimensional impedanceshaping control with enhanced harmonic power sharing for inverterbased microgrids[END_REF] tells that the stability of the system is much more sensitive to K vp than K vi . As a result, we set K vp as 0. Then [START_REF] Qi | Two-dimensional impedanceshaping control with enhanced harmonic power sharing for inverterbased microgrids[END_REF] will give an approximate requirement for stability as K vi < 1/L gi . Noticeably, there will be a slight amplitude sag around 50 Hz. By trials shown in Fig. 5, we found that the slight sag is highly related to the proportional gain K cp of the current controller G i (s). Also, the sag tries to prevent the bandwidth of the system from exceeding 50 Hz. With the amplitude sag, the maximum bandwidth can only reach around 30 Hz.

IV. ENHANCED VOLTAGE CONTROL

As revealed in Section III, 2 facts are preventing us from increasing the bandwidth of the voltage control. First, the resonance peak value is highly dependent on the voltage controller G v (s), so the gains of G v (s) are restricted. Second, the voltage amplitude sag around 50 Hz is holding back the bandwidth from exceeding 50 Hz. The two problems are addressed separately to improve the voltage control performance. In the first place, a band-stop filter, i.e.

G bs (s) = (s 2 + ω r 2 )/(s 2 + ω r /Q bs s + ω r 2 ), (12)
is employed around the resonance frequency to damp the resonance peak, where Q bs is the quality factor chosen to be 0.4. Second, a bandpass filter, i.e.

G bp (s) = (s 2 + ω 1 /Q bp s + ω 1 2 )/(s 2 + ω 1 2 ), (13)
is utilized to counteract the amplitude sag around 50 Hz, where Q bp is the quality factor selected to be 3. The expression of the designed voltage controller becomes

G v proposed = G bs (s) G bp (s) G v (s) . (14)
Fig. 6 shows the bode diagram before and after applying the proposed voltage controller, where the integral gain K vi increases from 200 to 1800 to increase the bandwidth. As observed from Fig. 6, with the proposed controller, the resonance peak magnitude is well damped below 0 dB, and the magnitude sag around 50 Hz is well compensated. Compared to the system shown in Fig. 4, the bandwidth of the system shown in Fig. 6 is increased from 20 Hz to 200 Hz without causing instability.

V. SIMULATION RESULTS

To verify the analysis and the proposed controller, simulation is conducted under the PLECS environment. Fig. 7 -9 illustrate the simulation results. Fig. 7 verifies the stability requirements deduced in [START_REF] Qi | Two-dimensional impedanceshaping control with enhanced harmonic power sharing for inverterbased microgrids[END_REF], i.e. f r < f s /6. The switching frequency is arranged at 10 kHz, which results in a critical resonance frequency of f s /6 = 1667 Hz. When C gf = 7 µF, the resonance frequency 1553 Hz is smaller than the critical frequency, which results in a stable system. When C gf changes from 7 µF to 6 µF, the resonance frequency becomes 1667 Hz > f s /6, which yields an unstable voltage waveform.

Fig. 8 and Fig. 9 compare the performance with the traditional dual-loop voltage controller and the proposed voltage controller. The parameters are the same as those used in Section V. When there is a step-up voltage reference change at 1s, the conventional dual-loop voltage controller takes more than one cycle to settle because the bandwidth is restricted within 50 Hz. On the contrary, with the proposed enhanced voltage controller, it only takes around 0.005 s to settle, and no stability problem is caused.

VI. CONCLUSIONS

This paper analyzed the stability of dual-loop-controlled grid-forming converters. First, the stability requirement for conventional dual-loop voltage controllers is pointed out, i.e. f r < f s /6. Second, it is revealed that in the synchronous dq frame, the control bandwidth is limited to be less than 50 Hz by the resonance frequency magnitude and a magnitude sag at 50 Hz caused by the coupling effect. Third, an enhanced controller is proposed, which addresses the bandwidth restriction and boosts system dynamics. The effectiveness of the proposed controller is verified in the simulation, showing a great improvement in the tracking speed.

APPENDIX

The state space model of the system block diagram shown in Fig. 2 can be expressed as

ẋ = Ax + Bu, (15)
where x = [x vd x vq x id x iq v gid v giq i gid i giq v gfd v gfq i ggd i ggq ] T is the states and u = [v dref v qref v sd v sq ] T is the system inputs. The expression of A and B is shown in ( 16) and ( 17), respectively.

A =                        0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 K vi 0 0 0 0 0 -1 0 -K vp 0 0 0 0 K vi 0 0 0 0 0 -1 0 -K vp 0 0 8KcpKvi 3T d 0 8Kci 3T d 0 -4 3T d 0 - 8Kcp 3T d 0 - 8KcpKvp 3T d 0 0 0 0 8KcpKvi 3T d 0 8Kci 3T d 0 -4 3T d 0 - 8Kcp 3T d 0 -
1 Lgg -ω 1 0                        (16)
B =                       1 0 0 0 0 1 0 0 K vp 0 0 0 0 K vp 0 0 8KcpKvp 3T d 0 0 0 0 8KcpKvp 3T d 0 0 - KcpKvp Lgi 0 0 0 0 - KcpKvp Lgi 0 0 0 0 0 0 0 0 0 0 0 0 -1 Lgg 0 0 0 0 -1 Lgg                       (17)

f

  = f r = ω r / 2π f = 50 Hz

Fig. 3 .

 3 Fig. 3. Open-loop gain under the stationary frame.

Fig. 4 .

 4 Fig. 4. Open-loop gain with and without the coupling terms under the synchronous frame.

Fig. 5 .

 5 Fig. 5. Open-loop gain under different current proportional gains.

Fig. 6 .

 6 Fig. 6. Open loop gain under original controller and newly designed controller.

C

  gf = 6 μF f r = 1677 Hz C gf = 7 μF f r = 1553 Hz f s = 10000 Hz f s / 6 = 1667 Hz

Fig. 7 .

 7 Fig. 7. Simulation voltage waveform under different resonance frequencies.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Simulation voltage v gfabc waveforms (red,yellow and blue ones) and their amplitude (black) under the traditional dual-loop voltage control.