Nicolas Baghdadi
email: baghdadi@teledetection.fr

Pierre Todoroff

Mehrez Zribi

Multitemporal observations of sugarcane by TerraSAR-X sensor

Keywords: Sugarcane crop, TerraSAR-X, Reunion Island I

The potential of TerraSAR-X (X-band) in monitoring sugarcane growth was investigated on Reunion Island. Multitemporal TerraSAR data acquired at various incidence angles (17°, 31°, 37°, 47°, 58°) and polarizations (HH, HV, VV) were analyzed in order to study the behaviour of SAR (synthetic aperture radar) signal as a function of sugarcane height. The potential of TerraSAR for mapping the sugarcane harvest was also studied. Radar signal increased quickly with crop height until a threshold height, which depended on polarization and incidence angle. Beyond this threshold, the signal increased only slightly, remained constant, or even decreased. TerraSAR data showed that after strong rains the soil contribution for the backscattering of sugarcane fields can be consequent for canes with heights of terminal visible dewlap (htvd) less than 50cm (total cane heights around 155cm). Finally, TerraSAR data at high spatial resolution were shown to be useful for monitoring sugarcane harvest when the fields are of small size or when the cut is spread out in time. The radar incidence of 37° is more suitable to monitor the sugarcane harvest.

INTRODUCTION

Sugarcane is one of the most important crops in the tropics, with a global production estimated at 1,250 million tons a year and a cropped area of about 20 millions hectares. One of the main needs expressed by sugarcane industries is to have information on the harvest progress throughout the harvest season. The dynamic mapping of sugarcane harvest on a large spatial scale allows optimized cutter deployment, transport operations, efficiency of factories, and finally permits a better estimation of the effective yield.

The use of optical images is sometimes limited because of atmospheric conditions and cloud cover. Indeed, the interval between two cloud-free images is sometimes too long (more than 2 months); this makes difficult the discrimination between a standing crop and the regrowth in a field harvested at the beginning of the harvest campaign. On the contrary, Synthetic Aperture Radar (SAR) provides measurements day and night, regardless of meteorological conditions. With their frequent revisits, SAR sensors are very useful remote sensing data sources for agriculture monitoring in tropical regions. The new generation SAR sensors, such as TerraSAR-X, allow the acquisition of images at very high spatial resolution (~1 m).

Moreover, its short revisit interval makes it possible to monitor the harvest with high temporal frequency (daily to weekly).

The SAR sensors currently operational are TerraSAR-X (X-band), COSMO-SkyMed (X-band), PALSAR/ALOS (Lband), ASAR/ENVISAT (C-band), RADARSAT-1/2 (Cband), and ERS-2 (C-band). New generation SAR sensors (TerraSAR-X, PALSAR/ALOS, RADARSAT-2, COSMO-SkyMed) should considerably increase the current capacities for monitoring continental surfaces by providing data more frequently than before and at high spatial resolution (metric for TerraSAR-X). The increasing number of SAR satellites and the high revisit interval of new sensors (TerraSAR-X and COSMO-SkyMed) makes it possible to follow the harvest season with high temporal frequency (daily to weekly).

This study examined the relationship between TerraSAR signal and sugarcane height as a function of instrumental parameters (polarization and incidence), and precipitation. In addition, the potential of TerraSAR-X for mapping harvested sugarcane crop was studied.

II. STUDY AREA AND DAT SET

The study site covers a sugarcane farm located at the south of Reunion Island, close the town of Saint Pierre (latitude: 21°19' S -longitude: 55°31' E; Figure 1). The study site is composed mainly of agricultural fields intended for growing sugarcane. Fifteen sugarcane fields of an average size of 9 ha were studied: {2, 3, 4, 5, 6.1, 6.2, 12. 1, 12.2, 123, 15, 16, 18, 191, 192, 20}. These training fields extend on 4.5km length approximately, between 100m to 500m altitude. TerraSAR-X images were acquired over our study site. The images belong to the KALIDEOS database set up by the CNES (French space Agency) [START_REF] Deboissezon | Reference remote sensing data bases: Temporal series of calibrated and ortho-rectified satellite images for scientific use[END_REF].

64 TerraSAR-X images (X-band ~ 9.65 GHz) were acquired between 14 th of December 2008 and 20 th of January 2010 with a great range of incidence angle (17°, 31°, 37°, 47° and 59°), and in mono-and dual-polarization modes (HH, VV, HH/VV, HH/HV, VH/VV). The imaging modes used were Spotlight and Stripmap. The pixel spacing of TerraSAR images was between 1 and 3m.

Radiometric calibration using MGD (Multi Look Ground Range Detected) TerraSAR images was carried out in using the following equation [START_REF] Fritz | TerraSAR-X Ground Segment Level 1b Product Format Specification[END_REF]:

DN Ks

). It represents the influence of different noise contributions to the SAR signal. The NEBN is described using a polynomial scaled with Ks. The polynomial coefficients are derived from the TerraSAR product file. All TerraSAR images were then georeferenced using GPS points.

The strong values of NEBN found for images acquired in Stripmap mode did not allow a calibration of many pixels because the term Ks.DN² was lower than the noise NEBN. This problem is very frequent for pixels corresponding to smooth areas (specular reflexion), such as harvested fields. Moreover, the results show that the influence of the noise is stronger for cross-polarizations than for co-polarizations because even if the NEBN is of the same order of magnitude for cross-and copolarizations, the term Ks.DN² is weaker for crosspolarizations. Many pixels impossible to calibrate was also observed at high incidences. These aberrant pixels (Ks.DN² < NEBN) were not used in the calculation of the statistics (for certain images, nearly 30% of aberrant pixels were found what represents a strong loss of information).

In practice, the mean backscattering coefficients were calculated from calibrated SAR images by averaging the linear σ° values of all pixels within reference fields or (sub-fields in the case where only a part of field is harvested).

Ground truth measurements of sugarcane height were performed on several reference fields from November 07, 2008 to June 06, 2009. On each reference field, two experimental areas of 1.5m x 1.5m were used to collect the sugarcane height, number of stems and leaves. Ground measurements showed that the sugarcane in our study site grows about 25cm per month during the five first months, 40cm between the 6 th and 9 th month, and then of about 10-20 cm per month until reaching the mature height of the cane. The ground measurements of the sugarcane height correspond to the height of terminal visible dewlap (htvd). They exclude the leafy tops which have heights of the order of 55cm for sugarcane with htvd of 20cm, 105cm for htvd of 50cm, and of 125cm for htvd between 100 and 180cm. Beyond htvd of 180cm, the leafy top height is about 135cm. For our reference fields, the mean number of stems and leaves was about 17 and 77 per m², respectively (with a standard deviation of about 7 and 30, respectively).

In addition, the farmer of our reference fields also provided the harvesting dates of each reference field. Daily precipitation data recorded at four meteorological stations located on the farm were also used: Bérive-2, Isautier-Bérive, Isautier-Foyer, and Isautier-Ringuin. The effect of soil moisture content was taken into account in this study using precipitation data. Indeed, soil moisture measurements were difficult to carry out because the terrain is inaccessible in rainy weather and the soil is covered with mulch (dead leaves).

III. RESULTS

Sensitivity of TerraSAR-X signal to sugarcane height

The sensitivity of TerraSAR-X signals has been analyzed as a function of sugarcane height (htvd). Results show that the radar signal increases with the sugarcane height for the fields at the beginning of growth (htvd and total cane height respectively lower than 50 cm and 155 cm, depending on incidence angle and polarization) (Figure 2) [START_REF] Baghdadi | Multitemporal observations of sugarcane by TerraSAR-X images[END_REF]. The growth of the sugarcane leads to increase of its height, number and size of leafs, and number and size of stems. This involves an increase of volume backscattering coefficient as well as attenuation of radar signal. However, the increase and decrease of backscatter caused by volume scattering and attenuation at the same time make radar signal reach saturation and then decrease when plant height is larger than 50 cm [START_REF] Baghdadi | Multitemporal observations of sugarcane by TerraSAR-X images[END_REF]. The dynamic of radar signal with the sugarcane height is slightly higher at 47° than at 31°. A dynamic of 5 dB for 47° and 2.5 dB for 31° is observed for cane heights between 0 and approximately 50 cm. Results show that σ° is strongly influenced by the soil moisture since a clear increase in the radar signal is observed after rainy episodes, in particular for young canes. Results showed that the radar signal is very dependent on the precipitation particularly at low and medium incidence angles and for young canes. Indeed, at low and medium incidences, the soil contribution (influenced by soil moisture) to total backscattering could be important for cane heights lower than 95 cm (Figure 2). The soil effects are small for images acquired at high incidence angles and for sugarcanes with vegetation well developed. The decrease in radar signal for harvested fields could be reduced of 3-4 dB on images acquired after rainy period.

Temporal backscatter and sugarcane harvest detection

This study also examined the potential of different TerraSAR-X incidence angles and polarizations for mapping sugarcane harvests. Harvested fields are easily detected on SAR images if the image acquisition date is close to harvest date (ideally less than two months). Indeed, the harvest involves a decrease in the signal that can reach 7 dB (VV-37°) if the observation radar is relatively close to the harvesting date (few days). The incidences of 17° and 58° allow only partially the detection of the harvest because the decrease of radar signal after the cut is about 3 dB (Figure 3). Figure 3b also shows that HH and VV polarizations are strongly correlated. The general trend is that the HH response is slightly higher than the VV (on the order of 1 dB). This confirms the effect of higher attenuation at the VV polarization for sugarcanes with a vertical structure [START_REF] Le Toan | Multitemporal and dual-polarization observations of agricultural vegetation covers by Xband SAR images[END_REF]. Figure 4 shows segments of TerraSAR-X images acquired between August 01, 2009 and December 11, 2009 in VV polarization and with an incidence angle of 37°. The interpretation of TerraSAR images shows that the difference between the backscatter of mature cane and of harvested cane is well pronounced at medium incidence angles (37°). The images show high σ° for mature canes and low σ° for harvested fields. The discrimination between harvested fields or young canes (less than two months old) and canes of more than two months old is better with TerraSAR at 31°-37°.

