Enhancing the accuracy of automatic eddy detection and the capability of recognizing the multi-core structures from maps of sea level anomaly

Automated methods are important for automatically detecting mesoscale eddies in large volumes of altimeter data. While many algorithms have been proposed in the past, this paper presents a new method, called hybrid detection (HD), to enhance the eddy detection accuracy and the capability of recognizing eddy multi-core structures from maps of sea level anomaly (SLA). The HD method has integrated the criteria of the Okubo–Weiss (OW) method and the sea surface height-based (SSH-based) method, two commonly used eddy detection algorithms. Evaluation of the detection accuracy shows that the successful detection rate of HD is ∼ 96.6 % and the excessive detection rate is ∼ 14.2 %, which outperforms the OW and those methods using SLA extrema to identify eddies. The capability of recognizing multi-core structures and its significance in tracking eddy splitting or merging events have been illustrated by comparing with the detection results of different algorithms and observations in previous literature.


Introduction
Mesoscale eddies play an important role in ocean circulation as well as in heat and mass transport (McWilliams, 2008;Nencioli et al., 2010). With the development of observation technology, automated detection methods are essential for understanding the complex movements and dynamic characteristics of mesoscale eddies in large volumes of altimeter data. In Eulerian schemes (distinguished from Lagrangian detection schemes using drifting trajectory data (Dong et al., 2011)), existing eddy detection methods can be categorized into three classes: (1) those based on a physical parameter, (2) those based on flow geometry, and (3) those based on sea surface height anomaly (SSHA).
Among the first category, the Okubo-Weiss (OW) approach (Okubo, 1970;Weiss, 1991) is the most widely used. The physical parameter W is computed from the horizontal velocity field as W = Sn 2 +Ss 2 −ω 2 , where Sn and Ss are the shear and strain deformation, respectively, and ω is the vertical component of vorticity. Applications of the OW method can be found in much literature (Chelton et al., 2007;Henson and Thomas, 2008;Isern-Fontanet et al., 2003b, 2006Morrow et al., 2004;Xiu et al., 2010). Despite its popularity, three major weaknesses were reported about OW. First, it is difficult to decide on the optimal or appropriate threshold value of the parameter W that makes identification more accurate (Chelton et al., 2011). Recent studies (Williams et al., 2011;Petersen et al., 2013) have developed a new R 2 algorithm to eliminate the dependency on threshold choosing, which judges the quality of an eddy based on the similarity of certain functional fits with an idealized Gaussian vortex. Second, the derivation of this physical parameter may bring extra noise that increases falsely detected eddies (Ari Sadarjoen and Post, 2000;Chaigneau et al., 2008;Chelton et al., 2011;Nencioli et al., 2010). Third, the physical criteria may fail to locate the eddies or underestimate the dimensions (Basdevant and Philipovitch, 1994;Doglioli et al., 2007;Henson and Thomas, 2008;Isern-Fontanet et al., 2003a).
To improve eddy detection, some novel approaches based on flow geometry characteristics were developed. The winding-angle (WA) approach first proposed by Ari Sadarjoen and Post (2000) is a representative that identifies eddies by clustering closed or spiral streamlines. Chaigneau et al. (2008) adapted this approach by using sea level anomaly (SLA) local extrema as potential eddy centers and the outermost streamlines as corresponding boundaries. Nencioli et al. (2010) argued that this adaption incorporated a physical quantity (SLA) and thus should be regarded as a hybrid approach. Consequently, they further developed a pure flow geometry-based, or vector geometry-based (VG) approach, which identifies eddies only by the geometry characteristics of velocity fields and is independent from parameters derived using velocity derivatives as well as from the SLA field (Nencioli et al., 2010). However, the specification and sensitivity test of two additional parameters that are used for searching eddy centers complicate the identification procedure.
The third category directly uses SSHA or SLA for eddy identification in which a threshold is nevertheless always required to delimit eddy dimensions. Fang and Morrow (2003) used a 10 cm threshold to identify the large eddies in the South Indian Ocean. Chaigneau and Pizzaro (2005) chose 6 cm to detect eddies in the region west of South America, and Wang et al. (2003) selected 7.5 cm for eddy identification in the South China Sea (SCS). In 2010, Chelton et al. (2011) proposed a sea surface height-based (SSH-based) method for global studies and eliminated the threshold dependency on SSHA. However, this method relies on other thresholds (e.g., area and horizontal scale) that were specified in its detection criteria to determine the eddy boundaries, and may yield eddies with more than one local extremum, which are referred to as multi-core structures in this paper.
Multi-core structures exist in the vicinity of mutual interactions between eddies, so that identifying and tracking them may allow one to understand how eddies interact with each other, how they split or merge, and how the energy transfers and exchanges. However, few studies have developed algorithms capable of detecting and tracking them. Chelton et al. (2011) tried splitting the multi-core structures, only to find extra induced problems in tracking them, so they eventually abandoned the splitting procedure. This paper presents a hybrid detection method (HD) that attempts to identify and characterize the multi-core structures and lays the foundation for further tracking the evolution processes correctly. A companion tracking method that addresses the difficulties in tracking multi-core structures will be elaborated in another paper. This paper is mainly focused on the discussion of the new detection method and its performance.
The detection accuracy defines the quality of eddy identification algorithms and is a necessary evaluation criterion for newly developed methods. Integrating different methods presents a promising avenue to improving the performance of detection algorithms, for it tries to make the best of the advantages and minimize the disadvantages of each single method. In previous literature, the adapted WA method (Chaigneau et al., 2008) incorporated SLA extrema to make the algorithm more simple and efficient. Morrow et al. (2004) combined additional SLA criteria to reduce excessive detections of the Q-parameter approach. In this paper, the HD method is developed with the integration of the OW method and the SSH-based method; it combines the W parameter criterion and SLA local extrema to detect eddy centers and uses SLA contours to delineate the dimensions of eddies. An objective validation protocol that has been used by Chaigneau et al. (2008) and Nencioli et al. (2010) is adopted to measure the detection accuracy of the HD method. Meanwhile, the detection results of different algorithms are compared to examine the performance of HD and the ability of characterizing multi-core structures. Nan et al. (2011)'s study of three anticyclonic eddies among which splitting and merging events had been observed is used to verify HD's detection results and to demonstrate the significance of recognizing multi-core structures.
The rest of the paper is organized as follows: Sect. 2 introduces the altimetry data and the detailed procedures of the HD method; Sect. 3 presents the results and discussion of HD's detection evaluation, method comparisons, and historical verifications; and the final section provides the conclusions.

