Rafael Amaya-Gómez

Franck Schoefs

Mauricio Sánchez-Silva

Felipe Muñoz

Emilio Bastidas-Arteaga

Matching of corroded defects in onshore pipelines based on In-Line Inspections and Voronoi partitions

Keywords: Corroded pipeline, New defects, matching, Probability of detection, Voronoi cells

Onshore pipelines are usually subjected to a corrosion attack. Regular inspections known as In-Line inspections (ILI) are commonly used with magnetic (MFL) or ultrasonic (UT) tools to prevent any failure. New defects will appear between consecutive inspections due to the aggressiveness of the surroundings and the detection thresholds associated with the defects' depth. This work focuses on the matching problem between two inspections, aiming to identify the degradation increments and the position of new defects. Typically, it is linked to the well-known point matching problem in pattern recognition, where the objective is finding the best affine transformation between two sets of points in a plane. This work presents an alternative using Voronoi cells to filter possible matches and an iterative approach to determine the best affine transformation, considering the uncertainty in any direction. The approach was implemented for a real pipeline 45 km long and for synthetic corrosion defects, allowing us to identify possible matches easily. Based on the new and old defects, some insights about the probability of detection and false alarm are deduced. For this purpose, experimental probability and results from recognized exponential and log-logistic functions were considered.

Introduction

Onshore pipelines are critical components in the chain of production and distribution of oil and gas. These pipelines are generally subjected to extreme operation and environmental conditions that would lead to significant corrosion failure risks [START_REF] Zhang | Bayesian network model for buried gas pipeline failure analysis caused by corrosion and external interference[END_REF][START_REF] Taleb-Berrouane | Corrosion risk assessment using adaptive bow-tie (ABT) analysis[END_REF][START_REF] Adumene | Dynamic risk analysis of marine and offshore systems suffering microbial induced stochastic degradation[END_REF][START_REF] Liu | Dynamic risk assessment model of buried gas pipelines based on system dynamics[END_REF]. Therefore, condition assessment though inspections is crucial to ensure adequate levels of serviceability and safety. Corroded pipelines are commonly inspected with Magnetic (MFL) or ultrasonic (UT) tools following In-Line Inspection (ILI). These inspections provide the location and the extent of the metal loss detected based on the pipeline abscissa (longitudinal direction) and a 12-hour clock analogy for the circumferential position. This information is used to monitor the evolution of the corrosion attack every 2 to 6 years to support future maintenance and repair decisions. The challenge lies in the fact that these measurements are subjected to different uncertainties associated with the existence, position, and magnitude of the corrosion defects; for instance, ILI tools have a detection threshold of the depth of the defects near 10% of the wall thickness. Although ILI tools include this threshold, the detection problem does not follow a binary output where every defect above it would be indeed detected. On the contrary, there is a detection error in terms of the probability of being detected (PoD) and being a false alarm (PFA) [START_REF] Rouhan | Probabilistic modeling of inspection results for offshore structures[END_REF]. These probabilities are intrinsically correlated with the reported accuracy for the defect depth measurements, where it follows that the depth of the defects reported by ILI (d ILI ) is related to the real depth (d real ) by d ILI = d real ± 0.1t at 80% certainty.

Other uncertainties include shifted positions and the appearance of new defects possibly caused by a noise signal, complicating any direct matching between consecutive inspections to monitor the pipeline's condition. These new defects can be associated, on the one hand, with the aggressiveness of the surrounding soil, fluid operation, and fluid properties [START_REF] Amaya-Gómez | Statistical Soil Characterization of an Underground Corroded Pipeline Using In-Line Inspections[END_REF]. Several authors agreed that corrosion is favored when the resistivity is reduced regarding soil aggressiveness. Besides, it can be affected by the soil porosity, and groundwater conductivity [START_REF] Arriba-Rodriguez | Methods to Evaluate Corrosion in Buried Steel Structures: A Review[END_REF]. Whereas the flow velocity, Zhang and co-workers have obtained that a low flow velocity (< 2 m/s) may lead to a competition between the corrosive species and the inhibitor, and at high velocities (> 6 m/s), an erosion-corrosion may be favored [START_REF] Zhang | Effect of flow velocity on pipeline steel corrosion behaviour in H2S/CO2 environment with sulphur deposition[END_REF]. On the other hand, the new corrosion defects can appear due to a miss-detection in the first inspection or a false alarm in the second one. The reason why new defects appear is unknown without field measurements, but there is certainly an "indication" of a new corrosion point based on the reported in the ILI measurement.

Matching these corrosion defects is associated with the well-known point-matching problem for pattern recognition. This problem seeks an affine transformation using translation, scaling, and rotation parameters to match the highest number of pairs between the two inspection results. There are some available alternatives in this regard for pattern recognition, such as the Iterative Closest Point (ICP) of Besl & McKay [9] or the Thin Plate Spline Robust Matching (TPS-RPM) of Chui & Rangarajan [START_REF] Chui | A new point matching algorithm for non-rigid registration[END_REF]. For corroded pipelines, Dann & Dann [START_REF] Dann | Automated matching of pipeline corrosion features from in-line inspection data[END_REF] proposed a modified version of the TPS-RPM approach, considering possible displacements of the corrosion between inspections based on an Annealing algorithm. This approach seeks an affine transformation matrix with defined rotation, translation, and scaling factors. Considering the uncertainty in the results of the inspection, regarding the longitudinal and circumferential location of corrosion defects, it might be possible that the final matching may not be completely described by this transformation. This work presents an alternative framework to identify which points are classified as new and old defects using Voronoi partitions to filter preliminary matches and an iterative matching transformation. This framework is applied to a real pipeline of 45 km and synthetic corrosion defects. The approach incorporates location uncertainties and correspondence analysis with outlier filtering. The results in the real and synthetic measurements highlight this framework to be considered by pipeline operators in further condition-based monitoring. The result also allows us to discuss some insights into how this probability can be estimated using the depths of the new defects.

The document is structured as follows: Section 2 describes the point-matching problem and the main matching alternatives for pattern recognition. Section 3 introduces the matching alternative for the corroded pipeline and some drawbacks. Section 4 presents the proposed matching approach using a Voronoi partition. Section 5 characterizes the spatial dependencies of the real case study. Section 6 discusses the matching results for the real and synthetic measurements. Section 7 presents some insights about the probability of detection and false alarm, and Section 8 gives the final conclusions.

Point matching problem: Identifying new defects between inspections

General problem formulation

The point matching problem considers two sets of points P = {p i ∈ W, i = 1, 2, . . . , n} and Q = {q k ∈ W, k = 1, 2, . . . , m}, which represent the defects from the first and second ILI measurement, respectively. In this case, p i and q k are vectors containing the location of the pipe and W = (A × P) ⊂ R 2 is the plane produced by the abscissa and clock-position once the pipeline is unrolled axially. The traditional point pattern problem search for an affine transformation T for p i := (x p i , y p i ) given by:

T x p i y p i = t x t y + S cos θ -sin θ sin θ cos θ x p i y p i , ( 1
)
where S is a scaling factor, θ is a rotating angle and t x and t y are translations in both main directions. Figure 1 illustrates this matching problem. According to Van Wamelen et al. [START_REF] Van Wamelen | A fast expected time algorithm for the 2-D point pattern matching problem[END_REF], a matching between the two sets of observations would take place, given a matching probability δ ∈ [0, 1] and a matching size γ ∈ R, if there exists a subset of δn points P M of P such that for each p ∈ P M it follows that T (p)q < γ for a q ∈ Q.

Matching alternatives

Point pattern matching is a problem that has been studied before in the field of pattern recognition through different approaches, including fuzzy relaxation, asymmetric neural networks, 2D clustering, and nearest neighbor search approaches [START_REF] Van Wamelen | A fast expected time algorithm for the 2-D point pattern matching problem[END_REF]. Most of these approaches implement a filter approach to determine the number of ≤ min(n, m) of possible matches. One alternative is the Iterative Closest Point (ICP) algorithm that uses the set of nearest neighbors of the points in P from the points Q, here denoted as O := {o 1 , . . . , o n }, in a quaternion-based algorithm that implements an optimal translation qT and rotation qR = [q o , q 1 , q 2 , q 3 ] t transformations [START_REF] Besl | A method for registration of 3-D shapes[END_REF]. The latter is associated with a unit quaternion (i.e., q 2 o + q 2 1 + q 2 2 + q 2 3 = 1). The objective is to minimize the following function:

f ( q) = 1 n n i=1 o i -R( qR )p i -qT 2
where q = [ qR qT ] and R is the rotation matrix that is produced by the unit quaternion rotation as follows:

R =           q 2 0 + q 2 1 -q 2 2 -q 2 3
2(q 1 q 2q 0 q 3 ) 2(q 1 q 3 + q 0 q 2 ) 2(q 1 q 2 + q 0 q 3 ) q 2 0 + q 2 2q 2 1q 2 3 2(q 2 q 3q 0 q 1 ) 2(q 1 q 3q 0 q 2 ) 2(q 2 q 3 + q 0 q 1 ) q

2 0 + q 2 3 -q 2 1 -q 2 2          
.

Following the closed-solution reported by Horn [START_REF] Horn | Closed-form solution of absolute orientation using unit quaternions[END_REF], the optimal rotation transformation qR is determined from the eigenvector of the maximum eigenvalue of the symmetrical matrix whose elements are combinations of sums of products of the point coordinates. Optimal translation qT corresponds with the difference of the centroid of points in Q (center of mass) and the rotated centroid of points in P. This process follows iteratively until the difference between iterations of the mean square error of the point matching converges to a predefined threshold [START_REF] Besl | A method for registration of 3-D shapes[END_REF], which can be linked to the accuracy of the location.

Other alternatives consider not only the nearest neighbor, but also the k nearest neighbors of each point in both sets. For instance, the approach reported by Van Wamelen et al. [START_REF] Van Wamelen | A fast expected time algorithm for the 2-D point pattern matching problem[END_REF] compares the k nearest neighbors of each point p i with the points in q k and determines its local affine transformation with the approach reported by Chang et al. [START_REF] Chang | Fast algorithm for point pattern matching: Invariant to translations, rotations and scale changes[END_REF]. Chang and co-workers proved that if there are pairs ( ≥ 2) of possible matching points a i = (x a i , y a i ) ↔ b i = (x b i , y b i ) for i = 1, . . . , , then the affine transformation T that gives the best least-square match that minimizes i T (a i )b i 2 , which is denoted by r = (t x , t y , S cos θ, S sin θ) t , is given by [START_REF] Chang | Fast algorithm for point pattern matching: Invariant to translations, rotations and scale changes[END_REF]:

r = 1 det               l A 0 -µ X A µ Y A 0 l A -µ Y A -µ X A -µ X A -µ Y A 0 µ Y A -µ X A 0                             µ x B µ Y B l A+B l A-B               , (2)
where

µ X A = i x a i , µ X B = i x b i , µ Y A = i y a i , µ Y B = i y b i , (3)
l A+B = i (x a i x b i + y a i y b i ), l A-B = i (x a i y b i + y a i x b i ), l A = i (x 2 a i + y 2 a i ), (4)
det = • l A -µ 2 X A -µ 2 Y A . ( 5
)
This transformation is used by Van Wamelen et al. [START_REF] Van Wamelen | A fast expected time algorithm for the 2-D point pattern matching problem[END_REF] considering as possible matches those points with T (a i )b k < γ = R Q /n for some a i and b k neighbors of p i and q k , with R Q being the minimum disk radius that covers the points in Q. This process continues until a higher number of matches is achieved.

Finally, some approaches implement correspondence matrices with the possibility of outliers from both sets. For instance, the Thin Plate Spline Robust Matching (TPS-RPM) reported by Chui & Rangarajan [START_REF] Chui | A new point matching algorithm for non-rigid registration[END_REF], and the revised version of Yang [START_REF] Yang | The thin plate spline robust point matching (TPS-RPM) algorithm: A revisit[END_REF] using a double-sided approach to manage the outliers under a reduced optimization approach, which use an annealing temperature process to solve the optimization problem. Consider a correspondence matrix C = [c ik ] with dimensions n×m that indicates the probability of matching each p i ∈ P with each q k ∈ Q. This approach includes two possible outliers vectors r and s with dimensions (n × 1) and (1 × m), respectively, associated with the defects that could not be matched in the optimization problem shown in Eq. 6. In this objective function, α and λ are controlling parameters to avoid classifying every point as an outlier (later used in the annealing heuristic), LT is the thin-plate splines smooth regularization, and T is the temperature parameter of the annealing process beginning at T o and decreasing in each iteration.

argmin

C,r,s,T n i=1 m k=1 c ik q k -T (p i ) 2 + λLT -α n i=1 m k=1 c ik + T m k=1 n i=1 c ik log c ik subjected to 0 ≤ c ik , r i , s k ≤ 1, ∀i = 1, . . . , n, ∀k = 1, . . . , m m k=1 c ik + r i = 1, ∀i = 1, . . . , n n i=1 c ik + s k = 1, ∀k = 1, . . . , m (6)
The annealing process helps to solve this optimization problem by approximating the values of the correspondence matrix and the outlier vectors, given a fixed transformation T following Eq 7 to 9. Here, u c and v c are the centers of mass of both sets. The approximated version of c ik ,s k , and r i are finally implemented to obtain the optimal affine transformation.

c ik = 1 T exp α T - (q k -T (p i )) t (q k -T (p i )) T , ∀i = 1, . . . , n, ∀k = 1, . . . , m (7)
r i = 1 T o exp - (u c -T (p i )) t (u c -T (p i )) T o , ∀i = 1, . . . , n (8)
s k = 1 T o exp - (q k -v c ) t (q k -v c ) T o , ∀k = 1, . . . , m (9)

Matching alternatives in corrosion pipelines: main drawbacks

The alternatives shown in the previous section address the general problem of point matching, which usually deals with image recognition like fingerprints or distortion for 3-D shapes; see, for instance, [START_REF] Besl | A method for registration of 3-D shapes[END_REF][START_REF] Yang | The thin plate spline robust point matching (TPS-RPM) algorithm: A revisit[END_REF]. This problem has also been reported for matching defects in consecutive inspections as an alternative of matching learning approaches using Support Vector Machines (SVM), decision trees, and random forest [START_REF] Liu | Matching pipeline In-line inspection data for corrosion characterization[END_REF]. In this direction, Dann & Dann [START_REF] Dann | Automated matching of pipeline corrosion features from in-line inspection data[END_REF] adopted the approach of Yang [START_REF] Yang | The thin plate spline robust point matching (TPS-RPM) algorithm: A revisit[END_REF] of the correspondence matrix and outlier vectors in Eq. 6.

The approach of Dann & Dann divides the pipeline into fixed segments that overlap every 0.3 m at the joints, using the points mapped into a 2D plane by unrolling the pipeline in the axial direction. The plane uses an extended version by duplicating half segments on each side, as shown in Fig. 2 to account for the pipeline continuity in the circumferential direction and to prevent any biased misclassification. This figure illustrates how those defects located from 06:00 hr to 12:00 hr (and from 00:00 to 06:00) duplicate at -06:00 hr to 00:00 hr (at 12:00 hr to 18:00 hr, respectively). Dann & Dann used an alternative version of Eq. 1 by projecting the defects in both sets to the plane z = 1, i.e., a = (x a , y a , 1) and b = (x b , y b , 1) for any a ∈ P and b ∈ Q, and using a 3 × 3 affine transformation T 3 , whose components t 11 , t 12 , t 21 , t 22 are associated with the scaling and rotation, t 13 , t 23 with the translation, and the last row vector is [0, 0, 1] just to maintain the points in the same plane.

The optimization problem acknowledges the duplicates of the defects by augmenting the matrix c ik and the vector r i to (2n × m) and (2n × 1), respectively. This approach seeks for affine transformations that are not very different from I 3 (the identity 3 × 3 matrix) by using a controlling parameter β c as follows β c Tr (T 3 -I 3 ) t (T 3 -I 3 ) , where Tr[A] is the trace of matrix A. The final optimization problem is shown in Eq. 10, which is solved using again the annealing process with the temperature decreasing by T = T o γ l with γ ≤ 1 and l being the iteration step.

This formulation seems adequate for the matching of consecutive inspections given that it considers the possibility of outliers in both sets, it seeks for a probability of matching (known as Soft-assign), and it prevents the trivial solution of [c ik ] = 0 by including the term α 2n i=1 m k=1 c ik . The matching approach depends on the selection of α and β c , where α could be linked to the outliers proportion in both sets, and β c is a penalization parameter associated with the transformation deviation [START_REF] Dann | Automated matching of pipeline corrosion features from in-line inspection data[END_REF]. This approach assumes that both parameters are known, and the authors proposed the following feasible values: α ∈ [0, 0.1] and β c ≥ 200T . However, it is noteworthy that possible matches can be misclassified as outliers by forcing a final transformation to be nearly the same as the identity matrix or assigning a lower correspondence. According to Dann & Dann, this correspondence should be greater than 0.9 to classify as a positive match. For illustrative purposes, consider the same example used in the work of Dann & Dann, and also the subset of pairs (p i , q k ) with p i ∈ P and q k ∈ Q such that the nearest neighbor of p i is q k and vice versa. Figure 3 depicts the corresponding results based on the approach of Dann & Dann with annealing, where ILI 1 and ILI 2 depict the position from the two consecutive inspections and ILI 1 mod is the results after applying the annealing method. Note that the approach achieves a clear and good matching result when there is a marked affine transformation (Figure 3a), in this case, a diagonal displacement in the upper right direction. However, Figure 3b shows a clear case when the transformation is barely different from the identity matrix producing a significant misclassification. The correspondence and outliers results, in this case, indicate that only 7 over 21 defects could be matched, while in the previous, case almost all of them could be paired. The point matching problem for corroding pipelines goes beyond only the affine transformation. The physical explanation of this mismatching comes from the fact that the reported defects are subjected to location uncertainties in both the axial and circumferential directions, which can be thought of as a calibrating uncertainty of the defect's position, as remarked by Pakrashi et al. [START_REF] Pakrashi | ROC dependent event isolation method for image processing based assessment of corroded harbour structures[END_REF]. Let (x, y) denote a defect, this defect would actually lie in the region (x ± δ x , y ± δ y ), where δ x and δ y depend on the inspection tool. For instance, following the reported by an inspection vendor [START_REF] Rosen | Magnetic Flux Leakage[END_REF], it can be assumed that δ x = 0.1m and δ y = 10 • (i.e., π D 36 , where D is the pipeline diameter). In this regard, two possible criteria can be used to identify possible matches using the region formed by (x ± δ x , y ± δ y ) or the ball centered at each defect with radius δ 2 x + δ 2 y , as shown in Figure 4. Note that ball radius can be computed only when δ x and δ y have the same units. Overall, the matching problem of corroding pipelines should also consider that some defects already have good pairing, while others relate to their nearest neighbors. Figure 5 illustrates this case using a real segment in the inner wall once the approach of Dann & Dann is applied again. This figure shows some defects with suitable preliminary matches (e.g., between km 5654 and 5655) that produce transformation with very low correspondence probabilities.

This matching may be improved by omitting already detected matches, but this preprocessing would require significant manual feature selection that seeks to be avoided.

Proposed approach using a Voronoi partition

This work proposes a simple approach to match the corrosion defects based on the nearest neighbor criterion between two inspections, including the possibility of outliers. For this purpose, this approach uses Voronoi cells to filter possible preliminary matches between the two inspections. Given a set of points

P = {p 1 , . . . , p n } in R 2 , the
Voronoi cell of p l (V l ) is defined by the subset of points in R 2 such that the nearest point in P is p l . For a more formal

definition consider the Euclidean metric in R 2 d(p, q) = (x p -x q ) 2 + (y p -y q ) 2 , then V l = {x ∈ R 2 |d(x, p l ) ≤ d(x, p i ), ∀i l}.
The use of the Voronoi cells allows identifying the nearest neighbors between the points of P and Q quickly, as is illustrated in Figure 6. This figure depicts the Voronoi polygons given a query, which is the closest point of P or Q. Based on the mentioned above, the approach proposes three stages to identify possible matches between two consecutive inspections (Figure 7). First, defects from two consecutive ILI measurements are preprocessed by determining their crossed nearest neighbors. Second, these nearest neighbors are implemented in an iterative process to estimate matching transformation until the transformed data converges in each iteration. Finally, the transformation determines the correspondence matrix and outliers vectors. This approach will be described in detail below1 .

Figure 7: Proposed methodology for defects matching

Data processing

The data processing starts by duplicating each defect of the first inspection (ILI 1), as depicted in Figure 2.

This duplication aims to account for the pipeline continuity in the circumferential direction, i.e., ensuring that those points at 12:00 hr continue in the 0:00 hr. Denote the number of final points in the first inspection by 2n and m