IV. CONCLUSIONS

The objective of this study was to analyze the behaviour of TerraSAR signal as a function of sugarcane height. The radar backscattering coefficient of sampled fields was studied using ground truth measurements of sugarcane height, SPOT images, and harvest dates. The increasing trend of σ° as a function of sugarcane height is observed until a height htvd around 50 cm, corresponding to total cane height around 155 cm (depends on incidence and polarization). The discrimination between young and mature canes is limited to fields harvested less than 2-3 months earlier (cane heights htvd between 0 and 50 cm).

This study also examined the potential of different TerraSAR-X incidence angles and polarizations for mapping sugarcane harvest. Harvested fields are easily detected on SAR images if the image acquisition date is close to harvest date (ideally less than two months). Indeed, the harvest involves a decrease in the signal that can reach 7 dB (VV-37°) if the observation radar is relatively close to the harvesting date (few days). The incidences of 17° and 58° allow only partially the detection of the harvest because the decrease of radar signal after the cut is about 3dB.

Results showed that the radar signal could be very dependent on the soil moisture particularly at low and medium incidence angles and for young canes. Indeed, at low and medium incidences, the soil contribution (influenced by soil moisture) to total backscattering could be important for cane heights lower than 95 cm. The soil effects are small for images acquired at high incidence angles and for sugarcanes with vegetation well developed. The decrease in radar signal for harvested fields could be reduced of 3-4dB on images acquired after rainy period.

The very high spatial resolution (metric) of TerraSAR-X offers great potential for mapping harvested sugarcane crop. This SAR provides a diagnosis suited to agricultural areas where the parcels are of small size. The spatial resolution of TerraSAR images, between 1 and 3 m (for Spotlight and Stripmap modes) are well suited for sugarcane production areas dominated by small farmers as in Reunion Island with fields areas of about 1 ha on average. These results appear promising for the development of simplified algorithms for monitoring sugarcane harvest regardless of meteorological conditions, which are the main limitation with optical sensors.

  the amplitude of backscattered signal for each pixel ( i DN ) into a backscattering coefficient ( i °σ ) in decibels. The calibration coefficient Ks (scaling gain value) varies within the range of 7.7 10 -6 to 1.2 10 -5 , depending on radar incidence angle (θ i ) and polarization (low values for cross-polarizations or high incidences). NEBN is the Noise Equivalent Beta Naught (2 i

Figure 1 .

 1 Figure 1. SPOT-5 image acquired over the study site in Reunion Island on October 21, 2008 (Red: band-3; Green: band-2; Blue: band-1). Reference sugarcane fields are outlined in blue.

Figure 2 .

 2 Figure 2. Radar backscattering coefficient as a function of plant height for VVand HV polarizations and incidence angle of 31° and 47°. The sugarcane height corresponds to the height of terminal visible dewlap (htvd).

Figure 3 .

 3 Figure 3. Temporal variation of TerraSAR signal for the reference sugarcane field 16. (a) 31 and 37°, (b) 17° and 59°. Field 16 was harvested on August 29, 2009.

Figure 4 .

 4 Figure 4. Comparison of several TerraSAR image segments for reference sugarcane fields (6.1, 6.2, and 16). All images were acquired at incidence of 37° and in VV polarization. Fields 6.1, 6.2, and 16 were harvested on September 10, October 03, and August 29, 2009, respectively.

Table I MAIN

 I CHARACTERISTICS OF TERRASAR-X IMAGES.

75 100 125 150 175 200 225 250 275 300 325 350 375 400 Height of sugarcane (cm) Backscattering coefficient (dB

			0
			VV-31°-
			-2
	Backscattering coefficient (dB)	-18 -16 -14 -12 -10 -8 -6 -4	20/12/2008 11/01/2009 24/02/2009 02/02/2009			Rainfall (mm)	60 0 20 40	20/12/08	Date (dd/mm/yy) 04/01/09 06/01/09 07/01/09 08/01/09 11/01/09 02/02/09	24/02/09
		20
			0				25			50	75	100	125	150
										Height of sugarcane (cm)
			0
					60
		Backscattering coefficient (dB)	-16 -14 -12 -10 -8 -6 -4 -2	Rainfall (mm)	0 20 40	20/12/08	Date (dd/mm/yy) 04/01/09 06/01/09 07/01/09 08/01/09 11/01/09 22/01/09 02/02/09 13/02/09	24/02/09	20/12/2008 11/01/2009 22/01/2009 02/02/2009
			-18							13/02/2009
			HV-31°-					24/02/2009
			20
			0				25	50		75	100	125	150	175	200
										Height of sugarcane (cm)
			0		30
		VV-47°-25 50 ) 20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 Date (dd/mm/yy) 0 10 20 14/12/08 23/12/08 24/12/08 25/12/08 27/01/09 14/02/09 15/02/09 18/02/09 08/06/09 27/06/09 Rainfall (mm)	29/06/09	30/06/09	14/12/2008 25/12/2008 27/01/2009 18/02/2009 08/06/2009 30/06/2009
			0
			-2
			-4
			-6
			-8
			-10
			-12
			-14
			-16
			-18
										HV-47°-
			20
			0				25			50	75	100	125	150
										Height of sugarcane (

cm) Backscattering coefficient (dB) 14/12/2008 25/12/2008 16/01/2009 27/01/2009 07/02/2009 18/02/2009 0 25 50 14/12/08 23/12/08 24/12/08 25/12/08 16/01/09 27/01/09 04/02/09 05/02/09 06/02/09 07/02/09 14/02/09 15/02/09 18/02/09 Date (dd/mm/yy) Rainfall (mm)

28/11/2008 28/12/2008 27/01/2009 26/02/2009 28/03/2009 27/04/2009 27/05/2009 26/06/2009 26/07/2009 25/08/2009 24/09/2009 24/10/2009 23/11/2009 23/12/2009 22/01/2010 21/02/2010 Date (dd/mm/yyyy)

				HH-31		VV-31			VV-37			Harvest date
		0
	Backscattering coefficient (dB)	-20 -16 -12 -8 -4		Field 16
			HH-17		VV-17		HH-59	VV-59		Harvest date
		0
	Backscattering coefficient (dB)	-20 -16 -12 -8 -4	26/02/2009	Field 16 28/03/2009	27/04/2009	27/05/2009	26/06/2009	26/07/2009	25/08/2009	24/09/2009	24/10/2009	23/11/2009	23/12/2009	22/01/2010	21/02/2010
							Date (dd/mm/yyyy)

ACKNOWLEDGMENT

The authors wish to thank CNES (French Space Agency) and DLR (German Space Agency) for kindly providing TerraSAR-X (proposal BOISSEZO_LAN0237). TerraSAR and SPOT images were obtained within the framework of Kalideos programme, set up by the CNES. Thanks are also due to Nathalie Boyer, Louis Paulin, and Raymond Nativel for their participation in the measurement surveys.


Frédéric Baup

Rémy Fieuzal

Benoît Mingam

Armand Lopes

Roger Fjortoft

Claire Marais- Sicre

Jean Francois Dejoux

Mehrez Zribi

Nicolas Baghdadi

Eric Ceschia

CONTRIBUTION OF ALOS, RADARSAT-2 AND TERRASAR-X RADAR DATA FOR MONITORING AGRICULTURAL SURFACES

ABSTRACT

In a changing climate context, it becomes increasingly important to accurately estimate the physical processes involved in the surface-atmosphere interactions in order to predict climate changes and its impact on ecosystems. Increase of human pressure and changes in land use management contribute to alter water and energy budgets as well as carbon sequestration in the soils [START_REF] Freibauer | Carbon sequestration in the agricultural soils of Europe[END_REF]. Therefore, it is essential 1) to improve our understanding of the different processes governing water, carbon and energy exchanges between the continental biosphere in anthropised areas and the atmosphere, 2) to monitor crop dynamics, soil and crop management and land use change in order to understand the effect of human pressure on biogeochemical cycles. In this study, we used complementary tools for this purpose including field data, remote sensing data and modelling at high spatial resolution to be compatible with the size of the agricultural plots in our area.

The satellite data analysis is performed in 2010 over a studied site located in the south west of France (Figure 1a) where meteorological conditions are steered by a temperate climate [START_REF] Béziat | Carbon balance of a three crop succession over two cropland sites in South West France[END_REF]. They are compared with ground data measurements acquired synchronously with all satellite acquisitions (Figure 1b). Wet and dry biomass distribution, vegetation height measurements are also performed during the crop phenological cycles. At regional scale, 350 fields are also monitored to identify crop species, soil management, soil tillage orientation… providing a wide dataset of different soil and vegetation status.