Eddy detection methodology
Two major procedures are involved in conventional eddy identification: the first is to recognize the center of the eddy, the second is to extract the boundary. However, to identify the multi-core structures, HD includes an additional procedure to determine if the eddies constitute a multi-core structure of which the closed boundary contains more than one local extremum of the same polarity. The methodology of the new method thus has three steps: (1) identify eddy centers; (2) find eddy boundaries; and (3) distinguish multi-core structures.

Step 1: a hybrid way to identify eddy centers
The velocity field within an ocean eddy is dominated by rotation, corresponding to the negative W values in the OW method that can be used to recognize where eddies exist. Among the studies using the OW method, many adopted the criterion W < −0.2 σ w (σ w is the spatial standard deviation of W ) for detecting ocean eddies from altimetry data, but this criterion usually brings undesirable noise into the data and causes false detections (Isern-Fontanet et al., 2006;Henson et al., 2008;Xiu et al., 2010).
Having an SLA maximum (or minimum) inside the eddydominant region is another essential character of ocean eddies, which has been adopted in the WA method, the SSHbased method as well as other eddy detection methods; hereafter they are all referred to as the SLA extremum-based methods. Following the definition that "a vortex exists when instantaneous streamlines mapped onto a plane normal to the vortex core exhibit a roughly circular or spiral pattern [. . . ]" (Robinson, 1991), Chaigneau et al. (2008) improved the WA method by incorporating the SLA extrema criterion. Besides, the threshold-free SSH-based method developed by Chelton et al. (2011) for global eddy detection also required qualified eddies to contain at least one local minimum or maximum. However, it should be noted that noise on SSHA or SLA fields may form false maxima and minima and therefore could induce spurious detection of extrema.
Based on these definitions of an eddy and the corresponding eddy detection methods developed in previous studies, we argue that an eddy should have a core area where the velocity field is dominated by rotation and the streamlines show a circular or spiral pattern, and an SLA maximum or minimum inside. So, this study developed a hybrid detection algorithm that has integrated both the W parameter criterion and SLA extrema criterion to effectively reduce excessive detections and enhance the detection accuracy. Here, to avoid confusion, the local extrema of SLA in eddies are referred to as "eddy centers" while the connected regions where W < −0.2 σ w are referred to as "eddy cores". Figure 1 shows the details of step 1. Each global SLA map is first clipped to the study region (e.g., the South China Sea). Then, the program calculates the W parameter and searches for local extrema within that clipped SLA. Regions where W < −0.2 σ w are extracted as eddy cores. Meanwhile, the local SLA extrema are searched out by running a 3 × 3 grid window (the smallest scale) on an SLA map. If the central point is a maximum or minimum of the window, then it represents a local extremum, and the extremum will be qualified as an eddy center only if it is located in a core area.
In this procedure, extra restrictions can be added based on the data quality of different study regions. In the South China Sea, for example, eddy centers located over water depths shallower than 100 m are removed due to alias from tides and internal waves contained in the SLA data over the shallow shelf area (Yuan et al., 2006). In addition, a smoothing algorithm, like a half-power filter (Chelton et al., 2011) or Hanning filter (Penven et al., 2005), can also be added to reduce the noise in the W field, but to avoid removing any physical information, this study applied no smoothing algorithms.