for the second. The Voronoi cells V 1 1 , . . . , V 1 2n and V 2 1 , . . . , V 2
m are determined iteratively following any available method. Some of the common methods include the Plane-sweep, the tree expansion & deletion algorithm, or iterative approaches evaluating the dual diagram, i.e., the Delaunay Triangulation [START_REF] Dobrin | A review of properties and variations of Voronoi diagrams[END_REF][START_REF] Lee | Two algorithms for constructing a Delaunay triangulation[END_REF]. This work considers the latter approach based on the functions included in the R-project package deldir for Delaunay Triangulation and Dirichlet (Voronoi) Tessellation. These Voronoi cells are required to be assigned to a unique tile list from 0 to 2n (0 to m, respectively) for the first (second) inspection to identify the tiles containing each mixed point, i.e., points sharing being the nearest neighbors, as illustrated in Figure 8. Denote by T l1 and T l2 the tile list of the mixed points, and

T V 1
and T V 2 as the 2n × m matrices version of these lists, as the form of an indicator matrix. For instance, if the first entry of T l1 is s then the first row of T V 2 = e s , where e s is the standard unit vector of m positions with 1 in the s th position and elsewhere zero. Similarly, if the first entry of T l2 is r, then the first column of T V 1 = e r with 1 in the r th position.

The preliminary mixed nearest neighbors are determined by the Hadamard product or the element-wise product, i.e.,

T V := (T V 1 T V 2 ) i j = (T V 1 ) i (T V 2
) j , and by considering only those pairs of defects where d(p i , q k ) < δ where δ is defined by δ = δ 2 x + δ 2 y or by the vector (δ x , δ y ) as illustrated in Figure 4. Denote the set of nearest neighbors by (p i , q i ) for i = 1 : ≤ min(2n, m), and select from the original points or the duplicate those with the lower nearest neighbor distance.

Matching transformation approach

The second stage seeks the best feasible transformation given the initial mixed nearest neighbors obtained from both Voronoi cells. For this purpose, an iterative method is implemented following a similar procedure as reported by

Besl & McKay in the Iterative Closest Point (ICP) [START_REF] Besl | A method for registration of 3-D shapes[END_REF]. First, the nearest neighbors are determined using the Voronoi cells, then the best affine transformation is determined following the optimal transformation reported by Chang et al. [START_REF] Chang | Fast algorithm for point pattern matching: Invariant to translations, rotations and scale changes[END_REF] and shown in Eq. 2. This transformation is applied to the data from the first inspection, and the iterative process finishes when the mean square error (MSE) of the matchings between two consecutive iterations is shorter to a predefined threshold, i.e., if d k = i=1 p iq i 2 for the k th iteration, then the stopping criterion is achieved when

τ > min(d k -d k+1 , d k+1 ).

Correspondence and outliers optimization approach

The final matches are identified following the correspondence matrix and vector outliers. Given the affine transformation T 3 in Eq. 10, the optimization problem turns into linear programming with linear constraints, as shown in Eq.

11, where w ik = q k -T 3 (p i ) 2 , which may be limited to a binary solution indicating the matches and the classification as outliers in both sets [START_REF] Yang | The thin plate spline robust point matching (TPS-RPM) algorithm: A revisit[END_REF]. Several methods can solve this linear optimization problem, but this approach uses linear programming following the package lpSolve of R-project due to its simplicity. Other approaches can be studied like a Benders Decomposition Method (BDM), which is an approach with an increasing interest in the field of optimization for both Linear (LP) and Mixed-Integer Linear Programming (MILP) problems; however, it is out of the scope of this work.

argmin

C,r,s 2n i=1 m k=1 c ik (w ik -α) subjected to 0 ≤ c ik , r i , s k ≤ 1, ∀i = 1, . . . , 2n, ∀k = 1, . . . , m m k=1 c ik + r i = 1, ∀i = 1, . . . , 2n 2n i=1 c ik + s k = 1, ∀k = 1, . . . , m (11)
5. Spatial dependencies of the case study

Main parameters

The case study concerns an API 5LX52 pipeline 45km long with six main vales; its height lies between 2560 to 2660m above sea level. Near kilometer 33, there is a river crossing, whereas the last 10km are close to urban zones.

The pipeline has a nominal wall thickness of 6.35mm and an external diameter of 273.1mm. The analysis presented here was based on data obtained from two consecutive ILI measurements two years apart. According to the ILI report, this diameter is maintained along the entire abscissa, while the wall thickness exhibits greater variability due to the location of welded covers, valves, dents, and manufacturing flaws. For further details, please refer to Amaya-Gómez et al. [START_REF] Amaya-Gómez | A condition-based dynamic segmentation of large systems using a Changepoints algorithm: A corroding pipeline case[END_REF]. Table 1 shows a broad classification of the soil along the pipeline following the taxonomy of the USDA (United States Department of Agriculture). The pipeline has a bituminous coating of coal tar and an impressed current cathodic protection (ICCP) system. The defects measuring tool was a Magnetic Flux Leakage (MFL). Based on information reported in Amaya-Gómez et al. [START_REF] Amaya-Gómez | Pattern recognition techniques implementation on data from In-Line Inspection (ILI)[END_REF] about the inspection vendor, it can be assumed a circumferential uncertainty of 5 • during the inspection. The measurement uncertainties of the defect depth, length, and width are given by d ILI = d real ± d , l ILI = l real ± l , and w ILI = w real ± w , where d ILI , l ILI , w ILI stand for the depth, length, and width reported by the ILI tool, and d , l , w are the measurement errors. The measurement errors can be assumed to follow normal distributions centered at 0 with standard deviations obtained from the inspection vendors [START_REF] Pandey | Estimation of parameters of degradation growth rate distribution from noisy measurement data[END_REF]. It is reasonable to assume that these standard deviations are as follows: σ d = 0.1 t with t the nominal wall thickness, σ l = σ w = 11.70mm, considering a length and width accuracy of 15mm with a confidence of 80% of the data. For confidential agreements, further details of the case study cannot be provided.

Main descriptors of corrosion defects

A summary statistics of these data sets is depicted in Table 2 with the mean value and the coefficient of variation (CoV) of all the reported defects per inspection and pipe wall. Because further information about defects shape is not available in ILI, the maximum rather than the average depth for each defect will be considered from now on. Note that Table 2 show a higher CoV for the outer wall than in the inner wall that can be explained by the higher variability of corrosion factors in the soil than inside in the pipe [START_REF] Amaya-Gómez | Statistical Soil Characterization of an Underground Corroded Pipeline Using In-Line Inspections[END_REF]. 6. Matching results and New defects identification

Matching implementation example

For a better understanding of the matching approach, the segment shown in Figure 6 is used as an example to illustrate how the matches are determined and what are the "farther matches". As mentioned in Section 4, data processing starts determining the Voronoi cells for each of the corrosion points in both inspections. Consider the example shown in Figure 6 with n = 13 and m = 29, the number of points in sets P and Q. Figure 9 depicts this example, including Voronoi tile identification from 1 to 13, for the first inspection and 1 to 29 for the second. In this example, duplication is unnecessary, considering that none of the points in ILI1 were located near the border. The tile list of the mixed points, i.e., the nearest neighbor in the other inspection, for the first thirteen positions are given as follows:

T l1 = {1, 3,
Consequently, the 13 × 29 matrix version of these lists, which are determined using row and column standard unit vectors (i.e., e i = [0, . . . , 1, . . . , 0] at i th position) follows:

T V 1 =                    , T V 2 =                       e 1 e 3
e 1 e 2 . . .

                     
Note that the each unit vector e i in T V 1 have dimension 1 × 13, whereas each unit vector in T V 2 29 × 1. The mixed nearest neighbors are determined by the element-wise product of T V 1 and T V 2 with their non-zero elements. In this example, the sparse version of this matrix, i.e., the non-zero positions, are given by T V = {(2,3), (3,1), (4,2), (5,4), [START_REF] Amaya-Gómez | Statistical Soil Characterization of an Underground Corroded Pipeline Using In-Line Inspections[END_REF][START_REF] Rouhan | Probabilistic modeling of inspection results for offshore structures[END_REF], [START_REF] Zhang | Effect of flow velocity on pipeline steel corrosion behaviour in H2S/CO2 environment with sulphur deposition[END_REF][START_REF] Besl | A method for registration of 3-D shapes[END_REF], [START_REF] Dann | Automated matching of pipeline corrosion features from in-line inspection data[END_REF][START_REF] Chui | A new point matching algorithm for non-rigid registration[END_REF], [START_REF] Van Wamelen | A fast expected time algorithm for the 2-D point pattern matching problem[END_REF][START_REF] Yang | The thin plate spline robust point matching (TPS-RPM) algorithm: A revisit[END_REF], [START_REF] Horn | Closed-form solution of absolute orientation using unit quaternions[END_REF][START_REF] Amaya-Gómez | A condition-based dynamic segmentation of large systems using a Changepoints algorithm: A corroding pipeline case[END_REF]}, where the first position is associated with the points in ILI 1 and the second in ILI 2 that are mixed nearest neighbors. For instance, note in Figure 9 that the seventh and eighth points of ILI 1 are the closest points to the ninth point in ILI 2, but the latter is the nearest neighbor, which leads to the mixed point [START_REF] Zhang | Effect of flow velocity on pipeline steel corrosion behaviour in H2S/CO2 environment with sulphur deposition[END_REF][START_REF] Besl | A method for registration of 3-D shapes[END_REF].

These nine pairs are initially used in the optimal transformation reported by Chang et al. [START_REF] Chang | Fast algorithm for point pattern matching: Invariant to translations, rotations and scale changes[END_REF] and shown in Eq. 2.

Afterward, the mixed nearest neighbors and the optimal transformation are determined iteratively until the stopping

criterion of τ = 0.001 > min(d k -d k+1 , d k+1
) is achieved, where recall that d k is the matching MSE for the k th iteration.

Considering a separation greater than δ = δ 2 x + δ 2 y = 0.11m, the final results correspond with all the mixed nearest neighbors reported at T V , except from the (2,3) pair.

Figure 10 shows the obtained results for the complete segment of seven meters long that was previously used for the annealing method (see Figure 5). This figure includes some additional "farther matches", defined as pairs with a separation greater than δ. In this case, Figure 10a depicts the complete segment shown in Figure 5 with the 21 and 39 defects reported for the first and second inspection, respectively. Figure 10b shows eleven matches determined with the proposed approach after transforming the data, which contrasts with the zero feasible matches in the annealing case. Finally, Figure 10c includes three possible additional matches with a larger separation, which should be further monitored. After the correspondence and outlier approach is implemented (see Eq. 11) with α = 0.01, the matches are confirmed, and the farther matches are considered as outliers. This parameter was considered because it obtained a higher true matching ratio for different synthetic samples.

Comparison of the proposed matching approach with the Annealing method

Section 4 proposed an alternative matching method based on the Voronoi tessellation and an iterative affine transformation. This approach raised as an alternative of existing methods like the Annealing method proposed by Dann & Dann [START_REF] Dann | Automated matching of pipeline corrosion features from in-line inspection data[END_REF]. This section compares both approaches using synthetic datasets. Let n = m = 30 and consider P := {(p xi , p yi )} 30 i=1 to be uniform random sampled in [0, 1] 2 and Q := {(p xi + ∆ x , p yi + ∆ y )} 30 i=1 where ∆ x and ∆ y are uniform random variables associated with the location uncertainties δ x = 0.1 and δ y = 10

• , i.e., ∆ x ∼ Uni f (-δ x , δ x )
and ∆ y ∼ Uni f (-δ y , δ y ). Consider 50 random P and Q datasets to compare both approaches using the true matching ratio. Considering that the annealing process depends on α, β c and λ (see Section 3), a previous sensibility analysis was implemented to determine which parameters achieve the highest prediction (Table 3). The results in Table 3 indicate that the matching ratio increases for α ≥ 0.01, while β c and λ do not affect the matching ratio drastically, but β c ∼ 400T and λ = 0.9 can be considered. For the proposed approach, an α ≥ 0.01 also obtained the highest true ratio, as is depicted in Figure 11. Based on the above, the correct matching ratio for the annealing process using these parameters and the proposed approach are compared in Figure 12. This figure includes the results with the annealing process with the parameters The results indicate that the proposed approach outperforms the annealing matching almost in every case. The true mean matching for the proposed approach was 0.848 against 0.607 for the annealing process. Figure 12 also displays the true matching ratio for the solver KNITRO, which is a solver specialized for non-linear optimization problems.

For this case, the optimization problem in Eq. 10 was solved using the NEOS server under an AMPL language with the parameters obtained from the sensitivity analysis. The results show similar results to those in the KNITRO solver.

However, it should be highlighted that the Voronoi approach does not require any specific optimization software as in the case of AMPL, nor using a specific server or having a software license to obtain the corrosion matches. The results were obtained using the open-source software R-project, but other software alternatives are available (e.g., Python), which facilitates the implementation in practical applications. The proposed approach represents an interesting alternative for corroded pipelines based on its independence from the controlling parameters in the annealing case.

Although the nearest-neighbor criterion could be affected by possible clusters (or closed neighbors), the obtained results suggest that the matching ratio was higher than in the annealing case. The latter, considering that the proposed approach focuses on an affine transformation and incorporates potential location uncertainties in both directions. This statement should be further confirmed with other synthetic and real datasets considering clustered and dispersed point patterns.

Matched defects and descriptors of new defects

The matching results allow classifying each defect as new or old, where the former refers to those defects reported only in the second inspection and the latter in both inspections. These new defects can be attributed to a miss-detection in the first inspection, a possible false alarm in the second inspection, or a purely new defect produced by the soil aggressiveness or the fluid properties. Straub [START_REF] Straub | Generic approaches to risk based inspection planning for steel structures[END_REF] remarked that determining which is the reason for the appearance of these new defects may be impossible if the real condition of the pipeline is unknown. However, something that is certain is that there is an "indication" of new defects. This section describes the main features of both new and old sets.

The iterative approach identified both likely and farther matches in the previous example, i.e., defects with a separation higher than δ. Let us focus first on how the matching ratio, #Matches/Total points per segment, changes when the farther pairs are also contemplated. For this purpose, consider only pipe segments between consecutive joints without purely new defects (or initial mismatching), i.e., where no defects were only reported in the second inspection, and segments with possible false alarms -only records in the first inspection-. If the farther are included, the matching ratio increases from 52% to 65% for ILI1-Int and 44% to 51% for ILI2-Int, respectively. The difference is shorter for the outer wall, with an average increment from 75% to 78% for ILI1-Ext and from 65% to 68% for ILI2-Ext. Note that the matching ratio decreased for the second inspection because the number of points increased between the two inspections. However, this difference can also be favored by the presence of clustering or "corrosion colonies". Considering the clustering criterion of ASME B31G (i.e., 3t), about 36% and 46% of defects not being matched in the inner and outer wall are potential clusters. It indicates that the "corrosion colonies" could hide the location of these corrosion points and classify them as part of a corrosion cluster, and the "indication" is no longer available. The matching ratio was also compared with the different soil categories to verify if there is some relationship with the number of points (Table 4). The results indicate that the matching ratio does not change drastically at the inner wall despite the number of defects varying significantly in each soil class. For the first inspection, results ranged from around 60 to 70%, where the S3 soil had a higher rate but fewer defects. The results for the second inspection range from 40 to 55%. The outer wall presents a more significant variation, which is explained by the low number of defects in some soil classes. In both cases, the matching ratio would increase around 4 to 5% on average if potential clusters were omitted, especially soils S1 and S7.

The correspondence matrix and outliers identification contemplate all the reported defects per segment joint; however, only the matches (including the farther distances) and those defects initially not matched with a separation shorter than δ = 0.11m were considered. This filter was implemented because the remaining records would not be physically matched, and it also reduces computational time. Table 5 summarizes the matching results, ratifying 84% of the preliminary matches, discarding almost all the farther matches, and replacing them with defects initially not matched. The remaining 16% of preliminary matches were initially classified as outliers. Recall that the optimization approach in Eq. 11 produces binary outputs, which can be affected by nearly cluster segments; for instance, a defect with two close neighbors in the second inspection. In this regard, about 30% (and 38%) of the preliminary matches discarded with the correspondence matrix at the inner (outer, respectively) wall were classified as possible clusters, considering the ASME B31G grouping criterion again. Bearing in mind this limitation, the old defects would include both the preliminary and the final matches to recognize also the possibility of new defects from the nearest neighbor perspective. Percentage in parentheses are the matches obtained by the correspondence and outliers stage.

The old and new defects tend to be close together. The majority of old defects are separated less than 0.2m from the new ones, as depicted in Figure 13, considering the nearest neighbor distance. The nearest new defects have an average depth within 11 to 12%t for both pipe walls. Table 6 shows the summary of the corrosion extent from both sets. It can be noticed that larger and wider defects appeared in the new defects' but with shallow depths, associated with uniform corrosion that was not initially detected in the first inspection.

Insights about the probability of detection and false alarm

Besides classifying each defect as new or old, the matches can be used as indicators of how the probability of detection (PoD) and the probability of false alarm (PFA) behaved for the inspection tool. Non-destructive Techniques (NDT) such as ILI deal with the possibility of detecting a corrosion defect given detection threshold, usually associated with a defect depth of d * = 10%t, where t stands for the pipeline wall thickness. This detection threshold is used as a filter to distinguish between defects that can be detected or not, and it is commonly used on one side to define the probability of detection as PoD(d) = P(d > d * ) [START_REF] Rouhan | Probabilistic modeling of inspection results for offshore structures[END_REF]. The filtering process could influence the burst capacity of the corroded pipeline; however, a threshold of d * = 10%t has been indicated to be negligible for this purpose [START_REF] Bao | Influence of depth thresholds and interaction rules on the burst capacity evaluation of naturally corroded pipelines[END_REF]. The probability of false alarm, on the other side, represents the case of detecting a defect, although it does not exist and is linked with the noise of the inspection tool. Different researchers have implemented these probabilities to acknowledge the uncertain number of real defects. For instance, Heidary & Groth [START_REF] Heidary | A hybrid population-based degradation model for pipeline pitting corrosion[END_REF] proposed a populationbased pitting degradation model that uses the probability of detection and the probability of false alarm (false call) to estimate the number of existing pits in the ILI measurement and those detected on a second consecutive inspection.

Wang et al. [START_REF] Wang | Failure probability assessment and prediction of corroded pipeline under earthquake by introducing in-line inspection data[END_REF] considered the probability of detection at a corroded pipeline rehabilitation process subjected to seismic activity.

Given a particular defect, the false detection would be related to a noise signal that produces a higher response than the detection threshold. According to Rouhan & Schoefs, these may be attributed to the surrounding conditions produced by humans' environmental conditions or the measuring process [START_REF] Rouhan | Probabilistic modeling of inspection results for offshore structures[END_REF]. That was confirmed by Schoefs et al. [START_REF] Schoefs | The αδ method for modelling expert judgement and combination of non-destructive testing tools in risk-based inspection context: application to marine structures[END_REF] with an analysis of human factors accuracy and the use of protocols. These two probabilities are frequently represented with two separate probability density functions from noise and defect signals, as shown in Figure 14, which sometimes are known as a pure noise signal, and a noise+defect signal [START_REF] Rouhan | Probabilistic modeling of inspection results for offshore structures[END_REF]. This figure depicts that the probability of a false alarm is the area under the noise signal's curve. Besides, the probability of detection is determined by the area under the defect signal. In both cases, after the detection threshold [START_REF] Zeng | Probability of detection model for gas transmission pipeline inspection[END_REF]. In the absence of the inspection signal, let us concentrate on the reported defects from the ILI measurements. If an ILI measurement occurs, the report would include a certain number of metal loss measurements roughly known as corrosion defects. These defects have a depth greater than the reporting threshold of 10%t; however, some questions remain unanswered: (i) Are those detected defects exist? (ii) Is there any undetected defect? These questions can be translated into the following conditional events given that a defect may be detected or not [START_REF] Rouhan | Probabilistic modeling of inspection results for offshore structures[END_REF]:

• E1: absence of a defect, conditioned to no defect detected;

• E2: absence of a defect, conditioned to defect detected;

• E3: presence of a defect, conditioned to no defect detected; and

• E4: presence of a defect, conditioned to defect detected.

Denote the detection output by d e (X) and the existence of a defect as X. These outputs can be seen in a binary point of view, i.e., d e (X), X ∈ {0, 1}; therefore, P(E 1 )+P(E 3 ) = P(E 2 )+P(E 4 ) = 1. Following this notation, the probability of detection and false alarm can also be defined as PoD(X) = P(d e (X) = 1 | X = 1) and PFA(X) = P(d e (X) = 1 | X = 0). These definitions were used with the Bayes' conditional theorem by Rouhan & Schoefs to express the previous four conditional events as shown below [START_REF] Rouhan | Probabilistic modeling of inspection results for offshore structures[END_REF]:

P(E 1 ) = P(X = 0 | d e (X) = 0) = (1 -PFA(X))(1 -γ) (1 -PoD(X))γ + (1 -PFA(X))(1 -γ) (14)
P(E 2 ) = P(X = 0 | d e (X) = 1) = PFA(X)(1 -γ) PoD(X)γ + PFA(X)(1 -γ) (15)
P(E 3 ) = P(X = 1 | d e (X) = 0) = (1 -PoD(X))γ (1 -PoD(X))γ + (1 -PFA(X))(1 -γ) (16)
P(E 4 ) = P(X = 1 | d e (X) = 1) = PoD(X)γ PoD(X)γ + PFA(X)(1 -γ) ( 17
)
where γ is the probability of presence of a defect over the inspection section, i.e., γ = P(X = 1), so (1 -γ) = 1 -P(X = 1) = P(X = 0) is the probability of the absence of a defect.

Normally, the PoD is a monotonically increasing function that depends on the defect depth, which can be approximated "experimentally" using the ratio of the actual number of defects over the total number of reported defects in a defect depth range [START_REF] Rouhan | Probabilistic modeling of inspection results for offshore structures[END_REF]. Different authors have reported this approach to describe the probability of indication in a given range of interest [START_REF] Timashev | Diagnostics and Reliability of Pipeline Systems[END_REF]. Consider a segmentation of the depth of the defects given by d i , i = 0, . . . , m where d m ≤ max(d ILI ) being d ILI the reported depths at the ILI measurement, whereas d 0 = 0 represents an intact pipeline.

Accordingly, denote by n i for i = 1, . . . , m the number of defects correctly identified by the inspection tool, over a total of N i , for each segment [d i-1 , d i ], then the probability of true indication (PTI) can be estimated by PT I ∼ n i /N i [START_REF] Timashev | Diagnostics and Reliability of Pipeline Systems[END_REF]. The problem of this approach for PoD without a field validation is that the defect's existence is assumed but cannot be guaranteed.

On the contrary, consider the event E 4 , which evaluates the probability that a defect exists given it was already detected. This event was used assuming detection in both inspections 2 it does exist. Note that this approach would estimate the real number of existent defects from the first inspection; however, this assumption would neglect those repaired defects, and it may include the unlikely case where a non-existent defect was reported in both inspections.

Nonetheless, this is the best approximation based on the available information considering the lack of field sampling.

The obtained results are shown in Figure 15, including exponential and log-logistic fits based on the equations depicted in Eq. 18 and Eq. 19, once the defect depth is segmented using equidistant quantiles. Similar expressions have been implemented by other authors to estimate the PoD such as Park et al. [START_REF] Park | Comprehensive framework for underground pipeline management with reliability and cost factors using Monte Carlo simulation[END_REF], Baskaran et al. [START_REF] Baskaran | Probability of detection modelling in eddy current NDE of flaws integrating multiple correlated variables[END_REF], or Yazdi et al. [START_REF] Yazdi | Operational subsea pipeline assessment affected by multiple defects of microbiologically influenced corrosion[END_REF].

f exp (d) = γ d 1 -exp -d λ d , (18)
f Logl (d) = exp (α d + β d ln(d)) 1 + exp (α d + β d ln(d)) . (19)
These fittings were determined using the modified nonlinear least-squares function nlsLM of the package minpack.lm in R, which incorporates the Levenberg-Marquardt algorithm. The obtained coefficients, as well as their confidence intervals and standard error, are displayed in Table 7. Overall the exponential and log-logistic functions agree very well with the "experimental" data obtained from the depth of the matched defects from the first inspection. The probability reaches 1 near 30%t deep for the inner wall, which means that any defect deeper than 30%t already detected would guarantee its presence on the inner wall. Timashev & Bushinskaya also reported a similar percentage by establishing that the probability of detection of an ILI tool is almost one [START_REF] Timashev | Diagnostics and Reliability of Pipeline Systems[END_REF]. For the outer wall, a depth of 30%t achieves a lower probability of 0.8, partially explained by the amount of data at the outer wall and the overall more profound defects than those located at the inner wall. Following a similar procedure, the probability of event E 2 was determined, considering the defects that could not be matched, as the counting variable n i . As it was remarked before, the definition of E 2 and E 4 follows that P(E 2 ) = 1 -P(E 4 ), so additional functions were not necessary, and the "experimental" results were compared with the complement of the functions obtained from Table 7 (see Figure 16). The results show a good agreement between the fitted functions and the data of the unmatched defects, although this result could be expected using the definition of the raw probability of no matches, i.e., the complement from each depth segment. 2 Initially detected in the first inspection and latter matched with another defect from the second inspection. Nevertheless, this pattern is obtained regardless of the data used to fit P(E 4 ) and to calculate P(E 2 ). To demonstrate this pattern, consider a cross-validation approach using 70% of the data as a training set to determine the exponential and log-logistic fits for P(E 4 ), and the remaining 30% to estimate P(E 2 ) (test set). This procedure was implemented 100 times using random sampling, and the Mean Standard Error (MSE) from the 1 -P(E 4 ) predictions and the estimates of P(E 2 ) were determined. The results shown in Figure 17 indicate that these probabilities fit quite well by obtaining an average MSE of 0.002 and 0.007 for the inner and outer wall, respectively, which confirms the results obtained in Figure 16.

Although a complete PoD description is unknown, it is commonly assumed to behave as in Eq. 18 and Eq. 19 with exponential or log-logistic fits [START_REF] Straub | Generic approaches to risk based inspection planning for steel structures[END_REF][START_REF] Zhang | Cost-based optimal maintenance decisions for corroding natural gas pipelines based on stochastic degradation models[END_REF][START_REF] Rodriguez | Part ii: Development of a general failure control system for estimating the reliability of deteriorating structures[END_REF]. Inspection vendors usually provide a reference PoD for depths from 10 to 15%t of 0.90, commonly used to fit these theoretical probabilities. For instance, Zhang & Zhou implemented the exponential version (Eq. 18) with γ d = 1 and q = 1/λ d , where q represents a detection capability inherent from the inspection tool. The objective is to avoid this theoretical PoD but take advantage of the probabilities of E 4 or E 2 to estimate the PoD. Recall from Eq. 17 that P(E 4 ) depends on the probability of detection, false alarm, and defect existence (γ). From these probabilities, inspection vendors commonly report a PoD reference point, whereas PFA is known to be independent and constant for the defect size [START_REF] Rouhan | Probabilistic modeling of inspection results for offshore structures[END_REF]. Therefore, if γ is estimated, and the reference point is used, an estimate of the probability of false alarm can be solved using Eq. 17. Following Timashev & Bushinskaya, the probability of true detection and thus the existence of a defect would be at most [START_REF] Timashev | Diagnostics and Reliability of Pipeline Systems[END_REF]:

γ = A d A p
where A d is the surface area occupied by all the existing defects, and A p is the total surface area being inspected. These surface areas would not contemplate the depth of the defects, and they depend on the inner and outer diameters 3 for the inner and outer walls, respectively. The existing defects were again assumed to appear in both inspections, but this probability could be higher by including additional defects. Regarding the PoD, a reference point of PoD(15%t) = 0.90 was chosen based on the information reported by the inspection vendor for different MFL inspection tools.

Based on the mentioned above, the detection capabilities in each inspection, the probability of existence, and the probability of false alarm were determined for both pipe walls (Table 8). The results indicate that γ for the inner wall is about five times the probability for the outer wall, which is consistent with the number of defects and the obtained matches. However, the difference is shorter between the two obtained PFA. These false alarm probabilities are not far away from those reported by other approaches [START_REF] Wang | A Bayesian model framework for calibrating ultrasonic in-line inspection data and estimating actual external corrosion depth in buried pipeline utilizing a clustering technique[END_REF][START_REF] Howard | On the probability of detecting wall thinning defects with dispersive circumferential guided waves[END_REF]; however, they are specific for the inspection tool. Finally, based on the results in Table 8, the probability of detection (PoD) was estimated using the PFA, γ, and the P(E 4 ) from both pipe walls. For this purpose, each "experimental" point from Figure 15 were contemplated to solve again Eq. 17. The results are depicted in Figure 18 with green dots. This figure also includes a log-logistic fit and the expected theoretical function using the exponential fit with the detection capability. Note that the expected PoD has a clear difference with the log-logistic fit, but this difference is almost diminished after a PoD of 0.9, considering the procedure to determine both false alarm probabilities. The estimated PoD would be more conservative considering a lower detection probability that requires further tests and field validations.

This figure demonstrates that determining the probability of detection might not be straightforward; the proposed approach seeks to provide insights about PoD and PFA based on the data from consecutive inspections. Although the 3 The inner diameter is given by D i = D -2 t i , where D is the outer diameter, and t is the wall thickness at the i th location. matching process may not identify all the existent defects in the pipeline, which could affect both P(E 4 ) and P(E 2 ), the obtained results are fair enough based on the available information. Figure 17 supports this result with the MSE outcomes under a cross-validation approach. The estimated PoD and PFA could support future intervention decisions considering whether a defect exists or not based on the defects detected by the ILI tool; for instance: 1) contemplating a predictive maintenance approach that dynamically adapts to imperfect repairs and replacements [START_REF] Huynh | An adaptive predictive maintenance model for repairable deteriorating systems using inverse Gaussian degradation process[END_REF]; 2) quantifying how would be the economic gain by doing an inspection or a preventive maintenance [START_REF] Yuan | Estimation of the value of an inspection and maintenance program: A Bayesian gamma process model[END_REF]; and 3) predicting the distribution of corrosion defects' depth under adaptive time-dependent environmental conditions [START_REF] Kim | Adaptive approach for estimation of pipeline corrosion defects via Bayesian inference[END_REF]. Besides, this information can be used to support reliability calculations under a noisy data context; for instance, considering a