Although this work is included in an assimilation scheme for multispectral data (optical, thermal and microwaves), the analysis of the satellite data is here mainly focused on high spatial resolution (meter) microwave data to complete all the analyses performed with lower resolution sensors from ERS scatterometer or ASAR coarse mode data (over crops or natural vegetation) [3]- [START_REF] Dente | Assimilation of leaf area index derived from ASAR and MERIS data into CERES-Wheat model to map wheat yield[END_REF]. Combined to the high spatial resolution data, the contribution of multi-frequency data (L-, C-and X-Band) is assessed by using quasi-synchronous data provided by ALOS, RADARSAT-2 and TERRASAR-X satellites for monitoring the spatio-temporal variations of the surface properties (soil roughness, soil moisture, and vegetation status) associated to crop management (tillage, irrigation, fertilisation application…). In the following, results only focuses on surface soil moisture and soil roughness estimated from empirical and physicals methods. Soil management (ploughing, tillage, sowing…) and vegetation monitoring purposes are not discussed in this paper. Figure 3 shows an example of color-composed images over the region of interest in April and May 2010 processed by combining X, C and L-band images respectively displayed on red, green and blue colors. All images are acquired in HH polarization. Author-produced version of the article published in IEEE International Geoscience and Remote Sensing Symposium -Proceedings, 2011.

Figure 1 :

 1 Figure 1: Presentation of the studied sites superimposed with the swaths of ALOS (AP-L), TERRASARX (TS-X) and RADARSAT2 (RS-C) images (a) and their temporal sampling (b) in 2010. Optical and thermal acquisitions provided by SPOT, FORMOSAT2 and Landsat5-7 satellites are not used in this study.

Figure 2 :

 2 Figure 2: Example of one of the 37 temporal evolution of SSM (mean and standard deviation are plotted). This site is corn cultivated (green curve), and the soil is mainly composed of clay (>50%).

Figure 3 :

 3 Figure 3: Color-composed images over the studied area in April (2010-04-14) and May (2010-05-01) processed by combining X, C and L band images respectively displayed on red, green and blue colors.

Figure 4 :

 4 Figure 4: Comparisons of simulated backscattering coefficients by using the IEM model (a) and Dubois Model (b) versus satellite observations. C-band results are here presented at HH polarization.

Author-produced version of the article published in IEEE International Geoscience and Remote Sensing Symposium -Proceedings, 2011.

This approach is consistent with previous studies performed in microwave multi-frequency analyses over crops ( [10]- [11]) and with the launch of future multi sensors satellites constellations which will combine optical and radar data (Sentinel and Pléiades projects [12], [13]).


Mark Schmidt

Karteek Alahari

Generalized Fast Approximate Energy Minimization via Graph Cuts: α-Expansion β-Shrink Moves

We present α-expansion β-shrink moves, a simple generalization of the widely-used αβswap and α-expansion algorithms for approximate energy minimization. We show that in a certain sense, these moves dominate both αβ-swap and α-expansion moves, but unlike previous generalizations the new moves require no additional assumptions and are still solvable in polynomial-time. We show promising experimental results with the new moves, which we believe could be used in any context where α-expansions are currently employed.

Introduction

We focus on the problem of finding the most probable configuration in a pairwise Markov random field over discrete variables, a fundamental problem in the study of graphical models. This is equivalent to minimizing the sum of a set of unary and pairwise energy functions defined over a set of discrete variables. Problems of this type arise in many applications, but in general it is NP-hard to solve these problems even in the case of binary variables (see Kolmogorov and Zabih, 2002, Theorem 4.2). A classical method for computing an approximate solution to this problem is Besag's iterated conditional mode (ICM) algorithm [START_REF] Besag | On the statistical analysis of dirty pictures[END_REF], but due to the local nature of this method it may become stuck in a poor local optimum.

In the special case where the variables are binary and the pairwise energy functions satisfy a submodularity condition, it is possible to solve this problem in polynomial-time [START_REF] Hammer | Some network flow problems solved with pseudo-boolean programming[END_REF][START_REF] Greig | Exact maximum a posteriori estimation for binary images[END_REF][START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF]. This has motivated the αβ-swap and α-expansion moves proposed by Boykov * INRIA/ENS/CNRS UMR 8548 et al. (1998,1999), which we review in §2. These methods find an approximate solution to a non-binary problem satisfying a submodularity condition over pairs of states (or triplets of states in the case of α-expansions) by solving a sequence of binary submodular problems. The experiments of [START_REF] Szeliski | A comparative study of energy minimization methods for markov random fields with smoothness-based priors[END_REF] show that these algorithms perform well compared to competing approximate minimization methods at minimizing the energy functions arising from problems in computer vision, such as stereo image matching, creating photo montages, and image restoration. Further, this study found that α-expansions are often faster than competing methods with similar performance for these problems (such as variational message-passing algorithms like tree-reweighted belief propagation). While originally introduced in the context of computer vision, these algorithms and their generalizations have also proved to be effective in other domains, such as protein structure prediction [START_REF] Gould | Alphabet SOUP: A framework for approximate energy minimization[END_REF].

In a sense that we formally define in §3, αβ-swaps and α-expansions dominate the classic ICM algorithm. However, although α-expansions often lead to better experimental results than αβ-swaps [START_REF] Szeliski | A comparative study of energy minimization methods for markov random fields with smoothness-based priors[END_REF], under our definition of dominance neither of these more advanced methods dominates the other. In this paper we propose a new type of move, α-expansion β-shrink moves ( §4). These move are a simple generalization of both αβ-swap and α-expansions, that dominates them both. Although we delay discussion of other generalizations of these moves to the end ( §6), we note that unlike previous generalizations the new moves require no additional assumptions beyond those needed to apply α-expansions, and the moves can be computed in polynomial-time by solving a binary submodular problem defined on the original graph structure. Thus, we can use them in place of α-expansions for applications where these moves are currently used. Our experiments on standard test problems from the field of computer vision ( §5) show that the new moves can lead to improved performance.

The problem that we formally address is

max x∈{1,2,...,N } p p(x) ∝ i∈V φ i (x i ) (i,j)∈A φ ij (x i , x j ),
where V and A are the vertices and arcs of a graph (V, A), while the potentials φ i (x i ) and φ ij (x i , x j ) map assignments of subsets of the discrete vector x to nonnegative values. If we define the real-valued energy functions

E i (x i ) = -log φ i (x i ), E ij = -log φ ij (x i , x j ),
then finding the optimal assignment is equivalent to the following energy minimization problem:

min x∈{1,2,...,N } p i∈V E i (x i ) + (i,j)∈A E ij (x i , x j ), (1)
In general, solving this optimization problem is NPhard (see Kolmogorov and Zabih, 2002, Theorem 4.2), and a common approach to finding an approximate minimizer is with an iterative descent algorithm. The input to each iteration of these algorithms is a particular configuration of the variables x, and at each iteration an iterative descent method finds a configuration that minimizes the energy among a set M(x) of possible 'moves', i.e. the next iteration is an element of argmin

y∈M(x) i∈V E i (y i ) + (i,j)∈A E ij (y i , y j ). ( 2
)
We have two conflicting desiderata on the set of moves M(x k ): we would like this set to be as large as possible, but we would like to able to efficiently find the optimal move.

Iterated Conditional Mode

In the classic ICM algorithm [START_REF] Besag | On the statistical analysis of dirty pictures[END_REF], a node j is selected and we replace x j with a value that maximizes the conditional probability p(x j |x -j ), where x -j are the states of all variables except j. In the framework of energy minimization, this can be viewed as an iterative descent method where the elements y of the set of possible moves have the form

y i ← γ ∈ {1, 2, . . . , N } if i = j, x i otherwise.
That is, the current state of node j can be replaced by any other possible state, a form of coordinate descent. We use M I j (x) to denote the set of all y of this form. With this definition of the move space, the iterative descent update (2) for the ICM move given x and j simplifies to argmin

y∈M I j (x) E j (y j |x -j ), (3)
where we will find it convenient to define the conditional energy of a variable i given a set of variables a as

E i (x i |x a ) = E i (x i ) + j|j∈a,(i,j)∈A E ij (y i , x j ) + j|j∈a,(j,i)∈A E ij (x j , y i ). (4) While exp(-E i (x i |x -i )) ∝ p(x i |x -i ),
this is a slight abuse of the conditioning notation since for other conditioning sets a it ignores factors that depend on variables besides i and those in a. Clearly, we can efficiently compute the optimal ICM move by simply testing each E j (y j |x -j ).

αβ-Swaps

We say that a pairwise energy function E ij defined on binary variables is submodular if it satisfies the inequality 1

E ij (1, 1) + E ij (2, 2) ≤ E ij (2, 1) + E ij (1, 2). (5)
This type of pairwise energy prefers the neighbouring variables to take the same state. In the special case of binary variables where all pairwise energies are submodular, the optimal solution to (1) can be computed in polynomial-time as a minimum-cut problem, see [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF]. 2

In some non-binary problems we have, for all combinations of states α and β, that the pairwise energies

E ij satisfy E ij (α, α) + E ij (β, β) ≤ E ij (β, α) + E ij (α, β). (6)
That is, when restricted to any two states α and β, the energy function is submodular. Though we can no longer guarantee that we can find the optimal solution in polynomial-time given only this restriction, it

1 Submodularity is normally defined as a property of functions on sets, and its use here is because it is equivalent to submodularity of a function that takes the set of variables labeled 2 and returns the corresponding Eij (Kolmogorov and Zabih, 2002, §7).

2 There also exist several other notable cases where it possible to compute the solution in polynomial-time, such as the case where the pairwise energies are convex [START_REF] Ishikawa | Exact optimization for Markov random fields with convex priors[END_REF], or in the case of general energies where the graph structure has low treewidth (Koller and Friedman, 2009, §13) or is outer-planar [START_REF] Schraudolph | Efficient exact inference in planar Ising models[END_REF]. However, these are not our focus and a full discussion of this extensive literature is outside the scope of the current work.