Step 2: definition and extraction of eddy boundaries
Existing algorithms lack a consistent definition of eddy boundaries, and detection results can vary significantly depending on the rules for extracting the boundary (Nencioli et al., 2010). The OW method defines enclosed regions as eddies where W values satisfy the criterion. However, this criterion is restricted to the core of the vortices and may underestimate eddy dimensions (Basdevant and Philipovitch, 1994). The WA method and the VG method represent eddy boundaries by streamlines, which make better approximation of eddy shapes, and given the good agreement with streamlines, the SSH-based method takes the outermost closed contours of SSHA as eddy boundaries. So, in the HD method, we also use SLA contours to represent eddy boundaries, but the "outermost" criterion is further refined to make the boundary extraction more sensible. The detailed procedures are shown in Fig. 2. First, a series of contours of an SLA map is generated with increments of 0.5 cm, although 1 cm increments had been confirmed good enough to resolve the eddies with compact interiors (Chelton et al., 2011). The smaller increment is preferred for it provides a better approximation of eddy dimensions. We carefully checked the detection results of different contour intervals and found that contours with 1 cm increments may underestimate the dimensions of eddies (Fig. 3a). We also found that 0.25 cm improves the boundary approximation, but the time consumption increases significantly when the increment becomes much smaller. So, the increment of 0.5 cm represents a compromise choice. Contours that are unclosed or have a diameter greater than 500 km are removed from the data set.
Second, the eddy centers and cores identified in step 1 are combined to decide which contours surrounding the eddy centers should be their corresponding boundaries. For each eddy center, the smallest contour that encompasses its core area is regarded as the qualified boundary. This definition improves the approximation of eddy shape while choosing the outermost closed contours, as boundaries tend to enlarge the dimensions. Figure 4a illustrates the difference between the two criteria.
It is noteworthy that some special situations may occur during the process of extracting eddy boundaries. For instance, it is possible that eddies may have no contours completely containing their core areas (Fig. 4b), or even no closed contours around the eddy center (Fig. 4c). In the former case, eddy boundaries are defined by the outermost closed contour that intersects the core areas, while in the latter case, the eddy dimensions are defined by the core, just as the OW method does. While some may argue about the existence of these weak eddies in the latter case, this study retains them during the identification process, as they may be in the immature or unstable stages (e.g., forming or disappearing) of their evolution processes and removing them will break a coherent process into discrete evolution pieces when they are tracked. Figure 4d illustrates another case that eddy boundaries contain extra local extrema but not qualified eddy centers. In such cases, the boundaries improperly enlarge the dimensions, and so are removed before the boundary extraction procedure. Moreover, contours containing two or more eddy centers of opposite polarities are also discarded in this study.
We summarize all the criteria or constraints used for detecting eddy centers and boundaries in the HD method as follows: 3. only those SLA maxima or minima that lie within core areas are qualified as eddy centers; 4. eddy centers located over shallow depths (taken to be 100 m in the SCS) are discarded; 5. the contours of SLA are generated with increments of 0.5 cm; 6. unclosed contours or diameter greater than 500 km are discarded; 7. the minimum closed SLA contours that encompass eddy core areas are defined as qualified boundaries; 8. if no closed contours completely contain the core areas, then take the outermost closed contours that intersect with the core areas as eddy boundaries; 9. if no closed contours intersect with the core areas, then use core areas to define eddy boundaries.