Gamma Process [START_REF] Hazra | Approximate Bayesian computation (ABC) method for estimating parameters of the gamma process using noisy data[END_REF], multiple correlated defects [START_REF] Yu | Subset simulation-based reliability analysis of the corroding natural gas pipeline[END_REF], or evaluating a Spatio-temporal analysis that benefits from FEM simulations [START_REF] Aryai | Time-dependent finite element reliability assessment of cast-iron water pipes subjected to spatio-temporal correlated corrosion process[END_REF]. Further inspections and field validations would provide more realistic results, aiming to deal with the uncertainties inherited from the inspection tool, which, in turn, would seek to avoid unnecessary or costly repairs. Unfortunately, the inspection vendor did not provide the PoD and the PFA, but it would be interesting to compare how far they are from the predictions obtained with the matching process. Further estimates could also help to evaluate the consequence of the detection in terms of repair decisions, considering, for instance, the probability of good and wrong assessments [START_REF] Sheils | Development of a two-stage inspection process for the assessment of deteriorating infrastructure[END_REF][START_REF] Sheils | Investigation of the effect of the quality of inspection techniques on the optimal inspection interval for structures[END_REF].

Non-destructive techniques are subjected to different types of randomness affecting the acceptance and rejection decisions of defects; PoD gives a path to evaluate the pipeline's condition, including these uncertainties. Field validations would not only allow evaluating conditional probabilities given a detection, but also the joint probabilities based on a correct detection and existence [START_REF] Timashev | Diagnostics and Reliability of Pipeline Systems[END_REF]:

• Probability of True Detection (PTD) P(d e (X) = 1, X = 1);

• Probability of True Non-Detection (PTND) P(d e (X) = 0, X = 0);

• Probability of False Non-Detection (PFND) P(d e (X) = 0, X = 1); and

• Probability of False Detection (PFD) P(d e (X) = 1, X = 0).

According to Timashev & Bushinskaya, these probabilities can be interpreted from a geometrical perspective by considering the area occupied by each case after a field validation [START_REF] Timashev | Diagnostics and Reliability of Pipeline Systems[END_REF]. Consider a pipeline segment with a length of L and diameter D, and let Figure 19 represent the output after the field validation. Consider that the black spots correspond to the defects that exist and were detected; the contoured spots with dash lines were defects not detected, but they do exist, and the light spots with continuous lines were non-existing defects reported by the inspection tool.

Finally, the remaining area corresponds with an intact pipeline that was correctly not reported any defect by the

Conclusions

Pipeline integrity evaluation requires monitoring the evolution of the corrosion attack in time. Location and detection uncertainties in ILI measurements complicate this task, so a matching approach is essential. This work proposed a new approach using a Voronoi partition to filter initial possible matches, which is later implemented under an iterative approach with closed-forms "best" affine transformations and a correspondence matrix and outlier vectors optimization approach. Although pipeline operators may have more than two ILI measurements, which are helpful for corrosion growth monitoring and modeling, the proposed approach is recommended for two consecutive inspections.

ILI measurements are implemented every 2 to 6 years, making the pipeline prone to maintenance or replacement procedures. This sequential approach could help follow the corrosion evolution more efficiently.

The main findings are summarized below:

1. The proposed approach focused on a transformation that contemplates the uncertain location of a corrosion defect in any direction by considering a nearest-neighbor perspective. This approach helps get competitive results than a traditional point matching problem with a defined affine transformation like the temperature annealing method.

2. Based on a segment of a real case study, this paper illustrated how the approach could identify 11 out of 21 possible matches. Based on the proposed matching, about 50% of the data could be matched between two inspections, obtaining a final number of matches of 14314 for the inner wall and 1714 for the outer.

3. The synthetic results indicated a true matching ratio for the proposed approach of 0.848 against 0.607 for the annealing case, after a sensibility analysis of their parameters. This matching ratio was also compared with an AMPL solver KNITRO in the NEOS server, obtaining similar (and even better) results. These results highlight the proposed approach as an interesting alternative for pipeline engineers that monitor the evolution of corrosion degradation.

4. Some insights about the probability of detection (PoD) and false alarm (PFA) were obtained based on the classification of new and old defects. For this purpose, the probabilities of the presence or absence of a defect given its detection from Rouhan & Schoefs [START_REF] Rouhan | Probabilistic modeling of inspection results for offshore structures[END_REF] were implemented from an "experimental" point of view with recognized exponential and log-logistic fitting functions. The results allow us to estimate the PoD and PFA based on the reported data, which, in turn, can be implemented for further reliability analysis of the pipeline.

Figure 1 :

 1 Figure 1: Scheme of the point matching problem.

Figure 2 :

 2 Figure 2: Scheme of the proposed matching by Dann & Dann [11].

  + β c Tr (T 3 -I 3 ) t (T 3 -I 3 ) subjected to 0 ≤ c ik , r i , s k ≤ 1, ∀i = 1, . . . , 2n, ∀k = 1, . . . , m m k=1 c ik + r i = 1, ∀i = 1, . . . , 2n; 2n i=1 c ik + s k = 1, ∀k = 1, . . . , m

Figure 3 :

 3 Figure 3: Point matching results with a) the example of Dann & Dann and b) a subset of points with the crossed nearest neighbors.

Figure 4 :

 4 Figure 4: Location uncertainty of a defect given an axial δ x and circumferential δ y uncertainties position

Figure 5 :

 5 Figure 5: Possible mismatching with the approach of Dann & Dann [11].

Figure 6 :

 6 Figure 6: Voronoi tessellations from the case shown in Figure 5.

Figure 8 :

 8 Figure 8: Scheme of the data processing.

Figure 9 :

 9 Figure 9: Voronoi tessellations and tile identification -Continuing Figure 6.

Figure 10 :

 10 Figure 10: a) Initial ILI data; b) obtained matches by the proposed approach, and c) matches considering farther nearest defects.

Figure 11 :

 11 Figure 11: Sensibility analysis of α in the correspondence linear optimization in the proposed approach

Figure 12 :

 12 Figure12: Comparison of the matching ratio of the proposed approach, the annealing method, and using an optimization solver

Figure 13 :

 13 Figure 13: Nearest neighbor distance defects in old to new set for the a) inner and b) outer walls.

Figure 14 :

 14 Figure 14: Defect and noise signals scheme

Figure 15 :Figure 16 :

 1516 Figure 15: Estimated P(E 4 ) with the exponential and log-logistic fits for the a) inner and b) outer wall matched records in ILI1.

Figure 17 :

 17 Figure 17: Mean Standard Error of P(E 2 ) and the 1 -P(E 4 ) predictions under a cross-validation approach.

Figure 18 :

 18 Figure 18: Estimated Probability of Detection with the complement log-logistic PoD fit and the theoretical exponential PoD for the a) inner and b) outer wall.

Figure 19 :

 19 Figure 19: Scheme of a corroded pipeline condition with detected and existence results. Adapted from Timashev & Bushinskaya [32].

Table 1 :

 1 Pipeline segmentation based on the USDA soil classification

	Segment*	Category	Classification

Table 2 :

 2 Summary of corrosion defects along the abscissa

	Parameter	Mean (Coefficient of Variation) ILI-1 Inner wall ILI-2 Inner wall ILI-1 Outer wall ILI-2 Outer wall
	Average depth (%t)	5.49 (0.26)	5.29 (0.27)	7.28 (0.49)	6.77 (0.46)
	Maximum depth (%t) 11.54 (0.21)	11.14 (0.19)	15.84 (0.46)	14.62 (0.43)
	Length (mm)	26.07 (0.49)	26.07 (0.43)	28.07 (0.48)	27.37 (0.44)
	Width (mm)	22.5 (0.40)	25.92 (0.53)	28.81 (0.67)	32.60 (0.75)
	Number of defects	23708	43399	2862	4264

Table 3 :

 3 Sensitibity analysis for the annealing matching approach

	Parameter	α	β c	λ	Matching ratio
		1e -05	200T	0.96	0.339
	Analysis of α	1e -04 1e -03 1e -02	200T 200T 200T	0.96 0.96 0.96	0.342 0.388 0.608
		1e -01	200T	0.96	0.614
		1e -04	300T	0.96	0.340
	Analysis of β c	1e -04 1e -04 1e -04	400T 800T 1000T	0.96 0.96 0.96	0.344 0.344 0.342
		1e -04	2000T	0.96	0.342
		1e -04	200T	0.9	0.346
	Analysis of λ	1e -04 1e -04 1e -04	200T 200T 200T	0.92 0.94 0.98	0.344 0.340 0.341
		1e -04	200T	0.99	0.341

Table 4 :

 4 Matching results based on the soil type

		Matches including farther pairs	Matches without farther pairs		Number of defects*
	Soil	Inner wall	Outer wall	Inner wall	Outer wall	Inner wall	Outer wall
		ILI1	ILI2	ILI1	ILI2	ILI1	ILI2	ILI1	ILI2	ILI1	ILI2	ILI1	ILI2
	S1	68.0% 46.2% 68.4% 63.5% 60.6% 41.1% 65.9% 61.2% 1550	2283	405	436
	S2	68.6% 46.3% NR	NR	63.8% 43.1% NR	NR	1584	2346	0	0
	S3	72.6% 48.3% 81.7% 53.8% 65.6% 43.6% 78.3% 51.6% 931	1400	60	91
	S4	64.6% 43.5% 72.3% 48.5% 56.0% 37.7% 64.6% 43.3% 2774	4120	65	97
	S5	66.3% 48.6% 75.0% 39.1% 59.3% 43.5% 75.0% 39.1% 5332	7268	12	23
	S6	64.4% 49.1% 67.9% 43.2% 57.9% 44.2% 67.9% 43.2% 2230	2922	28	44
	S7	57.6% 54.2% 59.7% 35.1% 53.8% 50.7% 56.1% 33.0% 936	994	1456 2478
	UZ	58.4% 50.5% 59.1% 53.0% 52.0% 45.0% 54.9% 49.3% 7069	8167	718	800
	Complete	63.1% 48.6% 61.7% 42.8% 56.6% 43.5% 58.1% 40.3% 23655 30726 2763 3985
	NR. No defects reported
	*Segments with possible false alarms or mismatches were not considered

Table 5 :

 5 Matching results after applying the correspondence matrix

	Dataset	Inner wall	Outer wall
		ILI1	ILI2	ILI1	ILI2
	Number of defects	23708	43399	2862	4264
	Preliminary matches	12980 (84%) 12980 (84%) 1592 (83%)	1592 (84%)
	Farther matches	1343 (0%)	1343 (0%)	95 (5%)	95 (5%)
	Correspondence matches	12200	12200	1455	1455
	Final matches	14314	14314	1714	1714

Table 6 :

 6 Summary corrosion extent of the new and old sets

	Set	Parameter			Inner wall					Outer wall
			Min Q1 Q2 Mean Q3 Max Min Q1 Q2 Mean Q3 Max
		Depth	10	10	11	11.64	12	36	10	11	14	16.8	20	70
	old	Length	10	18	22	25.31	30	85	10	19	24	27.47	33	109
		Width	14	18	20	26.82	30	255	16	19	28	35.69	41	270
		Depth	10	10	10	10.89	11	36	10	10	11	13.11	14	36
	new	Length	10	18	23	26.44	32	92	10	19	25	27.3	33	132
		Width	12	18	20	25.47	28	271	16	18	24	30.52	35	813