does allow us to take advantage of the ability to efficiently solve binary submodular problems in order to implement a more powerful descent move.

In particular, given two states α and β, the set of moves M S αβ (x) associated with the αβ-swap move introduced by [START_REF] Boykov | Markov random fields with efficient approximations[END_REF] are of the form:

y i ← α or β if x i = α or x i = β, x i otherwise.
That is, the move can simultaneously change any combination of nodes labeled α to β, and any combination of nodes labeled β to α. In this case, the iterative descent update is a solution of the problem argmin

y∈M S αβ (x) i∈V|xi∈{α,β} E i (y i |x -αβ ) + (i,j)∈A|xi,xj ∈{α,β} E ij (y i , y j ),
where we again make use of our definition of the conditional energy ( 4) and where we have used x -αβ to reference the states of the variables not labeled α or β. This is a binary problem over the subgraph induced by the nodes labeled α or β, and under condition (6) this update can be computed in polynomial-time because all edges in the induced subgraph are submodular.

α-Expansions

A closely-related set of moves later proposed by [START_REF] Boykov | Fast Approximate Energy Minimization via Graph Cuts[END_REF] are α-expansions. Here, we choose a state α and we can use α to replace the current state of any variable. Thus, this set of moves

M E α (x) is of the form y i ← α if x i = α, α or x i otherwise.
We can write the optimal α-expansion as the solution to the problem

argmin y∈M E α (x) i∈V|xi =α E i (y i |x α ) + (i,j)∈A|xi =α,xj =α E ij (y i , y j ),
where we use x α to reference variables labeled α. This is again a binary problem, this time on the subgraph induced by those nodes not labeled α. However, condition ( 6) is no longer sufficient to guarantee that the edges in the induced subgraph are submodular since each pairwise term involves the three states α, x i , and x j (which may all be different). Nevertheless, it is sufficient to ensure that

E ij (α, α) + E ij (γ 1 , γ 2 ) ≤ E ij (γ 1 , α) + E ij (α, γ 2 ), (7)
for all combinations of states α, γ 1 , and γ 2 . This is a stronger condition than (6), which corresponds to the special case where γ 1 = γ 2 . If E ij (α, α) = 0 for all α, this is the triangle inequality.

Local Dominance of Iterative Algorithms

Some authors analyze the approximation accuracy of local optima found by iterative descent methods such as α-expansions (Veksler, 1999, §4.3.4), but the guarantees on the approximation accuracy tend to be very conservative since these global bounds must hold over all possible optima. Further, they often do not provide guidance in choosing between iterative descent methods. For example, a small modification to the analysis of (see Veksler, 1999, §4.3.4) gives us an identical global bound for the ICM algorithm.

Instead of focusing on such a global analysis, we instead focus on a simple local analysis of the methods. Specifically, we are interested in the improvement in the objective function that can be obtained in a single iteration with different types of moves, given the same input to the iteration. From this viewpoint, we say that a move set M A (x) dominates a move set M B (x) if:

1. For any energy function E(x) and any input configuration x, optimizing over M A (x) never gives a higher energy than optimizing over M B (x):

∀ E,x , min y∈M A (x) E(x) ≤ min y∈M B (x) E(x). (8)
2. There exists energy functions E(x) and input configurations x where optimizing over M A (x) gives a strictly lower energy than optimizing over M B (x):

∃ E,x , min y∈M A (x) E(x) < min y∈M B (x) E(x). (9)
This is an intuitive notion of the relative local improvement that is made by different moves, but it can also be used as a measure of the relative 'strengths' of local optima found by different types of moves; if the move set M A (x) dominates the move set M B (x), then no move in the set M B (x) will ever be able to improve on a local minimum with respect to move set M A (x), but the optimal move in move set M A (x) may be able to escape from a local minima with respect to move set M B (x).

To establish dominance relations between various search moves, it will be convenient to define versions

α β γ α β γ β γ α α β γ
Figure 1: From left to right: Initial labeling, labeling after αβ-swap, labeling after α-expansion, labeling after α-expansion β-shrink. The optimal labeling of the α pixels is outlined by a white triangle, and is achieved from the initial labeling by one α-expansion β-shrink move.

of the moves that do not depend on any parameters except x. Thus, we make the following definitions:

M I (x) = ∪ j M I j (x), M S (x) = ∪ α,β M S αβ (x), M E (x) = ∪ α M E α (x).
That is, these expanded move spaces simply search over all possible values of the parameter(s). Under these definitions, it is straightforward to establish the following relationships between the methods:

Proposition 1. The αβ-swap move set M S (x) dominates the ICM move set M I (x). The α-expansion move set M E (x) dominates the ICM move set M I (x).

Proof. Assume an arbitrary energy function E and configuration x, and let y be the optimal move in M I (x), which must be an element of M I j (x) for some j. Then the optimal ICM move is in M S αβ (x) with α = x j and β = y j , which establishes (8). To establish (9), it is sufficient to construct an energy function E along with a configuration x where the optimal solution involves changing two elements of x. For example, take a two variable problem with x = (1, 1) where the optimal solution is (2, 2). Establishing that α-expansions dominate ICM can be done similarly.

Proposition 2. The α-expansion move set M E (x) does not dominate the αβ-swap move set M S (x). The αβ-swap move set M S (x) does not dominate the αexpansion move set M E (x).

Proof. To establish that α-expansions do not dominate αβ-swaps, it is sufficient to construct an energy function E and configuration x where the optimal solution is obtained by a single αβ-swap but can not be obtained by a single α-expansion. For example, take a two-variable problem with x = (1, 2) where the optimal solution is (2, 1). We can similarly show that αβ-swaps do not dominate α-expansions by taking a three variable problem with x = (1, 2, 3) where the optimal solution is (1, 1, 1).

Thus, based on the definition of dominance discussed in this section, we should prefer both αβ-swaps and αexpansions to ICM. However, local dominance does not guide us in selecting between these two more advanced methods.

α-Expansion β-Shrink Moves

Towards the goal of developing a method that dominates both these methods, consider a generalization of both αβ-swaps and α-expansions where the set of moves M G αβ (x) are of the form:

y i ← α or β if x i = α α or x i otherwise
That is, we allow any node not currently labeled α to take the value α, but in addition we allow the nodes currently labeled α to take the value β (α expands everywhere, but is shrunk by β). We call these αexpansion β-shrink moves, and Figure 2.3 illustrates an instance of the move.

Under the definition of the previous section, we have the following result for an analogously defined M G (x):

Proposition 3. The α-expansion β-shrink move set M G (x) dominates the αβ-swap move set M S (x). The α-expansion β-shrink move set M G (x) dominates the α-expansion move set M E (x).

Proof. To establish (8), we note that M G (x) includes all αβ-swaps and α-expansions as special cases. Further, it also includes moves that are not instances of either of these moves (i.e. cases where α expands into more than one state, but β shrinks into α). Thus, to establish (9) we can simply take a three variable problem with x = (1, 2, 3) where the optimal solution is (2, 1, 1).

Note that the argument above can also be used to show that the new move set M G (x) also dominates the move set M S ∪ M E , the union of αβ-swap and α-expansion moves.

We can write the optimal α-expansion β-shrink move as a solution to the problem argmin

y∈M G αβ (x) i∈V E i (y i ) + (i,j)∈A E ij (y i , y j ), (10)
Note that unlike the previous move sets we have discussed where the move is computed by solving a problem involving conditional energies over an induced subgraph, this move involves solving a problem with the original unary energies on the original graph structure. In some scenarios, these properties might make the generalized move simpler to implement than αβswaps or α-expansions. Further, if the original graph has a special structure, it allows the use of specialized codes for solving binary minimum-cut problems with this structure. For example, [START_REF] Delong | A scalable graph-cut algorithm for nd grids[END_REF] propose a method for the special case of problems with very large grid structures. The simple form of the subproblem may also simplify the implementation of dynamic graph cuts (Kohli and Torr, 2007), where the similarity between subproblems when using the same value of α (and β) is used to substantially speed up the computation.

There are many possible moves we could define that would dominate one or more of the methods discussed in §2, and we return to this topic in §6. However, an appealing property of α-expansion β-shrink moves is the following:

Proposition 4. If for each edge E ij condition (7) holds for all γ 1 and γ 2 , then all edge energies E ij in problem (10) are submodular.

Proof. It is sufficient to show this for an arbitrary E ij and any possible assignment to x i and x j . If x i = α and x j = α, then we require

E ij (α, α) + E ij (x i , x j ) ≤ E ij (x i , α) + E ij (α, x j )
as with α-expansion moves, which is (7) with γ 1 = x i and γ 2 = x j . If x i = α and x j = α, then we require

E ij (α, α) + E ij (β, β) ≤ E ij (β, α) + E ij (α, β) as with αβ-swap moves, which is (7) with γ 1 = β and γ 2 = β. If x i = α and x j = α, then we require E ij (α, α) + E ij (β, x j ) ≤ E ij (β, α) + E ij (α, x j )
, which is (7) with γ 1 = β and γ 2 = x j . The remaining possibility is similar.

This implies that an optimal α-expansion β-shrink move can be computed by solving a minimum-cut problem, using (for example) the construction of [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF]. Thus, an optimal αexpansion β-shrink move can be computed in polynomial time under the same condition required to use α-expansions. Further, computing an optimal α-expansion move and computing an optimal αexpansion β-shrink move have the same worst-case time-complexity. This would indicate that the new moves could be used in place of α-expansions in any of the many applications where these moves are currently used.