Step 3: multi-core recognition and boundary restoration
The multi-core structures of eddies, which have two or more closed eddies of the same polarity within the boundaries, represent the important transitional stages of their lives in which component eddies may experience splitting, merging or other energy-transferring interactions. They were little mentioned and had no clear definition in previous studies. The SSHbased criteria can yield multi-core structures, but they fail to characterize their shapes and track their evolution processes. The HD method developed in this study attempts to recognize eddy multi-core structures fully so that we can address the difficulties they bring to tracking algorithms. After step 2, every eddy center has been associated with a unique boundary. This step further defines that if the boundary of an eddy contains other eddy centers, it reveals a multicore structure, and the boundary is taken as the "composite" border of the whole structure (Fig. 5a). If there are multiple composite borders, the algorithm just retains the outermost one. The borders of inner component eddies, which we called "footprint" borders, are defined by the outermost closed contours that only contain one eddy center. In brief, HD characterizes the spatial domain of multi-eddy structures and included eddies by the composite borders and the footprint borders, respectively.
After recognizing multi-core structures, the procedure needs to check every eddy's boundary and extract the composite borders and footprint borders. For example, in Fig. 5a, after step 2, sub-eddy C's boundary was found to be the outermost boundary containing sub-eddies A and B, and therefore it should be recognized as the composite border. Subeddy B's boundary conflicts with the definition of a footprint border for it contains eddy A, so that subsequent restorations are needed to assign the correct boundary. Figure 5b shows the correct borders of the multi-core structure and its component eddies after the boundary restoration based on the definitions of composite borders and footprint borders.

Detection accuracy in the South China Sea
The accuracy is an important measure of the eddy identification quality. According to the validation protocol in Chaigneau et al. (2008), two quantities can be used to evaluate the accuracy of identification algorithms: the success of detection rate (SDR) and the excess of detection rate (EDR). Definitions are shown as follows: Where N c corresponds to the number of the eddies identified by both experts and the automated algorithm, N e corresponds to the number of the eddies only identified by the experts, and N om corresponds to the number of the eddies only identified by the algorithm. This study randomly selected ten SLA maps of SCS (100 • E-125 • E; 5 • N-26 • N) from 1992 to 2012 to compute SDR and EDR of the HD method. Five experts were invited to manually detect eddies on these sample maps. On each map, only eddies recognized by at least three experts were counted in the detection result used for accuracy evaluation. As subjective bias usually exists among the experts, estimating the uncertainty of each derived manual detection result is necessary to understand how difficult it may be to correctly identify eddies from a specific SLA map, and to what extent the detection error of automatic algorithms can be tolerated. Here, we introduce the ratio of inconsistency (RI) to measure the uncertainty of expert results. The definition is as follows: Where RI i denotes the inconsistency rate of expert i. N i is the number of eddies identified by expert i and N ce is the number of eddies that are commonly recognized by at least three experts. N i may be less than N ce , and so to avoid negative values we only compute the absolute value of RI. Figure 6a shows the comparison between expert detection (circles) and the HD result (dots) on 9 May 2007, SCS. There are 25 eddies (N e = 25) detected by the experts on this data set (12 cyclonic and 13 anticyclonic). Four eddies were overly detected (N om = 4) and one was missed (N c = 24) by the HD. So, according to the definition, SDR = 96.0 % and EDR = 16.0 % in this experiment. Table 1 presents the evaluation results of the ten randomly selected SLA maps. About 22 eddies were identified in each SLA map by the experts on average. The mean inconsistency rate is ∼ 20.7 %, indicating that such a number of eddies raises disagreement among experts, and when evaluating the detection accuracy, the rate value can be used as  an upper bound of the acceptable tolerance for detection errors. The detection accuracy of the OW method and the SLA extremum-based method were also presented in Table 1 for comparison.
The results show that all the detection methods have relatively high SDRs. The methods using SLA local extrema to identify eddies maintain 100 % SDR on every map, which confirms that the SLA local extrema are necessary evidence for the existence of eddies. The SDRs of the OW method and the HD method are both equal to ∼ 96.6 % on average; ∼ 3.4 % eddies failed to be recognized by either OW or HD. Further examination reveals that these overlooked eddies were not located in the intense vorticity regions (where W < − 0.2 σ w ) but nearby (Fig. 6b), and were hence excluded from the qualified eddy centers identified by the HD method.
Second, EDR of different methods in the results varies wildly. The EDR of OW is high at ∼ 70.3 % on average and even exceeds 100 % in single experiments (e.g., on 19 May 2007), which is evidence of the over-detection weakness of the OW method. The average EDR of the SLA extremum-based method (∼ 36.8 %) is much lower than that of OW, which agrees with the validation results in Chaigneau et al. (2008) and confirms that the WA method is more accurate than OW. The HD method has the lowest EDR in single experiments and on average. This result substantiates that the integration of the OW and SSH-based criteria effectively improves the detection accuracy. In addition, the average EDR of HD (∼ 14.2 %) is under the tolerance bound (20.7 %) of acceptable detection errors, which further confirms the feasibility as an automatic algorithm in practical applications. In-depth investigation shows that the ∼ 14.2 % excessively detected eddies are mainly ascribed to two reasons: most of them (1) constitute a part of multi-core structures (Fig. 6c); or (2) are located near the boundary of the study region (Fig. 6d). Both cases tend to be overlooked in manual detection.
Another analysis using the receiver operating characteristic (ROC) graphs (Fawcett, 2006) was also performed to compare the detection accuracy of different methods, and the analysis result (Fig. 7) reconfirms the good performance of the HD method over the OW method and the SLA extremumbased method. Most of the ROC points of the HD method are located near the top-left corner, which represents perfect detection accuracy. The OW method's performance is the poorest as most of the ROC points are clustered in the top-right corner.