Table 7 :

 7 Fitted parameters for the exponential and log-logistic function applied on the P(E 4 )

			Inner wall				Outer wall
	Parameters	Coefficient	Confidence interval	SE	Coefficient	Confidence interval	SE
			2.5%	97.5%			2.5%	97.5%
	λ d	0.836	0.623	1.180	0.172	0.838	0.5214	1.552	0.282
	γ d	1.115	0.966	1.365	0.085	0.887	0.7185	1.301	0.110
	α d	1.320	1.211	1.438	0.050	0.541	0.4384	0.648	0.049
	β d	2.272	1.896	2.667	0.171	1.011	0.6545	1.385	0.167

Table 8 :

 8 Detection and existence summary results

	Pipe wall Detection Capability (q)	γ	PFA
	Inner	2.42	0.2655 0.0971
	Outer	2.42	0.0503 0.0291

Note that this type of multistage approach has also been reported in defect detection[START_REF] Sheils | Development of a two-stage inspection process for the assessment of deteriorating infrastructure[END_REF] and spatial statistics[START_REF] Schoefs | Characterization of random fields from NDT measurements: A two stages procedure[END_REF]

Acknowledgments R. Amaya-Gómez thanks the National Department of Science, Technology and Innovation of Colombia for the Ph.D. scholarship (Grant No. 727 COLCIENCIAS in 2015) and the French Ministry for Europe and Foreign Affairs and Campus France in framework of the Eiffel Excellence Program (2018). The first author also acknowledges the valuable support of John Fontecha for implementing the KNITRO solver in the NEOS server and Juan José Torres for the discussions about the optimization resolution method.


Jiazi Yi

Eddy Cizeron

Salima Hamma

Benoît Parrein

Simulation and Performance Analysis of MP-OLSR for Mobile Ad hoc Networks

Keywords: ad hoc networks, link state protocol, multipath routing

come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

Nowadays, great demands for self-organizing, fast deployable wireless Mobile Ad hoc Networks (MANETs) come along with the advances in wireless portable technologies. Compared with the conventional cellular wireless mobile networks that rely on extensive infrastructure to support mobility, the MANETs do not need base stations and wired infrastructure. This future makes it useful in battlefields, emergency searches and rescue operations where fixed base stations are undesirable or unavailable. For commercial applications such as convention centers, electronic classrooms and conferences, a rapid deployment of all-on-air networks provides users with more flexible and cheaper ways to share information.