Problems with Many States

In the computer vision problems that originally motivated the moves based on minimum-cuts, the number of states N may be non-trivial, since the states may represent a discretization of a continuous value. For example, [START_REF] Szeliski | A comparative study of energy minimization methods for markov random fields with smoothness-based priors[END_REF] evaluate the performance of these moves on problems with 256 states, each state representing an intensity level in an image. Thus, in addition to their better empirical performance, α-expansions may be preferred over αβ-swaps simply because the number of possible values of α is N , while the number of (non-exchangeable) combinations of α and β is N (N -1)/2. This quadratic scaling in terms of N would also seem to be a problem for α-expansion β-shrink moves. However, if we consider any strategy for choosing a value of β given α, then we still obtain an algorithm that dominates α-expansions (assuming we do not always choose α = β, corresponding to the special case of αexpansions). Thus, the new moves can be modified to have the same scaling with N as α-expansions. That is, if N is very large we can consider variants where β is function of α so that we only consider N moves instead of N 2 , and these variants would still dominate α-expansions. We empirically evaluate three possible strategies for selecting β given α in §5.

Truncation for Non-Submodular Potentials

In many problems condition (7) is not satisfied. In these cases, a widely-used approach is to modify the potentials to be submodular, in such a way that an optimal move with the modified energy is guaranteed to not increase the original energy [START_REF] Rother | Digital tapestry[END_REF].

In the case of α-expansions, one way to construct such a modified energy is by replacing each

E ij (x i , x j ) with Ēij (x i , x j ) = min{E ij (x i , x j ), E ij (α, x j ) + E ij (x i , α) -E ij (α, α)}.
Condition (7) holds with this modified energy by construction. Further, the optimal α-expansion with this modified energy does not increase the original energy, since the modified energy simply decreases the energy of the current assignment (x i , x j ). As discussed in [START_REF] Rother | Digital tapestry[END_REF], we can alternately increase E ij (x i , α) or E ij (α, x j ) to make condition (7) satisfied maintaining the descent property of the moves.

We can define a modified energy function with similar properties in the case of α-expansion β-shrink moves, though maintaining the descent property requires a slightly more complicated construction. Although there are many possible constructions, we describe one here. If x i = α and x j = α, then as before we take:

Ēij (x i , x j ) = min{E ij (x i , x j ), E ij (α, x j ) + E ij (x i , α) -E ij (α, α)}.
If x i = α and x j = α, then we take:

Ēij (α, α) = min{E ij (α, α), E ij (α, β) + E ij (β, α) -E ij (β, β)}.
If x i = α and x j = α, then we take:

Ēij (α, β) = max{E ij (α, β), E ij (α, α) + E ij (x i , β) -E ij (x i , α)}.
Finally, if x i = α with x j = α we take:

Ēij (β, α) = max{E ij (β, α), E ij (α, α) + E ij (β, x j ) -E ij (α, x j )}.
The other terms in the energy function are unchanged. If ( 7) is already satisfied, then the modified energy under this construction is identical to the original energy. This construction maintains the appealing property that any move that does not increase the modified energy will not increase the original energy.

Experiments

To empirically evaluate the performance of the new moves, we performed several experiments on the nonbinary data sets examined by [START_REF] Szeliski | A comparative study of energy minimization methods for markov random fields with smoothness-based priors[END_REF]. These data sets are summarized in Table 1, and we extracted the terms in the energy functions from the code available online.3

Our experiments compared the following methods:

• αβ-Swap: performing αβ-swap moves in the order β = 1, 2, . . . , N in an outer loop and α = β+1, β+ 2, . . . , N in an inner loop.

• α-Expansion: performing α-expansion moves in the order α = 1, 2, . . . , N .

• Random β: performing α-expansion β-shrink moves in the order α = 1, 2, . . . , N with β selected randomly among {1, 2, . . . , N }.

• β = α -1: performing α-expansion β-shrink moves in the order α = 1, 2, . . . , N with β set to max{1, α -1}. • β = α + 1: performing α-expansion β-shrink moves in the order α = 1, 2, . . . , N with β set to min{N, α + 1}.

• All β: performing α-expansion β-shrink moves in the order β = 1, 2, . . . , N in an outer loop and α = 1, 2, . . . , N in an inner loop.

Note that the iterations of the first and last method are much more expensive for large N because they consider O(N 2 ) combinations of α and β, while the remaining methods only consider N moves.

In our first experiment, we initialized all variables to the first state and ran each method until the energy did not change between iterations. We used the truncation described in §4.2 for problems that did not satisfy (7) for all triplets of states. In Table 2, we show the energy of the local minima obtained divided by the energy of the local minimum with respect to α-expansion moves.

In this table, a value of 1 indicates that the energy was identical to the energy obtained by α-expansion moves, and we use 1.0000 if the energy is close but not identical. In this experiment, the new moves with all β obtained the lowest energy on 6 of the 7 data sets, and strictly so in 4 of these cases (the exceptions were Teddy where using a random β lead to a lower score, and the montage data sets where other methods reached the same energy). Among the remaining methods, the more computationally efficient strategy of simply setting β = α + 1 obtained the lowest energy on 5 of the 7 data sets.

In our second experiment we used the local minimum with respect to α-expansion moves as the initialization, and tested whether the methods that only consider N moves could escape from this local minimum. Table 3 shows the energies of the local optima obtained with this initialization. Here, we see that even with these simple choices of β that the new moves are able to escape the local minimum with respect to α-expansions for 5 out of the 7 data sets. In this experiment choosing β = α + 1 gave the lowest energy on all data sets, and gave a strictly lower energy on 4 of the 7. Selecting β = α + 1 is somewhat intuitive, the moves can Although an improved configuration was found on at least one data set for each of the three tasks, the largest improvements were seen in the image restoration/inpainting tasks. This is likely due to the 256 possible labels, which is much larger than the other tasks. In Figure 5, we show the restoration of the house data where the largest improvement was observed. In this figure we see a noticeable difference between the local minima with respect to the different moves. This visual difference is also reflected quantitatively; the new moves lead to a reconstruction error (ℓ 1 -norm distance to the original image) that is approximately two-thirds that of the local minimum with respect to α-expansions (reduced from 1303537 to 911207).

Discussion

There have been a variety of other generalizations of αβ-swaps and α-expansions proposed in the literature, and an exhaustive list of this literature would be outside our scope. However, among the generalizations we are aware of these other moves either (i) require additional assumptions, or (ii) can not be solved in polynomial-time. For example, rangeswap and range-expansion moves generalize αβ-swaps and α-expansions (respectively), but these can not be computed in polynomial-time without further assumptions [START_REF] Veksler | Graph cut based optimization for MRFs with truncated convex priors[END_REF][START_REF] Gould | Alphabet SOUP: A framework for approximate energy minimization[END_REF][START_REF] Kumar | Improved moves for truncated convex models[END_REF]. [START_REF] Kolmogorov | Minimizing nonsubmodular functions with graph cuts-a review[END_REF] consider implementing minimum-cut methods wihout restrictions on the energy functions, but in general their algorithm can not compute the optimal move. [START_REF] Lempitsky | Fusion Moves for Markov Random Field Optimization[END_REF] consider a general class of moves called fusion moves. These moves consider two full configurations of the variables x and z, and the set of moves are of the form y i ← x i or z i .

That is, we can replace an arbitrary number of elements of x by the corresponding elements from some alternative configuration z. The α-expansion β-shrink moves proposed here are a special case of this type of move, where z i is set to β if x i = α and z i is set to α otherwise. However, while Proposition 4 shows that this special case is solvable in polynomial-time under condition (7), [START_REF] Lempitsky | Fusion Moves for Markov Random Field Optimization[END_REF] do not give restrictions that would allow an optimal general fusion move to be solved in polynomial time. Indeed, because of the generality of fusion moves, any such restriction would need to be much more restrictive than (7).

Although in this paper we have focused on the case of minimizing unary and pairwise energies, it is possible to apply αβ-swaps and α-expansions to certain classes of higher-order energies (Kohli et al., 2007). We expect that the new moves discussed in this paper can also be extended to these scenarios, and indeed we believe that the new moves could be used in place of α-expansions in any of the many applications where these moves are currently used. Finally, we note that the implementation of the new moves will be made available online.

Table 1 :

 1 Data sets from[START_REF] Szeliski | A comparative study of energy minimization methods for markov random fields with smoothness-based priors[END_REF] used in the experiments.

	Name	Task	Nodes	Edges	States
	Family	Montage	425632 849946	5
	Pano	Montage	514080 1026609	7
	Tsukuba Stereo	110592 220512	16
	Venus	Stereo	166222 331627	20
	Teddy	Stereo	168750 336675	60
	Penguin Restoration 21838	43375	256
	House	Restoration 65536	130560	256

Table 2 :

 2 Energy of local minima with respect to iterative descent methods beginning from all variables set to state 1, divided by the energy of the local minima with respect to α-expansion moves.Name αβ-Swap α-Expansion Random β β = α -1 β = α + 1 All β

		Family	1.0203	1	0.9998	1	0.9998	0.9998
		Pano	1.3182	1	1.0006	1	1	1
		Tsukuba	1.0315	1	1.0012	1	1.0000	1.0000
		Venus	1.8561	1	1.0015	0.9992	0.9979	0.9968
		Teddy	1.0037	1	0.9998	1	1.0007	0.9999
		Penguin	1.1283	1	1.0037	0.9936	0.9793	0.9758
		House	0.7065	1	0.7841	0.9973	0.7038	0.7032
	Table 3: Relative energy of local minima for different
	choices of β when initializing with a local minima with
	respect to α-expansions.
	Name	Random β β = α -1 β = α + 1
	Family	0.9998	1	0.9998
	Pano	1	1	1
	Tsukuba	1	1	1
	Venus	1.0000	0.9992	0.9979
	Teddy	1	1	0.9999
	Penguin	0.9998	0.9902	0.9775
	House	0.8050	0.9971	0.7038
	be viewed as 'prematurely' expanding the next value
	of α into the region occupied by the current value of
	α.