Detection results in the Eastern South Pacific Ocean
Although Sect. 3.1 has provided a quantitative evaluation of the detection accuracy of different methods in the SCS, this section intends to further examine the detection results on a specific SLA map in a different region and to observe the detection capability of different algorithms. We employed different detection methods to identify eddies in the East- detection result as well as the automatic methods' results on that map.
By comparing with the expert result, the detection performance of each method within our test can be sorted from good to poor: SSH-based method, HD, WA, OW, and VG. The OW method has a 100 % SDR, but the EDR is relatively high (71.9 %). The detection accuracy of VG (SDR = 50.0 %, EDR = 6.3 %) is not as good as it was in Nencioli et al. (2010) in which high-resolution data sets were used. The limited resolution of the SLA maps could be the main reason. The HD, WA and SSH-based methods all have a good SDR (90.6 %), but the SSH-based method outperforms the other two for its EDR (6.3 %) is minimal. However, considering the capability of detecting eddies in multi-core structures, only HD has developed algorithms to characterize the boundaries of component eddies and composite wholes. WA was unable to recognize such mixing structures, and the SSH-based method, although it had recognized two composite structures located around 17 • S, 78 • W, failed to distinguish the contained eddies.
In brief, this analysis has demonstrated that HD, WA, and the SSH-based method are more suitable for detecting ocean eddies from SLA maps than OW and VG. Moreover, HD has the advantage of characterizing eddies in the multi-core structures over other methods, which provides the opportunity to explore eddy interactions.