However, because of the dynamic feature of the Ad hoc networks and the instability of the wireless medium, the routing Manuscript received September 19, 2007. This work is supported by the French program RNRT (Réseau National de Recherche en Télécommunications) under the project SEREADMO (http://www.sereadmo. org).

protocols used in the traditional wired networks seems not be able to be applied in Ad hoc networks. So there are a lot of unipath routing protocols being proposed for Ad hoc networks.

In recent years, more and more multipath routing protocols are also proposed. These protocols consist of finding multiple routes between a source and destination node. These multiple paths between source and destination node pairs can be used to compensate for the dynamic and unpredictable nature of ad hoc networks. The multipath routing could offer several benefits: load balancing, fault-tolerance, higher aggregate bandwidth, lower end-to-end delay, effectively alleviate congestion and bottlenecks [START_REF] Mueller | Multipath Routing in Mobile Ad Hoc Networks: Issues and Challenges[END_REF] and security.

Optimized Link State Routing protocol [START_REF] Jacquet | Optimized Link State Routing Protocol for Ad Hoc Networks[END_REF] [3] is a proactive unipath routing protocol specially designed for Ad hoc networks. As a table-driven protocol, it changes control messages periodically to enable the individual nodes be aware of the topology of the whole network. It inherits the stability of the link state algorithm and minimizes the flooding of the control traffic by using only certain selected nodes, called MPR (Multi Point Relays).

In this study, we discuses a new multipath routing protocol called MP-OLSR based on OLSR to provide fault-tolerance, higher aggregate bandwidth and load balancing. It exchanges control messages periodically as OLSR to get the topology information of the whole networks. Based on this topology information, our Multipath Dijkstra algorithm is used to obtain the multiple paths for the routing. With the algorithm, we can get node-disjoint routes or path-disjoint routes as necessary by adjusting distinct cost functions. In the network, the packets are forwarded from the source to the destination by employing a semi-source routing mechanism (source routing with route recovery). In addition, to meet the need for the reliable transmission, multiple description coding strategy is used in the data transmission. Thus the contribution of the paper is double. First, we propose a multipath routing algorithm based on the Dijkstra algorithm that allows different multipath approaches. Routing recovery from intermediate nodes is also included.

Simulation and Performance Analysis of MP-OLSR for Mobile Ad hoc Networks

Jiazi Yi, Eddy Cizeron, Salima Hamma, Benoît Parrein Université de Nantes, Nantes Atlantique Universités IRCCyN, CNRS UMR 6597, Polytech'Nantes, rue Christian Pauc -BP50609 44306 Nantes cedex 3 -France {jiazi.yi, eddy.cizeron, salima.hamma, benoit.parrein}@univ-nantes.fr Second, the multiple routes are exploited via an original multiple description coding based on a discrete Radon transform.

The remainder of the paper is organized as follows. Section II gives a brief review of several existing routing algorithms. The specification of the protocol is presented in Section III in detail. In Section IV we introduce the application of Multiple Description Coding in the MP-OLSR. The simulation model and performance results are demonstrated in Section V, and we conclude the paper in Section VI.

II. ROUTING ALGORITHMS FOR MANET

A. Proactive Routing and Reactive Routing

Proactive Routing and Reactive Routing are two main kinds of routing protocol for Ad hoc networks [START_REF] Hamma | Performance evaluation of reactive and proactive routing protocol in IEEE 802.11 ad hoc network[END_REF].

For Proactive Routing, also called table driven routing, each node maintains a routing table containing routes to all nodes in the network. Nodes must periodically exchange messages with routing information to keep routing tables up-to-date. The routing table is calculated before needed. So it has minimal latency but also has high control overhead. The OLSR protocol mentioned above is a typical proactive routing protocol.

For Reactive Routing, also called on-demand routing, a node only tries to find a route when necessary. However, because sometimes the route could not be get immediately, the network using reactive routing usually has longer delay. The DSR (Dynamic Source Routing) [START_REF] Johnson | DSR: The Dynamic Source Routing Protocol for Multi-Hop Wireless Ad Hoc Networks[END_REF] is a well-known reactive protocol. In DSR, the source explicitly lists this route in the packet's header, identifying each forwarding "hop" by the address of the next node to which to transmit the packet on its way to the destination host.

Our approach is to get the topology information proactively and compute the routes on-demand.

B. Multipath Routing

In the literature, multipath routing protocols are often used for backup routes. Otherwise, if the goal is the repartition of information, the implementation is generally based on pure source routing. In [START_REF] Lee | AODV-BR: Backup Routing in Ad hoc Network[END_REF], the author proposes a scheme called AODV-BR to improve existing on-demand routing protocols by creating a mesh and providing multiple alternate routes. But it just provides backup routes, which means during the transmission, there is still one route is used. In [START_REF] Lee | Split Multi-Path Routing with Maximally Disjoint Paths in Ad Hoc Networks[END_REF], the author proposes an on-demand routing scheme called Split Multi-path Routing (SMR) that establishes and uses multiple routes of maximally disjoint paths. However, as we will see in the following, a pure source routing strategy seems to be not suitable for dense network. The route maintenance by spending Route Error messages will cause longer delay than route recovery. A source routing multipath OLSR is presented in [START_REF] Zhou | A novel routing protocol for ad hoc sensor networks using multiple disjoint paths[END_REF] by using the shortest path algorithm. However, the suppression of nodes in multiple calls of Dijkstra algorithm could not work for sparse networks. Furthermore, strict node-disjoint multiple paths are not suitable for partition or fusion of group of nodes that can imply temporary a single link for connection. And one more time, our work calls into question pure source routing for dense MANET when multiple paths are used.

III. MP-OLSR SPECIFICATION

The MP-OLSR can be regarded as a hybrid multipath routing protocol. It sends out HELLO messages and TC messages periodically to be aware of the network topology, just like OLSR. However, MP-OLSR does not always keep a routing table. It only computes the routes when there are data packets need to be sent out.

The core functioning of MP-OLSR has two main parts: topology sensing and routes computation. The topology sensing is to make the nodes get to the topology information of the network, which includes link sensing, neighbor detection and topology discovery. This part gets benefit from MPRs as well as OLSR. The routes computation uses the Multipath Dijkstra Algorithm to populate the multipath based on the information get from the topology sensing. The source route (the hops from the source to the destination) will be saved in the header of the data packets. The intermediate nodes just read the packet header and forward the packet to the next hop. Furthermore, to overcome some drawbacks of the source routing, the route recovery is introduced.

A. Topology Sensing

The topology sensing is to make the nodes get to know the topology information of the network, which includes link sensing, neighbor detection and topology discovery. This part gets benefit from MPRs as well as OLSR to minimize the flooding of broadcast packets in the network by reducing duplicate retransmissions in the same region.

By sending the routing control messages proactively, the node could get to know the topology of the network: its neighbors, 2hop neighbors and other links, which are saved in neighbor set, 2-hop neighbors set and topology set respectively.

B. Routes Computation

Contrary to classical OLSR, routes are not determined by nodes each time they receive a new Topology Control message, but only if need be, in order to avoid the loud computation of several routes for every possible destination. When a given source must send packets, the route computation procedure uses the algorithm shown in Figure 1.

The general principle of this algorithm is at step i to look for the shortest path P i to the destination d. Then the edges in P i or pointing to P i have their cost increased in order to prevent the next steps to use similar path. f p is used to increase costs of arcs belonging to the previously path P i (or which opposite arcs belong to it). This encourages future paths to use different arcs but not different vertices. f e is used to increase costs of the arcs who lead to vertices of the previous path P i . We can choose different f p and f e to get link-disjoint path or node-disjoint routes as necessary.

Figure 1 The Multipath Dijkstra Algorithm

In Figure 1, Dijkstra(G,n) is the standard Dijkstra's algorithm which provides the source tree of shortest paths from vertex n in graph G; GetPath(ST,d) is the function that extracts the shortestpath to n from the source tree ST; Reverse(e) gives the opposite edge of e ; Head(e) provides the vertex edge e points to.

Another possible solution is instead of increasing the cost of the arc, we delete the related node from the node set so that we can get totally node-disjoint routes, which is the ideal case [START_REF] Zhou | A novel routing protocol for ad hoc sensor networks using multiple disjoint paths[END_REF]. However, in the real scenario, especially the cases that the nodes are sparse, the delete of some "key" nodes might prevent from finding new route. As shown in Figure 2, A is trying to send a packet to E. Initially, the costs of all the links are set to 1. It first find route A→B→C→D→E according to the Dijkstra algorithm. And to get the second route, it will increase the costs of used links, according to f p , which are A→B, B→C, C→D and D→E in this case and also the links that lead to these nodes according to f e (not shown in the figure) . Then we use the Dijkstra algorithm again to get A→F→G→D→E. If we delete node B, C and D after we find the first route, there is no way to find the second one. As shown in the figure, by using the Multipath Dijkstra Algorithm, we can get three node-disjoint routes. When the number of routes required comes to ten, we get ten optimal routes (some nodes might be shared when node-disjoint routes are unavailable) according to the cost functions. In this article, only three or four routes are used.

C. Route Recovery

In the classical OLSR, the hop-by-hop routing is used, which means when a packet reaches an intermediate node, the protocol will check the routing table of the local node and then forward the packet to the next hop.

In contrast, in MP-OLSR, we use the semi-source routing approach. It will help the source node keep good control of the packets which will be forwarded in the multipath. However, in the mean time, the pure source routing might cause two problems: Firstly, the information in the source node might be not new enough because it needs time to flood the topology control messages to the whole network. It means when computing the routes, the source node might use the links that does not exist anymore. Secondly, even when the information in the source node is updated, the topology might change during the forwarding of the packet. Both of them will cause the failure of the packets forwarding.

To solve these problems, the route recovery is used: before a medium node trying to forward the packet, the node first check if the next hop in the source route of the packet is one of its neighbors. If yes, the packet is forwarded as it should be. If no, then it's possible that the "next hop" has moved out of the transmission range of the node. Then it is necessary to recompute the route and forward the packet through the new route.

For example, as show in Figure 4, node A is trying to forward packets to D. The original multiple source routes are A→B→C→D and A→E→F→G→D. However, node G moves out the transmission range of F and makes the second route unavailable anymore. When a packet is being forwarded through the second route, according to the source route stored in the packet header, and get to node F, it will first check if it's next hop G in node F's neighbor set or not. If not, node F will recompute the route to D, and get F→H→D, then forward the packet through the new route.

Figure 4 Route Recovery

Because the Route Recovery does not need to check the network information but just the topology information saved in the local node, it will not cause much extra delay and effectively improve the delivery ratio.

IV. THE MULTIPLE DESCRIPTION CODING

By adding redundancy to information streams and splitting them up into several sub-streams, we can improve the integrity of data, especially by sending these sub-streams along different paths from the source to the destination. This kind of transformation is called Multiple Description Coding (see [9] for a detailed review).

Given a piece of information I, a multiple coding method generates N independently communicable packets (D , D 2 ,…, D N ). Each description D i is generally much smaller than the original information. In the MDC scheme, the more descriptions are received, the closer to I is the reconstructed information Î. In our case, we just assume that it exists an integer M (0<M≤N) such that every subset of descriptions containing at least M different descriptions is sufficient to rebuild entirely I. Thus, the higher is M, the lower is the redundancy. In particular, M = 1 (respectively M =N) corresponds to the case where In MP-OSLR, the Mojette transform [START_REF] Servières | The Mojette transform: Discrete angles for tomography[END_REF] is used to produce different projections of the original information, each one being sent along a specific path. This discrete form of the Radon transform only requires the addition operation and is exactly invertible. For a detailed review, please refer to [START_REF] Pierre Guédon | Internet distributed image information system[END_REF].

However, applying MDC coding to original packets may significantly increase the total number of packets that are transmitted in the network, which may results in new congestion.

A possible solution consists in setting up a sending buffer and then performing the MDC on its content (that corresponds to a group of original packets).The procedure is shown in Figure 5.

V. SIMULATION AND PERFORMANCE ANALYSIS

A. Environment and Assumption

The proposed algorithm is simulated on NS2. The channel capacity of mobile hosts was set to 11Mbps. A two-ray ground reflection model, which considers both the direct and a ground reflection path, was used as radio propagation model. We use the DCF (distributed coordination function) of IEEE 802.11 for wireless LAN as the MAC layer protocol. It has functionality to notify the network layer about link breakage.

In the simulations, there are 50 nodes move in a 1000m×1000m square region for 200 seconds simulation time. The random waypoint mobility model is used and there is no data packet transmission in the first 20 seconds to make sure the nodes are well distributed and there is sufficient time for the nodes to finish the initialization process. All the nodes have the same transmission range of 250 meters.

In each simulation, there are 30 CBR (Constant Bit Rate) sources, generating 10 packets per second with the size of 512 bytes.

Figure 5 The Multiple Description Coding

B. Performance metrics

The following metrics are used to evaluate the performance of the protocols:

Packet delivery ratio: the ratio of the data packets delivered to the destination.

Routing load: gives the number of routing packets over the number of received data packets. Each routing packet sent or forwarded by a mobile is counted.

Average end-to-end delay: The end-to-end delay is averaged over all surviving data packets form the sources to the destinations. It includes queuing delay and propagation delay.

Load balancing: We use a graph G=(V,E) to denote the network, where V is the node set and E is the link set. We define a state function of f: V→I where I is the set of positive integers. f(v) represents the number of data packets forwarded at node v. Let coefficient of variation CoV(f) = standard deviation of f / mean of f. We use CoV(f) as a metric to evaluate the load balancing. The smaller the CoV(f), the better the load balancing.

C. Simulation results

We compared the performance of OLSR and MP-OLSR in the simulations.

Except for the difference of the protocols, two different kinds of strategy to discover the link failures are considered. The first one is the proactive way: when a node does not receive HELLO messages continuously from its former neighbors, it considers that the wireless link between them is broken and will delete the corresponding node from the neighbor set. The second strategy is to use link-layer feedback. When a node is trying to forward a packet to the next hop that has been out of its transmission range, its link layer will drop the packet because of MAC_RET (REtry Timeout). In the mean time, the link layer will give a feedback to the routing layer and inform the lost of the link. We also compared both of the strategies in the simulation.

The simulations are taken in 4 protocols:

The 6 shows the delivery ratio of the simulations. As we can see from the figure, the OLSR with feedback has better delivery ratio than the original OLSR. This is because the protocol with feedback could detect a link failure as soon as the first packet is dropped and recomputed the routing table. In the contrast, the original OLSR tends to continue to send the packets through a failed link until it finds out that it is not able to receive a HELLO message from the original neighbor. Both with the link layer feedback, the SR-MPOLSR's delivery ratio is about ten percent lower than the OLSR. This is because of the drawbacks of source routing that have been mentioned in the previous section and compared with OLSR, which always sends the packets through the best route, the SR-MPOLSR sends the packets in multiple routes, some of which might be more unreliable. However, the RE-MPOLSR gives better delivery ratio thanks to the route recovery.

Figure 7 gives the simulation results of routing load. Because in all the scenarios, the number of the generated control messages is the same, so the protocol with higher deliver ratio tends to have lower routing load.

Figure 8 shows the average end-to-end delay. It includes the queue delay in every node and the propagation delay from the source to the destination. The multipath routing could reduce the queue delay because the traffic is distributed in different routes. On the other hand, it might increase the propagation delay because some of the packets are sent through the sub-optimal route. As we can see from the figure, compared with the single path protocols, although the multipath routing might result in more propagation delay, it could effectively reduce the queue delay, and tends to have lower end-to-end delay. The figure also shows that although RE-MPOLSR might spend more time in recomputing and recovering the path, it still has the lowest delay because the recovery mechanism could avoid sending data packets through a failed route, which will result in more retry time and the network congestion. The multipath protocols also provide more stable end-to-end delay. The simulation of MP-OLSR with MDC is also taken. The gain of delivery ratio is very little (one or two percent) compared with the RE-MPOLSR in the scenario of 50 nodes. This is because given the low density of the nodes, sometimes the multiple paths will share a lot of links or even there is just no route at all. Then it is hard to make benefits from the MDC. In fact, the MDC tends to be more suitable for large and dense networks.

VI. CONCLUSIONS AND FUTURE WORKS

In this paper, we discussed a multipath extension to OLSR. Our Multipath Dijkstra Algorithm is used to discover the disjoint routes and the source routing with route recovery is performed to forward the packets. The new multipath protocol could offer better end-to-end delay and load balancing. Because of the overlapping radio-coverage of neighbor nodes and that the limitation of the MAC protocol can result in strong interdependence between multi-routes, the multipath protocol performance gains achieved in ad hoc networks is not as much as in the wired Internet. Moreover, a pure source routing strategy, can appear as a sub-utilization of a good topology knowledge inherent in proactive behavior. We observe that the SR-MPOLSR has worse packet delivery ratio than the OLSR with link layer feedback. However, with the routing recovery, the MP-OLSR can achieve the best performance. To meet the requirement for a reliable transmission, the multiple routes are exploited by a multiple description coding based on Mojette Transform.

The future research includes refining the incremental functions of f p and f e to make them adaptive to the specific network and optimize the redundancy allocation for MDC in the data transmission. More largely, more precise simulations of physical layer within NS2 are engaged. Model of radio propagation [START_REF] Escarieu | Outdoor and indoor channel characterization by a 3D simulation software[END_REF] in the middle of buildings could validate the importance of the multipath routing for ad hoc networks.

Figure 2

 2 Figure 2 Multiple Dijkstra Algorithm in sparse case

Figure 3

 3 Figure 3 Multiple Dijkstra Algorithm with f p (c) = f e (c) = 2c (a) Three routes (b) Ten routes

Figure 6 Figure 7

 67 Figure 6 Delivery Ratio

Figure 8 Figure 9

 89 Figure 8 End to end delay

Figure 10 Delivery

 10 Figure 10 Delivery Ratio in the scenario of 100 nodes

  original OLSR: It is the original single path routing protocol. It discovers the link failure in a proactive way.

	The implementation of UM-OLSR [12] is used here;
	OLSR with link layer feedback [12];
	SR-MPOLSR: MP-OLSR with link layer feedback, use
	source routing only. The incremental functions are: f p (c)
	= f e (c) = 2c ;
	RE-MPOLSR: MP-OLSR with route recovery and link
	layer feedback. The incremental functions are: f p (c) =
	f e (c) = 2c.
	Figure