http://vision.middlebury.edu/MRF/

Acknowledgements

We would like to thank Olivier Duchenne and Daniel Tarlow for valuable discussions, and the anonymous reviewers for helpful comments that improved the paper. Mark Schmidt is supported by the SIERRA grant from the European Research Council (SIERRA-ERC-Figure 2: From left to right: Initial degraded image with missing areas, local minimum of energy function with respect to α-expansions, improved local minimum found by α-expansion β-shrink moves with β = α + 1. 239993). Karteek Alahari is supported by the Quaero Programme, funded by the OSEO.


Marc Lambert

Frédéric Nouguier

Riadh Zorgati

Eddy-current modeling of a continous conductivity profile resulting from a diffusion process

Keywords: coils, eddy current testing, electrical conductivity measurement, inverse problems, probes

This paper deals with the modelisation of Eddy-Current measurements over combustion turbine blade coatings affected by depletion of aluminum. First, we model the response of an Eddy-Current coil over a layered metallic structure with a top over-aluminized coating by extending the analytical Uzal-Rose's model for one hyperbolic tangent conductivity profile to a conductivity profile using two hyperbolic tangents for taking inward and outward depletion of aluminum inside the coating into account. Results obtained with this model are similar to those obtained with a numerical multi-layer model but with a reduced computing time.

I. INTRODUCTION

The modelization of the response of a sensor in the presence of conductive material of different geometry is of great interest as exemplified in [START_REF] Theodoulidis | Interaction of an eddy-current coil with a right-angled conductive wedge[END_REF] where the interaction of an Eddy-current coil with a right-angled conductive wedge is dealt with, in [START_REF] Skarlatos | Impedance calculation of a bobbin coil in a conductive tube with eccentric walls[END_REF] in which the impedance calculation of a bobbin coil in a conductive tube with eccentric walls is treated. Some others works proposed some ideas to work with arbitrary shaped coil as in [START_REF] Ishida | Magnetic field analysis of an arbitrary shaped coil using shape functions[END_REF] whereas in [START_REF] Mayer | Inductance of coil in nonhomogenous surrounding[END_REF] the case of the inductance of a coil in a nonhomogenous surrounding is dealt with. Here the eddy-current modeling of the interaction of a coil with a continous conductivity profile resulting from a diffusion process is proposed.

In the electrical industry, power production using gas turbines has to be developed in a significant way. In order to achieve a good economical performance, the thermal efficiency of gas turbines has been improved by increasing the firing temperatures. Large increase by hundreds of degrees of gas temperature has been made possible using high technologies regarding both cooling of the hot section components and protection against hot corrosion and high temperature oxidation of the surface of some components, like blades in high-pressure section.

This protection involves specific coatings using Ni-based super-alloys where aluminum and chrome, contained in precise proportion, form protective Al 2 O 3 or Cr 2 O 3 surface oxides against corrosive attacks. Yet, advanced blades protected by a top over-aluminized coating are affected by a significant inward and outward depletion of aluminum. Inward depletion is due to migration of aluminum from the coating to the blade substrate, induced by differences of conductivities between the Marc Lambert and Frédéric Nouguier and Riadh Zorgati are with Département de Recherche en Electromagnétisme, Laboratoire des Signaux et Systèmes UMR8506 (CNRS-Supélec-Univ Paris-Sud), 91192 Gif-sur-Yvette cedex, France two layers, in order to restore the thermodynamical equilibrium. Outward depletion is due to migration of aluminum in the coating towards the surface to reform the protective oxides layer which spalls off during hostile conditions operation. Unfortunately, when losing aluminum, advanced top overaluminized coating degrades and, if not timely repaired, this might lead to blade failure. Benefits of a higher operating temperature is then counteracted by possible limited operating life and high replacement cost. It is hence necessary to predict remaining blade coating life by assessing condition and quality of the coating. Expected reduction of shut-downs caused by blade coating failure then should contribute to optimize the maintenance of combustion turbines.

By providing information related to the electrical conductivity depth profile, the multi-frequency Eddy Current (EC) non-destructive technique can be useful for reaching that goal. The paper is organized as follows. The configuration at hand is described in section II. Section III is devoted to the solution of the direct problem which aims at modeling the response of an EC coil above a layered metal structure protected by a top over-aluminized coating degraded by inward and outward diffusion of aluminum. To take such diffusion process of aluminum inside the blade into account, we develop a model where the conductivity profile, which varies continuously with depth in the blade, is described via two hyperbolic tangents. We hence extend the Uzal-Rose's model [START_REF] Uzal | The impedance of eddy current probes above layered metals whose conductivity and permeability vary continuously[END_REF]- [START_REF] Uzal | Experimental determination of the near-surface conductivity profiles of metal from electromagnetic induction (eddy current) measurements[END_REF] based on one hyperbolic tangent conductivity profile. Comparison of results obtained with this model and with those obtained with a numerical multi-layer model, based on discretization of the conductivity profile into a number of layers of constant conductivity are presented in section VI.

II. EDDY CURRENT MODEL FOR TURBINE BLADE COATING

Save its convex and concave parts, a blade turbine made of a substrate and a coating is considered as a doubly layered isotropic, linear, homogeneous, non-magnetic and conductive half-space. The top over-aluminized coating protects the substrate by forming a thin adherent oxide and degrades in-service by depletion of aluminum as commented upon before. The conductivity of the coating layer is a continuous function of depth z while the conductivity of the substrate layer is a given constant (cf. Fig. 1).

Above the blade, an EC probe with a circular coil of rectangular cross-section with axis perpendicular to the blade surface is operated at frequencies in the range of a few hundreds of kHz to a few tens of MHz using a scanning frequency procedure. In this section, we aim at solving the related direct problem consisting of modeling the EC impedance coil response over a layered metallic blade structure and taking the double diffusion of aluminum into account when the conductivity depth profile is given.

Following the approach in [START_REF] Uzal | The impedance of eddy current probes above layered metals whose conductivity and permeability vary continuously[END_REF]- [START_REF] Theodoulidis | Analytical solutions in eddy current testing of layered metals with continuous conductivity profiles[END_REF], based on the Cheng-Dodd-Deeds model [START_REF] Dodd | Analytical solutions to eddy-current probe-coil problems[END_REF], [START_REF] Cheng | General analysis of probe coils near stratified conductors[END_REF], the two main steps for establishing the model are determining the vector potential and then calculating the coil impedance from Z = V /I after computing the electrical field by E (r) = -jωA (r). From the fundamental result which gives the impedance of a coil located above a two-layer metal whose upper layer conductivity varies continuously with depth, one can calculate the variation of impedance for a large number of configurations.

III. THEORETICAL FORMULATION OF THE PROBLEM

A. General formulation for a continuous profile

The first step consists of calculating the magnetic vector potential in air expressed as A (r) = A Source (r)+A ec (r), where A Source (r) is the primary potential due to the time-harmonic source current of the coil I exp(jωt) in an unbounded space and A ec (r) is the secondary potential originated from the EC induced within the conductive material. Assuming axial symmetry and choosing the cylindrical coordinates centered on the axis of the coil leads to a scalar problem, all vector quantities being expressed with the only non-null component, the azimuthal one, as A (r) = A (r, z) u θ . The expression of the primary field A Source (r, z) due to a N turn coil can be found in [START_REF] Cheng | General analysis of probe coils near stratified conductors[END_REF] and the secondary field A ec (r, z) is solution of the diffusive equation (see [START_REF] Schmidt | Estimating the eddy-current modeling error[END_REF] for a justification of eddy currents model for the Maxwell equations):

∆A ec (r, z) = jωµ 0 σ (z) A ec (r, z) (1)
with

∆ = ∂ 2 ∂r 2 + 1 r ∂ ∂r + ∂ 2 ∂z 2 - 1 r 2 , (2)
where the continuously varying conductivity profile is given by σ(z), µ 0 being the permeability of the air and ω = 2πf the angular frequency and f the frequency (in Hz). Let us defined the skin depth δ (z) as

δ (z) = 2 ωµ 0 σ(z) . (3)
Equation ( 1) can be solved by separation of variables letting A ec (r, z) = R(r)W (z). The two ordinary differential equations for the r and the z dependence, respectively, follows as:

∂ 2 ∂r 2 R(r) + 1 r ∂ ∂r R(r) + a 2 - 1 r 2 R(r) = 0 (4) ∂ 2 ∂z 2 W (z) = a 2 + jωµ 0 σ(z) W (z). ( 5
)
Whilst the solution of ( 4) is classically obtained using first-and second-kind Bessel functions of the first order, the solution of ( 5) constitutes a more challenging task and can be performed by an analytical method for some particular conductivity profiles. Following the path presented in [START_REF] Theodoulidis | Analytical solutions in eddy current testing of layered metals with continuous conductivity profiles[END_REF] the general solution of ( 5) is given by

W (z) = CF 1 (f (z)) + BF 2 (f (z)) (6)
where F 1 and F 2 are functions related to the choice of the profile σ (z), f (z) originating from possible transformations of the independent variable that take place during the solution of ( 5). In so doing the expression of the azimuthal component of the potential vector is given in a general form as