Detection verification with previous literature
Eddies that have been studied in previous literature provide good references for verifying the detection results of automatic algorithms. Among the extensive research on mesoscale eddies in the South China Sea (Wang et al., 2008; Chelton et al. (2011). The VG result is derived by the MATLAB code provided by F. Nencioli. The results in these subgraphs are superposed on SLA field of 8 December 2004. Wang et al., 2006;Li et al., 1998;Jia and Liu, 2004;Yi et al., 2013), Nan et al. (2011) investigated three long-lived anticyclonic eddies (AE) in the northern SCS in 2007, two of which (named ACE2 and ACE3) experienced splitting and merging events during their evolution processes. So, this paper cites the two AEs to verify the detection results of the HD method and discusses the roles that multi-core structures play in eddies' evolution changes. Figure 9 shows four consecutive snapshots of ACE3 during the merging event. According to Nan et al. (2011), one small eddy lasted for about 2 weeks and then merged into ACE3 on 22 August. In Fig. 9, a two-core structure was found near Luzon Island on 8 August and the smaller component eddy close to the shore shrunk sharply but still remained in the composite structure after one week. On 22 August, the shrinking eddy disappeared and the two-core structure eventually turned into a single eddy. Such an evolution event of ACE3 recognized by HD is basically consistent with the description in Nan et al. (2011). Additionally, the multi-core structure is delineated, providing detailed information (e.g., how the energy tranfers) about the merging effect.
The other long-lived anticyclonic eddy, i.e., ACE2, which split off a small eddy on 11 July according to Nan et al. (2011), is also examined, and Fig. 10 shows the identification results of HD during the split. On 4 July, ACE2 constituted a two-core structure in a loose connection with ACE3, but one week later, they separated and ACE2 shrunk significantly. Afterwards, ACE2 was strengthened judging by the core areas, and contained in a more complex four-core structure on 18 July. In addition, between ACE2 and ACE3 a small newborn eddy can be observed expanding rapidly in the multi-core structure, which suggests that the evolution of ACE2 and ACE3 had contributed much to the generation and the development of this internal newborn eddy.
Some differences can be observed, relative to Nan et al. (2011), especially that ACE2 was not found splitting off a small eddy but rather separated from ACE3 during 4-11 July. Careful examination reveals how this discrepancy results from the different detection criteria used. Nan et al. employed the OW method to identify eddies so that they observed the splitting event since ACE2's core was indeed divided into two parts during that period. However, one part after the separation contained no SLA local maximum and was consequently excluded from HD's detection results. Despite some differences between the results of our algorithm and the description of Nan et al. (2011), our results demonstrate the feasibility of the automated detection even when eddy structures contain multiple cores, presenting a means to algorithmically identify splitting and merging events, enabling one to locate an entire class of events involving multiple eddy cores and a change in the number of eddies.

Conclusions
This study proposed a new hybrid method, HD, which is designed to automatically identify mesoscale eddies from SLA data sets. The combination of detection criteria of the OW method and the SSH-based method effectively enhances the accuracy of eddy detection, and procedures are implemented in the HD method to recognize the multi-core structure of mesoscale eddies, which serves as the foundation of tracking eddies' dynamic events (e.g., splitting or merging) and eddies' mutual interactions during the evolution processes.
To evaluate the detection accuracy of the HD method objectively, this study adopted the validation protocol used by Chaigneau et al. (2008) and made a comparison with the manual results detected by 5 experts. Ten random experiments in the South China Sea showed that the average successful detection rate was ∼ 96.6 % and the excessive detection rate was ∼ 14.2 %, which confirmed HD's improvement in detection accuracy over either OW or the SLA extremum-based method individually. Second, the comparisons between different detection algorithms in the ESP study area illustrated every method's performance as well as HD's capability of recognizing eddy multi-core structures. Finally, by comparing with two long-lived anticyclonic eddies investigated by Nan et al. (2011), this study verified HD's detection results and found that the merging event involving ACE3 was well recognized by HD, while the discrepancy of ACE2's splitting event was mainly due to the different detection criteria between OW and HD. This detection verification also demonstrated the additional capability of recognizing the multi-core structures, which facilitated the identification of the splitting and merging changes and unveiled details of eddy evolution that cannot be seen with other detection algorithms.
Limitations of the HD method are also obvious. The HD method can only detect ocean eddies that produce a sufficiently strong SLA. Although HD is able to address complex multi-eddy structures, the extent of mass or energy exchanges and interactions under the surface require further careful studies and validation with in situ data. So, having presented a method for the identification of multi-core structures, we expect and welcome subsequent refinement to our method.
This study has developed a hybrid method for identifying single eddies and composite structures containing multiple cores from an instantaneous map of SLA. In the next paper, we will present the companion tracking method that is able to extract the continuous evolution processes as well as complex dynamic changes and interactions of the identified eddies and composite structures.