A 1c (r, z) = +∞ 0 µN turn II (ar 1 , ar 2 ) e az J 1 (ar) 2a 3 (h 2 -h 1 )(r 2 -r 1 )
× e -ah1 -e -ah2 da

+ +∞ 0 C 1 e -az J 1 (ar)da, ∀z ∈ [0; +∞[ (7)
A 2 (r, z) = +∞ 0 [C 2 F 1 (f (z)) + B 2 F 2 (f (z))] J 1 (ar) da ∀z ∈ [-r; 0] (8)
A 3 (r, z) = +∞ 0 B 3 F 3 (g (z)) J 1 (ar) da ∀z ∈] -∞; -r] (9)
where f (z) and g (z) are two different functions of z and where I (x 1 , x 2 ) is given by

I (ar 1 , ar 2 ) = x2 x1 xJ 1 (x) dx. ( 10
)
Once A 1c (r, z) is known, Z is deduced as

Z = K +∞ 0 I (ar 1 , ar 2 ) 2 a 6 2 e -a(h2-h1) -1 + a(h 2 -h 1 )
+ e -ah2 -e -ah1 2 φ(a) da

(11) with φ(a) = C 1 K , (12)
K = µN turn II (ar 1 , ar 2 ) e -ah1 -e -ah2 2a 3 (h 2 -h 1 )(r 2 -r 1 ) (13)
K = jωπµN 2 turn (h 2 -h 1 ) 2 (r 2 -r 1 ) 2 . ( 14
)
The expressions of C 1 , C 2 , B 2 and B 3 are obtained from the continuity conditions of the quantities and/or their derivatives with respect to depth and/or their cancellation at ±∞. Denoting all such quantities as

L = F 1 (f (z = 0)) ; (15) M = F 2 (f (z = 0)) ; (16) N = F ′ 1 (f (z))| z=0 ; (17) O = F ′ 2 (f (z))| z=0 (18) P = F 1 (f (z = -r)) ; (19) Q = F 2 (f (z = -r)) (20) R = F ′ 1 (f (z))| z=-r ; (21) S = F ′ 2 (f (z))| z=-r (22) T = F 3 (g(z = -r)) ; (23) U = F ′ 3 (g(z))| z=-r (24)
where ′ means derivative with respect to z and solved the corresponding linear system, the expression of φ(a) is obtained as

φ(a) = (aM -O) (RT -P U ) + (aL -N ) (QU -ST ) (aM + O) (RT -P U ) + (aL + N ) (QU -ST ) (25)
Applying such an approach, analytical solutions of (5) for a continuous conductivity function σ(z) have been proposed when σ(z) varies linearly, quadratically or exponentially [START_REF] Theodoulidis | Analytical solutions in eddy current testing of layered metals with continuous conductivity profiles[END_REF] and when σ(z) is varying as an hyperbolic tangent [START_REF] Uzal | The impedance of eddy current probes above layered metals whose conductivity and permeability vary continuously[END_REF]. In [START_REF] Kolyshkin | Analytical solution in eddy-current testing of double-layer media with depthvarying magnetic properties[END_REF] the authors are interested in the case of double-layer media with depth-varying magnetic properties where the permeability is continuously varying as a function of the exponential of the depth whereas the case of a layered medium where both the conductivity and the permeability are continuously varying as a power of the depth is adressed in [START_REF] Kolyshkin | Analytical solutions to eddy-current testing problems for a layered medium with varying properties[END_REF].

B. Series expansion of the formulation (11)

According to the method described in [START_REF] Theodoulidis | Series expansions in eddy current nondestructive evaluation models[END_REF] the integral expression of the impedance [START_REF] Schmidt | Estimating the eddy-current modeling error[END_REF] can be replaced by a series expansion (26) which allows a better control of the convergence and a faster evaluation at least at low frequency. This method is based on a reduction of the solution region in the radial direction at a distance b where a homogeneous Dirichlet condition is imposed onto the radial component of magnetic field. The integrand of ( 11) is slightly modified and only evaluated at discrete values a i of a, the latter being obtained from J 1 (a i b) = 0. The number M of a i to be taken into account depends on the sought accuracy on the value of Z.

Z M ≈ 2K M i=1 I(a i r 1 , a i r 2 ) 2 [a i bJ 0 (a i b)] 2 a 5 i × e -aih2 -e -aih1 2 φ(a i ) + 2 e -ai(h2-h1) -1 + a i (h 2 -h 1 ) (26)
(note that Z M → Z when M → +∞ for sufficiently large b).

C. Particular formulation for a two-tanh-continous profile

For the applications which we are interested in, the inward and outward diffusion are to be taken into account, so we suggest to apply the analytical method for dealing with (5) with a conductivity depth profile defined by two hyperbolic tangents, one for the inward diffusion and the other for the outward one, as follows:

σ(z) =        σ 12 + ∆σ 1 + tanh z + c 1 2v 1 ∀z ∈ [-r, 0] σ 2 -∆σ 1 + tanh z + c 2 2v 2 ∀z < -r (27)
with ∆σ = (σ 1 -σ 12 ) /2 and where σ 1 is the conductivity limit value related to the first profile when z tends to +∞, σ 2 is the conductivity limit value related to the second profile when z tends to -∞ and σ 12 is the joining conductivity limit value related to the first profile when z tends to respectively -∞ and to the second profile when z tends to +∞.

The joining point between the two sub-profiles, denoted r, is obtained by equating expressions of σ(z) in ( 27) for z = -r. v 1 and v 2 can be seen as "the transition speed" of the conductivity in the [-r, 0] and [-∞, -r] and c 1 and c 2 the depth corresponding to the inflection point. A typical profile of such a kind is shown in Fig. 1. In fact, such a modeling allows us to define a large class of various profiles with a small number of parameters adequately chosen. Let us notice that with such a formulation a continuous increasing (or decreasing) profile cannot be easily obtained even if, at the boundary of the two layers, the conductivity, function of depth, is allowed to be discontinuous.

For the case of a single tanh-profile all details can be found in [START_REF] Uzal | The impedance of eddy current probes above layered metals whose conductivity and permeability vary continuously[END_REF] and in the following only the two tanh-profiles will be dealt with. In such a case it can be shown that the expression of the special functions F 1 , F 2 and F 3 are

F 1 (y 2 (z)) = y µ 2 (z) [1 -y 2 (z)] ν × 2 F 1 (µ + ν, µ + ν + 1; 2µ + 1; y 2 (z)) F 2 (y 2 (z)) = y -µ 2 (z) [1 -y 2 (z)] ν × 2 F 1 (ν -µ + 1, ν -µ; 1 -2µ; y 2 (z)) F 3 (y 3 (z)) = y λ 3 (z) [1 -y 3 (z)] τ × 2 F 1 (λ + τ, λ + τ + 1; 2λ + 1; y 3 (z)) (28) with y 2 (z) = 1 + e - z+c 1 v 1 -1 , y 3 (z) = 1 + e -z+c 2 v 2 -1 , µ = v 1 a 2 + jωµ 0 σ 12 , ν = v 1 a 2 + jωµ 0 σ 1 , λ = v 2 a 2 + jωµ 0 σ 2 , τ = v 2 a 2 + jωµ 0 σ 12 .
(29)

2 F 1 (α, β; γ; x) is the Gauss hypergeometric function defined as [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF] 2 F 1 (α, β; γ; x)

= +∞ k=0 (α) k (β) k (γ) k x k k! (30)
where (α) k is the Pochhammer symbol given by

(α) k = α(α + 1)(α + 2) • • • (α + k -1). (31)
Let us define the derivative of 2 F 1 (α, β; γ; x) with respect to

x as [15] ∂ ∂x 2 F 1 (α, β; γ; x) = αβ γ 2 F 1 (α + 1, β + 1; γ + 1; x) (32)
Introducing the expressions of (28) in (24) and using the derivative's rules (32) allows to express the analytical formulation of the functions L, M, N, O, P, Q, R, S, T, U (see § A for their analytical expressions) and then to obtain the expression of φ (a) through (25).

IV. BRIEF DESCRIPTION OF THE MULTI-LAYER MODEL

Proper results for numerically validating our approach in the two-tanh profile configuration do not appear in the literature. A classical scheme consisting in the approximation of the continuously varying conductivity profile by slicing the coating layer into a sufficiently high number of planar homogeneous layers of given thickness and constant conductivity is considered. It can be sketched as follows 1) the lower half-space (z < 0 in Fig. 1) is divided in two parts:

• the upper one (zone 2 in Fig. 1) is a nonhomogeneous slab of thickness d. It is divided into N homogeneous sub-layers of constant thickness ∆(= d/ (N -1)) and conductivity σ i , i = 1 . . . , N -1, such as to provide a good discretized description of σ (z); • the lower one (zone 3 in Fig. 1) is taken as a halfspace of constant conductivity σ (z = d); 2) the Fresnel reflection coefficient for a TE-polarized electromagnetic wave of the multi-layer structure is computed using the algorithm proposed in [START_REF] Cheng | General analysis of probe coils near stratified conductors[END_REF]. It can be shown (no detail is provided here) that such a coefficient is nothing but φ (a) in [START_REF] Schmidt | Estimating the eddy-current modeling error[END_REF] or in (26); 3) for each frequency f , a given N and a given d are chosen and the corresponding impedance Z is computed using either ( 11) or ( 26); This method should ensure, in theory, high accuracy of the solution. In practice, the high number of discretization layers requested leads to a method which might be unsuitable, despite its generality, due to increase of computing time with respect to the number of layers considered.

For our purpose and for sake of simplicity, let us consider a configuration for which we would like to compute the impedance Z for a range of frequencies from f min to f max . Then d is chosen large enough to entirely include the varying part of σ (z) and is kept constant for all frequencies. Here d has been arbitrarily chosen as:

d = 2 ω min µ 0 σ min ( 33
)
with σ min = min (σ 1 , σ 12 , σ 2 ) and ω min = 2πf min , f min being the smallest frequency which we are interested in. Once d has been fixed then N is chosen large enough by trial-error scheme so that the corresponding impedance Z N -obtained for N layers-tends to the true impedance Z at the higher frequency f max . Such a choice is clearly not optimal and could be efficiently improved.

1) The thickness ∆ i of the i th sub-layer should be adapted to the local curvature of the profile σ (z) -small when σ (z) is rapidly varying, large when σ (z) is almost constant. 2) Only the M first sub-layers, those which are able to correctly describe the attenuation phenomena, have to be taken into account. As a matter of fact, knowing that eddy currents only penetrate at few skin depths in a metal plate and that the skin depth varies with respect to frequency (see ( 3)) the higher the frequency the smaller the number of sub-layers to take into account. However for shake of simplicity such an adaptive scheme has not been applied in the following.

V. NUMERICAL TREATMENT

Numerically speaking the computation of ( 11) is difficult. As a matter of fact, even if we know that such an integral is convergent, -I(r 1 , r 2 ) is decreasing and φ (a) plays the role of a reflection coefficient and then is bounded by 1 when a goes to +∞the theoretical formulation of φ (a) (25) involves the hypergeometric function (30) or its derivative (32) in its numerator and denominator.

The fast computation of the Gauss hypergeometric function 2 F 1 with complex parameters is a difficult task and is still the subject of significant research as exemplified in [START_REF] Perger | A numerical evaluator for the generalized hypergeometric series[END_REF] and [START_REF] Michel | Fast computation of the Gauss hypergeometric function with all its parameters complex with application to the Pöschl-Teller-Ginocchio potential wave functions[END_REF] with two different approaches. Both of them provide Fortran subroutines which have been used in our numerical simulations, [START_REF] Perger | A numerical evaluator for the generalized hypergeometric series[END_REF] proposed a very general approach to compute all hypergeometric functions p F q (subroutine PFQ) whereas [START_REF] Michel | Fast computation of the Gauss hypergeometric function with all its parameters complex with application to the Pöschl-Teller-Ginocchio potential wave functions[END_REF] focuses on the computation of the Gauss hypergeometric function 2 F 1 (subroutine HYP 2F1) we are interested in.

The numerical integration of ( 11) is performed using the DQAGI subroutine provided by IMSL.

VI. NUMERICAL VALIDATION

A. A single tanh profile [START_REF] Uzal | The impedance of eddy current probes above layered metals whose conductivity and permeability vary continuously[END_REF] The validation of our approach is made by comparison with the results given in [START_REF] Uzal | The impedance of eddy current probes above layered metals whose conductivity and permeability vary continuously[END_REF] for a single tanh profile. With reference to Fig. 1 the parameters of the measurement configuration are taken as: r 1 = 1.3 mm, r 2 = 3.3 mm, h 1 = 0.5 mm, h 2 = 7.8 mm and N turn = 580. According to (27) the conductivity profile has the following parameters: σ 1 = 1.509 × 10 7 S m -1 , σ 12 = 3.766 × 10 7 S m -1 , c 1 = 0.3 mm and v 1 = 0.1857 mm.

A comparison of the results at three frequencies is shown in Table I for different values of b. A brief study of the convergence of the results with respect of the choice of the hyper-parameter b is presented. A very good agreement is observed with [START_REF] Uzal | The impedance of eddy current probes above layered metals whose conductivity and permeability vary continuously[END_REF]. To illustrate the ability of the series expansion to provide good results at a lower computional cost than the numerical integration a comparison of the CPU times needed to compute the impedance for 100 frequencies between 1 kHz to 100 kHz every 1 kHz is proposed. As exemplified in Table II the subroutine HYP 2F1 seems to be much more efficient than PFQ for our application.

B. A two-tanh profile case

As already said ( § IV), in the two-tanh profile configuration we have not been able to find proper results in the literature to numerically validate our theoretical approach. A multilayer model (described § IV) have then been developped and numerical results compared. Again with reference to Fig. 1 Let us define the discrepency between a reference ∆Z i , i = 1, . . . , N freq and an approximated value ∆Z i taken as

Err = 1 N freq Nfreq i=1 |∆Z i -∆Z i | |∆Z i | (34)
Comparison of variations of impedance obtained with the analytical model and those obtained with a numerical multi- [START_REF] Schmidt | Estimating the eddy-current modeling error[END_REF] WITH THE ANALYTICAL AND THE MULTI-LAYER MODELS. layer model using the numerical integration is very good as reported in Fig. 3, corresponding CPU time and error being given in Table 3. Let us emphasize that best matching with the multi-layer model requires a number of, at least, N = 100 000 layers, which can be explained by the fact that, as described in § IV, the two multi-layer model parameters d and N are kept constant for all frequencies and whatever the variation of the conductivity profile.

The analytical formulation and numerical integration being validated let us now focus on the series expansion approach and its advantages. In the following only the analytical formulation using the subroutine HYP 2F1 to compute 2 F 1 will be used and a comparison between the results using the numerical integration and the series expansion approach is provided in Fig. 4 and the corresponding CPU times and errors are reported in Table IV. It can be seen that b -which has been chosen as an integer multiple of R 2 as b = N R 2 -should be carefully chosen and a too small value, here N = 2, leads to a large error whereas a higher value, N = 10, gives very good results. As expected when N increases the error between numerical integration and the series expansion decreases and tends to 0 which is paid for by an increase in computing time even if the "series" CPU time is still ten times lower (for N = 20) than the "numerical integration" one.

VII. CONCLUSION

In the context of development of decentralized power generation, it is of economical importance to improve the thermal efficiency of gas turbine by increasing the firing tempera- ture. To optimize the maintenance cost of such advanced gas turbine, it is crucial to predict remaining coating life for reducing unexpected shut-downs caused by blade coating failure due to depletion of aluminum. Eddy Current Non-Destructive Evaluation can play a key role in assessing the conditions of turbine blade coating by estimating conductivity depth profiles. In this paper, we have reported on efforts towards that goal. We have developed the direct model of the impedance of an EC coil above a metal structure protected by a top over-aluminized coating degraded by inward and outward diffusion of aluminum, assuming that the conductivity profile follows a two-hyperbolic-tangent law. The analysis of the approach as suggested is quite interesting, since the results obtained with this model are similar to those obtained with a numerical multi-layer model with a significant reduction of computing time.

Fig. 1 .

 1 Fig. 1. Typical configuration of measurement with a theoretical conductivity profile with inward and outward diffusion of aluminum (left); example of a real diffusion process (right).

  008 28 -0.225 71 -1.497 19 Numerical Integration -0.008 28 -0.225 65 -1.496 37 N = 10 -0.008 267 -0.225 57 -1.496 45 N = 20 -0.008 28 -0.225 66 -1.496 77

  the parameters of the measurement configuration are: r 1 = 2 mm, r 2 = 4 mm, h 1 = 0.5 mm h 2 = 7.3 mm and N turn = 200. According to (27) the conductivity profile has the following parameters: σ 1 = 7 × 10 5 S m -1 , σ 12 = 5 × 10 5 S m -1 , σ 2 = 8 × 10 5 S m -1 , c 1 = 0.2 mm, v 1 = 0.03 mm, c 2 = 0.8 mm and v 2 = 0.1 mm. The variation of impedance ∆Z i is computed at N freq frequencies equally distributed between f min = 100 kHz and f max = 10 MHz every 100 kHz.

  Fig. 2. Conductivity Profile.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Comparison of results obtained with analytical and multi-layer models for different number of layers N .

TABLE I COMPARISON

 I OF THE RESULTS FROM [5] AND THE ONES OBTAINED USING THE NUMERICAL INTEGRATION AND THE SERIES EXPANSION WITH b = N R 2 , N = 10, 20 USING HYP 2F1.

		Numerical	Series expansion
		integration	N = 10	N = 20
	PFQ	-	37.41	66.22
	HYP 2F1	1.36	0.06	0.11

TABLE II COMPARISON

 II OF THE CPU TIMES (IN SECOND) NEEDED TO COMPUTE Z FOR 100 FREQUENCIES BETWEEN 1 kHz TO 100 kHz FOR THE TWO FORMULATIONS (INTEGRAL EQUATION (11) AND SERIES EXPANSION (26)) USING DIFFERENT SUBROUTINES TO COMPUTE 2 F 1 .

TABLE III COMPARISON

 III OF THE CPU TIMES (IN SECOND) AND THE ERROR (IN %) FOR THE COMPUTATION OF Z USING A NUMERICAL INTEGRATION OF

ACKNOWLEDGEMENTS

The authors would like to thank Dominique Lesselier and Theodoros Theodoulidis and the two anonymous reviewers for their helpful comments.

APPENDIX

In the following and in the case of the two tanh-profil the analytical expressions of the general functions presented in (24) are given.

+ [1 -y 3 (-r)] τ y -λ 3 (-r) (λ + τ )(λ + τ + 1) 2λ + 1 × 2 F 1 (λ + τ + 1, λ + τ + 2; 2λ + 2; y 3 (-r))

µ, ν, y 2 , y 3 are defined in (29).