Remote Sensing Image Subpixel Mapping Based on ... - IEEE Xplore

1306

IEEE TRANSACTIONS ON SYSTEMS, MAN, AND CYBERNETICS—PART B: CYBERNETICS, VOL. 42, NO. 5, OCTOBER 2012

Remote Sensing Image Subpixel Mapping Based on Adaptive Differential Evolution Yanfei Zhong, Member, IEEE, and Liangpei Zhang, Senior Member, IEEE

Abstract—In this paper, a novel subpixel mapping algorithm based on an adaptive differential evolution (DE) algorithm, namely, adaptive-DE subpixel mapping (ADESM), is developed to perform the subpixel mapping task for remote sensing images. Subpixel mapping may provide a fine-resolution map of class labels from coarser spectral unmixing fraction images, with the assumption of spatial dependence. In ADESM, to utilize DE, the subpixel mapping problem is transformed into an optimization problem by maximizing the spatial dependence index. The traditional DE algorithm is an efficient and powerful population-based stochastic global optimizer in continuous optimization problems, but it cannot be applied to the subpixel mapping problem in a discrete search space. In addition, it is not an easy task to properly set control parameters in DE. To avoid these problems, this paper utilizes an adaptive strategy without user-defined parameters, and a reversible-conversion strategy between continuous space and discrete space, to improve the classical DE algorithm. During the process of evolution, they are further improved by enhanced evolution operators, e.g., mutation, crossover, repair, exchange, insertion, and an effective local search to generate new candidate solutions. Experimental results using different types of remote images show that the ADESM algorithm consistently outperforms the previous subpixel mapping algorithms in all the experiments. Based on sensitivity analysis, ADESM, with its self-adaptive control parameter setting, is better than, or at least comparable to, the standard DE algorithm, when considering the accuracy of subpixel mapping, and hence provides an effective new approach to subpixel mapping for remote sensing imagery. Index Terms—Differential evolution (DE), remote sensing, subpixel mapping, superresolution mapping.

I. I NTRODUCTION

R

EMOTE SENSING sensors with a spatial resolution larger than the extent of the classes on the ground yield mixed pixels, i.e., pixels whose spectral signature is a composite of signatures of different classes [1]. These mixed pixels pose a difficult problem for land-cover mapping as their spectral

Manuscript received May 20, 2011; revised August 30, 2011 and December 1, 2011; accepted February 9, 2012. Date of publication April 11, 2012; date of current version September 12, 2012. This work was supported in part by the National Basic Research Program of China (973 Program) under Grant 2009CB723905, by the National Natural Science Foundation of China under Grants 40901213 and 40930532, by the Foundation for the Author of National Excellent Doctoral Dissertation of P. R. China under Grant 201052, by the Program for New Century Excellent Talents in University under Grant NECT-10-0624, and by the Fundamental Research Funds for the Central Universities under Grant 3103006. This paper was recommended by Associate Editor H. Ishibuchi. The authors are with the State Key Laboratory of Information Engineering in Surveying, Mapping and Remote Sensing, Wuhan University, Wuhan 430079, China (e-mail: [email protected]; [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TSMCB.2012.2189561

characteristics are not representative of any single pure landcover class [2]. Spectral unmixing techniques [3] and fuzzy classifiers [4], as effective unmixing algorithms, can obtain the fraction images with different pure classes or “endmembers.” The fraction images indicate the percentage of each pixel that is composed of each class, but they do not provide any indication about the subpixel spatial distribution of the classes within the coarse pixel; for example, they cannot provide the exact subpixel location of the interested target within the pixel. Subpixel mapping techniques [5], [6], also termed superresolution mapping [1] or downscaling [7], can be used to obtain the subpixel spatial distribution of any class, based on the fraction images, by decomposing the pixel into smaller subpixels, based on the spatial dependence phenomenon in which observations close together are more alike than those further apart [8]. One of the earliest approaches to subpixel mapping is that of Verhoeye and De Wulf [5], who proposed a deterministic solution based on linear programming. Their algorithm maps land cover within pixels from the proportions output by assuming the maximum spatial dependence within and between pixels. This simple algorithm works efficiently on simulated satellite sensor imagery and Système Probatoire d’Observation de la Terre high-resolution visible imagery, recreating the initial scenes more accurately than traditional hard classification did. Mertens et al. [9] proposed a subpixel mapping algorithm based on a genetic algorithm, using the same objective functions, as Verhoeye and De Wulf [5] did. A simple pixel-swapping strategy, as used in spatial simulated annealing, was adopted by Atkinson [10] to construct subpixel or superresolution maps. This work applied the swapping algorithm to an initial purely random subpixel map comprising the correct class fractions within each coarse pixel. Thornton et al. [8], [11] proposed a linearized pixel-swapping method for mapping rural linear land-cover features from fine–spatial-resolution remotely sensed imagery. Artificial neural networks (ANNs), as powerful tools for nonlinear prediction, have also been applied to subpixel mapping. Tatem et al. [12] trained a Hopfield neural network to optimize an initial subpixel map used for further iterations, with the simultaneous objectives of coarse fraction reproduction and spatial autocorrelation maximization; this method was then successfully tested on an actual case study [13]. Their neural network approach was further extended to account for indicator variogram models [14]. A back propagation (BP) neural network has also been used to improve subpixel mapping accuracy [15], [16]. Mertens et al. [17] combined wavelets with neural network models to account for the resolution difference between fine class labels and coarse fractions. Geostatistics provides another solution for subpixel mapping.

1083-4419/$31.00 © 2012 IEEE

ZHONG AND ZHANG: REMOTE SENSING IMAGE SUBPIXEL MAPPING BASED ON ADAPTIVE DE

Fig. 1.

1307

ADESM approach for remote sensing imagery.

Boucher et al. [1], [18] built up a sequential simulation framework based on a prior model of spatial structure or texture for generating alternative superresolution mapping of class labels. Atkinson et al. [7] proposed a superresolution mapping algorithm based on downscaling cokriging and demonstrated the performance of the method using Landsat Enhanced Thematic Mapper (TM) Plus (ETM+) images. In addition, Mertens et al. [19] proposed a direct-neighboring subpixel mapping (DNSM) algorithm, which exploits spatial dependence in a simple manner, giving accurate results in a limited computation time. The direct subpixel mapping algorithm was enhanced by using the spatial-attraction model [spatial-attraction subpixel mapping (SASM)] [20] and showed increased accuracy when compared with hardened soft classifications. More recently, a Markov random field model had been used to represent the spatial dependence within and between pixels, to give a subpixel mapping [21]. Ge et al. [22] proposed an advanced subpixel mapping algorithm to provide detailed information on the spatial distribution of land cover within a mixed pixel, but this algorithm cannot be applied to pixels falling into the boundary region of the degraded image [22], thereby losing the image information. Although the current subpixel mapping algorithms have obtained relatively satisfactory results, they can be further improved by hybrid methods or by exploring new strategies. For instance, to improve the subpixel mapping accuracy, a hybrid method combining a spatial-attraction model and a pixel-swapping algorithm [23] has been used to solve the subpixel mapping problem. In this paper, a new subpixel mapping algorithm based on the differential evolution (DE) algorithm, namely, the adaptiveDE subpixel mapping (ADESM) algorithm, is proposed. In ADESM, the subpixel mapping problem is considered as an optimization problem, maximizing the spatial autocorrelation within the pixel and the image. DE, as an effective optimization method, is used to solve the subpixel mapping optimization

problem in order to obtain the optimal subpixel mapping result. The DE algorithm, as proposed by Storn and Price [24], is a simple yet powerful population-based stochastic search and optimization technique [25]. DE uses simple real mutation and crossover operators to generate new candidate solutions in the continuous search space and applies a one-to-one competition scheme to greedily decide whether the new candidate or its parent will survive in the next generation, until finding the optimal result [26]. Due to its simplicity, ease of implementation, fast convergence, and robustness, the DE algorithm has gained much attention and a wide range of successful applications, such as numerical optimization [25], [27]–[29], mechanical engineering [30], feedforward neural network training [31], digital filter design [32], image processing [33], [34], and pattern recognition [35], [36]. Although DE has been effective in global optimization problems over continuous space, few applications of DE have been reported in the field of subpixel mapping. This is because the subpixel mapping problem needs to be solved in the discrete space as a permutation-based combinatorial optimization problem, in which each subpixel within the pixel is allocated to different classes. The effectiveness of DE for combinatorial discrete optimization problems is still considered limited [37] as the major obstacle in applying DE to combinatorial problems is caused by its working mechanism, which is based on real vectors. In addition, there are three crucial control parameters involved in DE: the population size N P , scaling factor F , and crossover rate CR. These parameters are often kept fixed throughout the optimization process and may significantly influence the optimization performance of DE [38], [39]. If we want to obtain the appropriate values of the parameters for a specific optimization problem, DE is run multiple times with different settings of the parameters [40]. This process requires high computational costs, and the time for finding these parameters is often unacceptably long for a remote sensing image process. Therefore, an adaptive

1308

IEEE TRANSACTIONS ON SYSTEMS, MAN, AND CYBERNETICS—PART B: CYBERNETICS, VOL. 42, NO. 5, OCTOBER 2012

mechanism should be considered in the optimization process. To remedy these drawbacks and utilize the global optimal ability of DE, this paper proposes an adaptive-DE approach to solve the subpixel mapping problem for coarse remote sensing images. The flowchart of the proposed algorithm is shown in Fig. 1. Accordingly, several issues should be addressed. 1) Discrete-encoding and reversible-conversion strategies. It is necessary to encode the candidate solution for the subpixel mapping in the population using a discrete value because it consists of integer numbers, e.g., the class attribute of the subpixel. However, traditional DE cannot be used in the discrete space. To solve this problem, a reversible-conversion strategy is utilized to change the candidate subpixel mapping solution from integer to real numbers, as a real transformation, and then back to integer numbers after crossover as an integer transformation. 2) Improved DE operators. After the conversion from integer numbers to real numbers, these discrete-encoding individuals may be evolved to obtain a better solution by standard mutation and crossover operators. However, the conversion process will also bring a high frequency of infeasible solutions, such as out of range. Three enhanced DE operators (repair, exchange, and insertion) were developed to repair the infeasible solutions and improve the subpixel mapping results. In addition, there is often a possibility of stagnation in DE. To move away from the point of stagnation, a local search [41], [42] operator is applied to find a more feasible solution in the local neighborhood after the selection process. 3) Self-adaptive parameter control. A self-adaptation strategy for the parameters F and CR is designed to improve the intelligence and practicability of the proposed subpixel mapping algorithm for different remote sensing images. Here, the parameters to be adapted are encoded into the individuals and undergo the evolution process using genetic operators, such as crossover, mutation, and selection. The better individuals, with better values of these encoded parameters, will be more likely to survive and produce offspring. This method reduces the time for finding the appropriate parameters, the parameter control technique, and can produce a flexible DE for remote sensing subpixel mapping. This paper proposes a novel adaptive-DE-based subpixel mapping approach for remote sensing imagery and addresses some important issues for this application, including discrete encoding, enhanced operators, and self-adaptive parameter control. We provide comparisons with several subpixel mapping methods for an artificial image and three remote sensing images. The rest of this paper is organized as follows. Section II provides the mathematical formulation of the subpixel mapping algorithm. Section III briefly introduces the basics of DE. In Section IV, we describe the proposed adaptive-DE approach for solving the subpixel mapping problem. Experimental results and analysis are given in Section V. Section VI discusses the main properties of the ADESM, in theoretical and empirical terms. Finally, the conclusions are provided in Section VII.

Fig. 2.

Example of subpixel mapping (3 × 3 subpixels).

II. S UBPIXEL M APPING P ROBLEM A. Basic Principles Spectral unmixing or fuzzy classification techniques may be used to obtain the fraction images or membership images of each class [3], [4], but they do not give the subpixel attribution of each class in the coarse pixel. A subpixel mapping method may transform the fraction images to obtain the most suitable subpixel locations for the different class fractions within a pixel. During this transformation, every pixel is divided into a number of subpixels by defining a scale factor s. The scale factor describes the difference of the resolution between the fraction images and the subpixel mapping results. That is, each pixel will consist of s2 subpixels. The basic principle of subpixel mapping is shown in Fig. 2 by a simple example with two classes (class 1 and class 2) in a raster grid of 3 × 3 coarse–spatial-resolution pixels. The fraction image of land-cover class 1 is shown in Fig. 2(a), and the proportion of another class, i.e., class 2, may be obtained by subtracting one from the fraction of the corresponding class 1 in the same pixel. If a coarse-resolution pixel zm,n is divided into nine subpixels, the scale s will be set to three, where m and n represent the mth lines and the nth samples in the image, respectively. Each subpixel corresponds to 11.11% coverage of the coarse-resolution pixel. As shown in Fig. 2(a), the fraction value of class 1 for zm,n , with 66.67%, represents six subpixels belonging to class 1 (black circle). Fig. 2(b) and (c) shows the two possible subpixel mapping solutions for Fig. 2(a). The spatial dependence principle is utilized to compare the two possible solutions, with the tendency for spatially proximate observations of a given property to be more alike than is the case for more distant observations. Land cover is spatially dependent, both within and between pixels, on the condition that the intrinsic scale of variation is not less than the sampling scale imposed by the image pixels [6]. The higher the spatial dependence, the better the subpixel mapping solutions. In Fig. 2, the solution in Fig. 2(c) is better than that in Fig. 2(b). B. Problem Formulation According to spatial dependence principles, subpixel mapping can be formulated as a maximum optimization problem, with a given appointed scale s. Suppose that the spectral unmixing model, or a soft classifier, yields fraction images for C land-cover classes and the coarse-resolution pixels are to be divided into D subpixels, where D is equal to s2 . The number of subpixels for the land-cover class i is calculated by (1). The

ZHONG AND ZHANG: REMOTE SENSING IMAGE SUBPIXEL MAPPING BASED ON ADAPTIVE DE

1309

number of subpixels that have to be assigned to land-cover class i is N Ci and has been derived from the fraction images N Ci = round(D ∗ F ractioni )

(1)

where D is the number of subpixels within a pixel, F ractioni represents the fraction value of class i for the pixel, which describes the percentage of the class, and round() is the operator that rounds its argument toward the closer integer. A measure for the spatial dependence index (SDI) SDIij will be calculated for land-cover class i and each subpixel j, considering the spatial correlation with neighboring pixels. Assuming that Nn neighboring pixels in a coarse-resolution fraction image are considered, SDIij can be expressed as a weighted linear combination of Nn neighboring fractions of land-cover class i using SDIij =

Nn 

wk · F ractioni,k

(2)

k=1

where wk is the weight and F ractioni,k represents the fraction value of class i for the kth neighboring pixel around the subpixel j. To satisfy the spatial dependence principle, wk is often calculated as the inverse of the distance of the subpixel to the coarse pixel center [6]. The subpixel mapping problem now becomes one of assigning land-cover classes to the subpixels while maximizing the SDI. The mathematical model of the subpixel mapping problem can be denoted as follows: Maximize f =

D C  

yij · SDIij

(3)

i=1 j=1

where yij is an attribute value described whether the subpixel j is assigned to land-cover class i as follows:  1, if subpixel j is assigned to land-cover class i (4) yij = 0, otherwise C 

yij = 1,

j = 1, 2, . . . , D

(5)

i=1 D 

yij = N Ci ,

i = 1, 2, . . . , C

(6)

j=1

where D is the number of subpixels within a pixel, C is the number of land-cover classes, and N Ci is the number of subpixels for the land-cover class i. The example with the scale setting of three, shown in Fig. 2, is used to describe the subpixel mapping problem. Each pixel is divided into D subpixels, where D is equal to 9 (3 × 3). For a pixel Zm,n = {zt |t = 1, 2, . . . , 9}, its fractions are equal to 66.67% and 33.33% for classes 1 and 2, respectively. The number of subpixels for classes 1 and 2 may be obtained according to the following equations: N C1 = 9 × 66.67% = 6 and N C2 = 9 × 33.33% = 3. It is worth noting that the sum of N Ci values should be equal to the number of subpixels D. Following the aforementioned description, the subpixel mapping problem can be described as the problem of how to locate this

Fig. 3. Permutation-based subpixel mapping problem.

land-cover class in D subpixels. Fig. 3 shows the representation of a possible subpixel mapping solution in Fig. 2(b), where the first, second, third, fourth, eighth, and ninth subpixels are allocated to class 1 and the other subpixels are allocated to class 2. For Fig. 2(c), the first, second, fourth, fifth, seventh, and eighth subpixels are allocated to class 1, and other subpixels are allocated to class 2, as shown in Fig. 3. According to the aforementioned analysis, the subpixel mapping problem can be formulated as a combinational optimization problem which maximizes the spatial dependence model f (y, SDI). In this paper, we utilize DE as a powerful optimal algorithm to solve the optimization problem of remote sensing subpixel mapping.

III. DE A LGORITHM The DE algorithm, as invented by Storn and Price [24], is one of the latest population-based stochastic global optimizers and is a simple yet powerful heuristic method for solving nonlinear, nondifferentiable, and multimodal optimization problems. It adopts a floating-point encoding scheme and takes advantage of the differentiation information among the population to find the global optimum in the continuous search space. The technique combines simple arithmetic operators with the classical events of crossover, mutation, and selection, to evolve from a randomly generated initial population to the final individual solution [26]. Successive populations are generated by adding the weighted difference of two randomly selected vectors to a third randomly selected vector. There are several variants of DE based on different mutation and crossover strategies [24], [43], [44]. In this paper, we use classical DE (DE/rand/1/bin) [24], [40] because this strategy is the most often used in practice. It can be described as follows. To solve the minimal optimization problem min f (x1 , x2 , . . . , xD ), DE starts with an initial population vector, which is randomly generated when no preliminary knowledge about the solution space is available. Let (k = 1, 2, . . . , N P for each Xk,G = {x1k,G , . . . , xD k,G } generation G) represent a candidate solution vector with D dimensions, where N P represents the population size. Step 1: Mutation. For each vector Xk,G in generation G, its corresponding mutant vector Vk,G = 1 2 D , vk,G , . . . , vk,G } may be calculated by the {vk,G following equation: Vk,G = Xr1 ,G + F · (Xr2 ,G − Xr3 ,G )

(7)

1310

IEEE TRANSACTIONS ON SYSTEMS, MAN, AND CYBERNETICS—PART B: CYBERNETICS, VOL. 42, NO. 5, OCTOBER 2012

where k = 1, 2, . . . , N P ; Xr1 ,G , Xr2 ,G , and Xr3 ,G represent the randomly selected vectors; and the indices r1 , r2 , and r3 are mutually exclusive integers randomly generated within the range [1, N P ], which are also different, r1 = r2 = r3 = k, so that N P ≥ 4 is required. The scaling factor F ∈ [0, 2] is a positive control parameter for scaling the difference vector. Larger values for F result in higher diversity in the generated population, and lower values cause faster convergence. According to the study, the value of F is often set to 0.5 [24], [45]. Step 2: Crossover. After the mutation process, to increase the potential diversity of the population, a crossover operation is applied to each pair of the target vector Xk,G and its corresponding mutant vector Vk,G to generate a trial vector Uk,G+1 = {u1k,G , u2k,G , . . . , uD k,G }. In the basic version, DE employs the binomial (uniform) crossover, defined as follows: utk,G+1 =

⎧ t ⎨ vk,G , ⎩ xt

k,G ,

if (randt [0, 1] ≤ CR) or (t = trand ), t = 1, 2, . . . , D otherwise (8)

where randt [0, 1] is a uniform random number generator and the crossover rate CR is a user-specified constant within the range [0, 1], which controls the fraction of parameter values copied from the mutant vector. In [24] and [45], the experiential value of CR is set to 0.9, and trand is a randomly chosen integer in the range [1, D] to ensure that at least one parameter of Uk,G+1 is selected from the mutated vector Vk,G . Step 3: Selection. This approach must decide which vector (Uk,G+1 or Xk,G ) should be a member of the next (new) generation G + 1. For a minimization problem, the vector with the lower fitness value is chosen by a greedy selection scheme, as follows:  Xk,G+1 =

Uk,G+1 , if f (Uk,G+1 ) ≤ f (Xk,G ) otherwise. Xk,G ,

(9)

Step 4: Stopping condition. Once the new population is obtained, the process of mutation, crossover, and selection is repeated until the optimal solution is located or a prespecified termination criterion is satisfied, e.g., the number of generations reaches a preset maximum Gmax . Compared with most other evolution algorithms, DE is much simpler and significantly more straightforward to implement. The main body of the algorithm takes four to five lines to code in any programming language. Due to its simplicity of coding, DE has been applied to various engineering problems, including numerical optimization [25], [27]–[29] mechanical engineering [30], feedforward neural network training [31], digital filter design [32], image processing [33], [34], and pattern recognition [35], [36].

IV. ADESM A LGORITHM The subpixel mapping problem should be considered as a permutation-based combinational problem in the discrete space, which allocates a land-cover class to the position of subpixels (see also Section II-B and Fig. 3). However, the canonical DE algorithm cannot be used to directly generate a discrete subpixel mapping problem since it was originally designed to solve continuous optimization problems. Therefore, in this section, we will present an enhanced DE approach, namely, an ADESM algorithm, to solve the subpixel mapping problem. In ADESM, a reversible-conversion strategy is used to transform the individual encoding of a subpixel mapping solution from integer to real numbers. That is, the discrete subpixel mapping individuals are changed to real encoding. By this transformation, traditional DE operators for continuous space may be retained and directly used, with the advantage of DE. In addition, a self-adapting parameter selection method is proposed to adaptively choose the appropriate parameters during the course of ADESM. To describe ADESM, the following notations are used. 1) Let Z denote the input coarse fraction image, where Z = {zmn : 1 ≤ m ≤ M, 1 ≤ n ≤ N } and M and N represent the row number and column number of the image, respectively. 2) Let zmn denote a pixel in the coarse fraction image i i : i = 1, . . . , C}, zmn denotes the Z, where zmn = {zmn fraction for the ith land-cover class, and C is the number of the land-cover classes or endmembers. 3) Let s represent the scale of the subpixel mapping and D represent the number of the subpixels in a pixel. Here, D is equal to s2 . ADESM consists of the following steps for each pixel zmn .

A. Calculating the Number of Subpixels for Each Land-Cover Class i In the first step, the number of subpixels in land-cover class i, N Ci , is derived according to the following equation:   i  , N Ci = round D ∗ zmn C 

i = 1, . . . , C

i zmn =1

i=1 C 

N Ci = D

(10)

i=1

where round() denotes the operator that rounds its argument toward the closer integer. Based on (10), the sum of N Ci values should be equal to D. By Step A, the number of subpixels for each land-cover class is calculated. The next process utilizes DE to obtain the optimal positions of subpixels for each land-cover class within the pixel.

ZHONG AND ZHANG: REMOTE SENSING IMAGE SUBPIXEL MAPPING BASED ON ADAPTIVE DE

1311

To construct the mathematical model, choice variables yij are defined such that [5] Minimize J = − f = − 

C  D 

1,

if subpixel j is assigned to land-cover class i otherwise

yij = Fig. 4.

0,

Initial individual.

B. Initial Population and Discrete Encoding After ADESM obtains the number of subpixels for each land-cover class i, N Ci , it needs to generate the initial population PG = {X1,G , . . . , Xk,G , . . . , XN P,G }, where Xk,G represents the kth individual in the Gth generation and N P is the size of the population. Each individual Xk,G = {x1k,G , . . . , xtk,G .., xD k,G } denotes a candidate subpixel mapping solution. Each individual xtk,G in ADESM is discretely encoded to describe the position of the subpixel with the corresponding class attribute, where xtk,G is unique and systematic, as shown in Fig. 4. For example, for Fig. 2(a), where N C1 = 6 and N C2 = 3, six and three subpixels belong to classes 1 and 2, respectively. According to the aforementioned encoding principle, for the candidate solution shown in Fig. 2(b), Xk,G is encoded to {1, 2, 4, 5, 7, 8, 3, 6, 9} (see also Fig. 3). That is, the first, second, fourth, fifth, seventh, and eighth subpixels are allocated to class 1, and other subpixels are allocated to class 2. The initial population PG=0 is obtained according to the aforementioned discrete encoding using (11) and should cover the entire search space as much as possible by uniformly randomizing individuals within the search space constrained by the prescribed minimum and maximum parameter bounds Xmin = {x1min , . . . , xtmin , . . . , xD min } and Xmax = t t {x1max , . . . , xtmax , . . . , xD max }, where xmin and xmax represent the bounds of the tth position. For the subpixel mapping problem, xtmin and xtmax are equal to one and D, respectively      xtk,G=0 = (int) randt (0, 1) · xtmax + 1 − xtmin + xtmin , t = 1, . . . , D;

xtk ∈ x1k , . . . , xt−1 k

(11)

where xtk,G=0 represents the tth position in the kth individual Xk,G .

C. Calculation of the Objective Function Using an SDI After Step B, the objective function of the new solution is calculated. It should be noted that the purpose of solving the subpixel mapping problem is to obtain the maximum spatial dependence using (3) while the purpose of DE is to obtain the minimum value of the objective function. Therefore, we must change the form of the objective function (3) using minimize J in (12).

yij · SDIij

i=1 j=1

(12)

where D is the number of subpixels within a pixel, C is the number of land-cover classes, and SDIij is the measure for spatial dependence assigned to land-cover class i in the jth subpixels. The value of yij is obtained according to xtk,G . For example, assuming that N C1 = 6 and N C2 = 3, a candidate subpixel mapping solution {1, 2, 4, 5, 7, 8, 3, 6, 9} in Figs. 3 and 4 ensures that the first, second, fourth, fifth, seventh, and eighth subpixels are allocated to class 1 y1j = {1¯1 , 1¯2 , 0, 1¯4 , 1¯5 , 0, 1¯7 , 1¯8 , 0} and other subpixels are allocated to class 2 y2j = {0, 0, 1¯3 , 0, 0, 1¯6 , 0, 0, 1¯9 }. The fitness of the individual Xk,G is calculated as follows: J(Xk,G ) = − [(SDI1,1 +SDI1,2 +SDI1,4 +SDI1,5 +SDI1,7 + SDI1,8 ) + (SDI2,3 + SDI2,6 + SDI2,9 )] . (13) The SDI can be calculated by (13). Unlike (2), the SDIij is the average SDI, not the sum of the SDI values. The reason is that each pixel has the nine neighboring pixels in the original version, while the number of neighboring pixels will change in the boundary pixel in the proposed algorithm. To better compare the SDI at the same scale, we use the average SDI value to replace the sum of the SDI values N n

SDIij =

wk · F ractioni,k

k=1

Nn

(14)

where wk is the weight, Nn is the number of the neighboring pixels, and F ractioni,k represents the fraction value of class i for the kth neighboring pixel around the subpixel j. It is worth noting that the value of the weight wk is inversely proportional to the Euclidean distance (dis) between the subpixels and the neighboring pixels. In this paper, the weight is equal to the reciprocal of the distance by (15), and the distance is calculated by (16). The shorter the distance, the higher the weight. In this step, it is necessary to project the coordinate of the subpixels (g, h) and the pixels (m, n) to the same coordinate systems using (17) and (18) and adaptively select the neighboring pixels to avoid losing boundary information. Coordinate Conversion: The aforementioned process is shown in Fig. 5 by a simple example with a scale of four. The original coordinate of the pixel (m, n) is (0.5, 0.5) in the XY coordinate system, and the coordinate of the subpixel (g, h) is (1.5, 0.5) in the xy coordinate system. To calculate the distance between subpixels and pixels, the new Euclidean coordinate system, with vertical X-axis and horizontal Y -axis, is defined. In the new coordinate system, the pixel (m, n) with (0.5, 0.5) is converted to (m , n ), (2, 2), using (17), and the subpixel (g, h) with (1.5, 0.5) is converted to (g  , h ), (5.5, 4.5), using (18).

1312

IEEE TRANSACTIONS ON SYSTEMS, MAN, AND CYBERNETICS—PART B: CYBERNETICS, VOL. 42, NO. 5, OCTOBER 2012

Fig. 6.

Self-adapting encoding strategy.

D. Real Transformation To apply DE to solve the subpixel mapping problem, ADESM transforms the subpixel mapping solution Xk,G = {x1k,G , . . . , xtk,G , . . . , xD k,G } in the discrete space to the continuous real space. In this paper, if the precision of the value to be generated is set to a defined decimal place a, an integer number xtk,G is transformed by the following equation: xtk,G = Fig. 5. Coordinate conversion and the distance calculation between the pixels and subpixels.

According to (15) and (16), the weight of the neighboring pixel w in the upper left corner of the subpixels is equal to 0.233

xtk,G 10−a

(21)

where the range of the variable xtk,G is between 10a and D∗ 10a . For example, when a is equal to two, then the value of the individuals in Fig. 2 is within [100, 900]. E. Adaptive Mutation and Crossover

w = 1/d

d = (m − g  )2 + (n − h )2 (m , n ) ← (s × m), s × n) (g  , h ) ← (s + g, s + n).

(15) (16) (17) (18)

Adaptively Selecting Neighboring Pixels: The traditional neighboring pixel zi,j for the pixel zm,n using the surrounding neighborhood [20] is selected as follows: N [zm,n ] = {zi,j |(i = m ∨ j = n) ∧ [(m − 1 ≤ i ≤ m + 1) ∧(n − 1 ≤ j ≤ n + 1)] ∧ [(1 < m < M ) ∧ (1 < n < N )]} . (19) In (19), the boundary pixels zm,n (m = 1, M or n = 1, N ) are not considered, to guarantee that each pixel has nine neighboring pixels. Thus, in some previous methods, such as BP subpixel mapping (BPSM) and the subpixel mapping algorithm proposed by Ge et al. [22], the subpixel mapping results do not contain the boundary pixels and lose some subpixel information. To avoid this problem, the proposed method improves (19) to adaptively select valid neighboring pixels, as in the following equation: N [zm,n ] = {zi,j |(i = m ∨ j = n) ∧ [(max(1, m − 1) ≤ i ≤ min(m + 1, M )) ∧ (max(1, n − 1) ≤ j ≤ min(n + 1, N ) ] ∧ [(1 ≤ m ≤ M ) ∧ (1 ≤ n ≤ N )]} (20) where M and N represent the row number and column number, respectively. The functions max() and min() return the largest and smallest values from the numbers provided, respectively.

By Step D, the population changes from integer to real number. These values may be used in the traditional DE mutation and crossover operators. There are two control parameters that need to be adjusted by the user: the scaling factor F and the crossover rate CR, according to (7) and (8), respectively. Suitable control parameters are different for different real problems. In some cases, the time for finding these appropriate parameters can be unacceptably long, particularly for large-volume remote sensing images. To solve the problem, a self-adaptive strategy for the control parameters is used. The control parameters F and CR are encoded to each individual, as shown in Fig. 6. That is, each individual has its corresponding F and CR and can be adjusted during the process of evolution [40]. The subpixel mapping solution (Fig. 6) is represented by the D-dimensional vector Xk,G and two control parameters Fk,G and CRk,G in the Gth generation, k = 1, . . . , N P . For each transformed vector Xk,G at the generation G, its 1 2 D , vk,G , . . . , vk,G } can associated mutant vector Vk,G = {vk,G be generated via the strategy DE/rand/1/bin, which is the strategy most often used in practice [24], [40], [44]. The mutation operators are as follows: Vk,G = Xr1 ,G +F · (Xr2 ,G −Xr3 ,G ),

k = 1,. . ., N P

(22) where the indices r1 , r2 , and r3 are mutually exclusive integers randomly generated within the range [1, N P ], r1 = r2 = r3 = k. The better objective values in (12) of these encoded control parameters are more likely to survive and produce offspring, which leads to better individuals and improves the probability of finding the optimal solution. To adaptively determine the

ZHONG AND ZHANG: REMOTE SENSING IMAGE SUBPIXEL MAPPING BASED ON ADAPTIVE DE

mutation rate pm according to the objective values of each individual, the process is as follows: pm =

J(Xk,G ) − min (J(Xk,G )) . max (J(Xk,G )) − min (J(Xk,G ))

(23)

1313

the solutions are within the bounds. If a solution is detected to have violated a bound, it is dragged to the offending boundary  t xmin , if utk,G < xtmin (28) utk,G = xtmax , if utk,G > xtmax

New control parameters in the G + 1th generation Fk,G+1 and CRk,G+1 are updated as follows, with probability pm :  b G (1− Gmax ) 1 − rand , if rand2 < pm Fk,G+1 = (24) 1 Fk,G , otherwise

where xtmin and xtmax represent the bounds of the tth position. For the subpixel mapping problem, xtmin and xtmax are equal to one and D, respectively.

where randt , t ∈ {1, 2}, denotes the uniform random values within the range [0, 1], G is the iteration number, Gmax is the maximal iteration number, b is a parameter to decide the nonconforming degree, and the experiential value is set to three [46]. In the initial generations, ADESM tends to search the solution space uniformly, and in the later generations, it tends to search the space locally, i.e., closer to its descendants, to guarantee convergence. After the mutation phase, a crossover operation is applied to generate a trial vector Uk,G = {u1k,G , u2k,G , . . . , uD k,G } for the mutant vector Vk,G as follows:  t vk,G , if (randt [0, 1] ≤ CR) or (t = trand ) utk,G = t otherwise, xk,G ,

Although the value utk,G in Uk,G ranges from xtmin to xtmax , the discrete individual of unique values cannot be guaranteed. To solve this problem, a repair operator is added. To repair the individual, we need to record the positions of the repeating values and the missing values, denoted by RV = {rvi } and M V = {mvi }, respectively. The position of the repeating values rvi is selected randomly, and the rvi th value is replaced by the missing value mvi . Since each value is randomly selected, the value has to be removed from the array after selection, to avoid duplication.

t = 1, 2, . . . , D;

k = 1, . . . , N P.

CRk,G+1 is updated using the following [40]:  if rand3 < pm rand4 , CRk,G+1 = CRk,G , otherwise

(25)

(26)

where randt , t ∈ {3, 4}, denotes the uniform random values within the range [0, 1]. Fk,G+1 and CRk,G+1 are obtained before the mutation is performed. Therefore, they influence the mutation, crossover, and selection operations of the new vector Xk,G+1 .

H. Enhanced Operators The traditional DE algorithm should be able to obtain a satisfactory result by mutation and crossover operators because it operates in real continuous space. However, to solve the discrete subpixel mapping problem using the aforementioned process, ADESM may still be unable to find an appropriate solution, due to some replicated solutions. In this section, enhanced operators (exchange and insertion) are added to extend the search space and improve the quality of the individuals. Exchange Operator: The exchange operator, as an improvement technique, may explore random regions in the searching space to find a better subpixel mapping solution by simply exchanging two values in an individual [41]. Two unique random values are selected: r1 and r2 ∈ [1, D], where r1 = r2 . The values indexed by these values are exchanged as exchange

F. Integer Transformation After the mutation and crossover processes, ADESM obtains the set of the solution Uk,G = {u1k,G , u2k,G , . . . , uD k,G } [see also (25)] in continuous space. The data type of these solutions is floating. To find the possible subpixel mapping solution in the discrete space, integer transformation is used to convert the real value back to the integer, as given in (27), assuming utk,G to be the real value obtained from (25)  utk,G  . int utk,G = 10a

G. Repair

(27)

The value utk,G is rounded to the nearest integer using an inverse process of (21). Real transformation and integer transformation effectively allow DE to optimize permutative solutions for discrete subpixel mapping problems. Once the set of the solution Uk,G = {u1k,G , u2k,G , . . . , uD k,G } is obtained after the transformation, it is then checked for feasibility. Feasibility, measured using (28), only refers to whether

1 2 ↔ urk,G , and the individual is evaluated. If the oburk,G jective value improves, then the new solution is accepted in the population. Insertion Operator: Insertion may provide greater diversity to the solution than the standard mutation in discrete space, which can be regarded as a more complicated mutation method [41]. In this step, two unique random numbers are selected: r1 1 and r2 ∈ [1, D]. The r1 th value urk,G is removed to behind the r2 other value uk,G . To explain the process from Step G to Step H, a simple example is shown in Fig. 7. For example, assuming an infeasible subpixel mapping solution obtained by integer transformation, Uk,G = {1, 3, 8, 1, 5, 7, 2, 9, 9}, D = 9. 1) Repair: As shown, the values 1 and 9 are repeated in the solution, so the positions of the repeated values are recorded: rv1 = 4, and rv2 = 8. We can find that the values 4 and 6 are missing: mv1 = 4, and mv2 = 6. ADESM randomly selects the missing value from M V to the located rvi th position in Uk,G . It assumes

1314

IEEE TRANSACTIONS ON SYSTEMS, MAN, AND CYBERNETICS—PART B: CYBERNETICS, VOL. 42, NO. 5, OCTOBER 2012

Fig. 8.

Local search process.

Fig. 7. Simple example using repair and enhanced operators.

K. Stopping Condition that the missing values 4 and 6 are arranged to the fourth and eighth positions, respectively, and Uk,G = {1, 3, 8, 4, 5, 7, 2, 6, 9}. 2) Exchange: Two random numbers are generated within the bounds: Rnd = {2, 7}. u2k,G and u7k,G in Uk,G are exchanged to obtain the new solution Uk,G = {1, 2, 8, 4, 5, 7, 3, 6, 9}. 3) Insertion: Two random numbers are generated within the bounds: Rnd = {3, 6}. Then, u3k,G is removed to behind u6k,G , Uk,G = {1, 2, 4, 5, 7, 8, 3, 6, 9}.

I. Selection After the calculation of the objective function using (12), a selection operation is performed. The objective function value of each trial vector J(Uk,G ) is compared with that of its corresponding target vector J(Xk,G ) in the current population. If the trial vector has a lower or an equal objective function value, compared with the corresponding target vector, the trial vector will replace the target vector and enter the population of the next generation. Otherwise, the target vector will remain in the population for the next generation. The selection operation can be expressed as follows:  Uk,G , Xk,G+1 = Xk,G ,

if J(Uk,G) ≤ J(Xk,G ), k = 1, 2,. . ., N P otherwise. (29)

J. Local Search There is always a possibility of stagnation in DE [41]. To avoid stagnation, a local search operator is used to find a more feasible solution in the local neighborhood. The basic process of a local search is shown in Fig. 8. The set of the selected number A is set to empty in the initialization. On each iteration, a random number I is chosen, with its range within [1, D]. Another random number j is chosen in another set B. The values in the individual indexed by i and j are then swapped. The objective function of the new individual is calculated, and only if there is an improvement is the new solution accepted.

If the generation G does not meet the maximum generation number Gmax , go to Step C. Otherwise, output the best individuals as the subpixels of the pixel zmn . The aforementioned process is repeated from Step A until the subpixels of all pixels in the image Z are located. Finally, the proposed algorithm outputs the subpixel mapping image. The flowchart for the ADESM framework is shown in Fig. 1. V. E XPERIMENTS AND A NALYSIS The aforementioned ADESM was coded in Visual C++ 6.0 and tested on different images. To demonstrate the effectiveness of the proposed approach, three different types of data sets were used to test its performance. Consistent comparisons between ADESM and the previous subpixel mapping algorithms—the DNSM algorithm, the subpixel mapping algorithm based on spatial-attraction models (SASM), and the BP neural network subpixel mapping algorithm (BPSM)—were carried out. To compare ADESM with other evolutionary algorithms (EAs), the subpixel mapping algorithm based on a genetic algorithm [genetic-algorithm subpixel mapping (GASM)] was also used in all the experiments. The surrounding neighboring type was used for all subpixel mapping algorithms, and the estimations of the subpixel mapping accuracy for the several subpixel mapping algorithms are provided. A. Experiment 1: Synthetic Artificial Images In experiments 1 and 2, synthetic imagery was used. Synthetic images are images that have been degraded to a coarser scale using an averaging filter. Degradation of a hard classification yields fraction images for each class in the classification. These images can be interpreted as a soft classification. Synthetic imagery has the advantage of lacking coregistration errors between the lower and higher resolution images. The resulting soft classification does not contain any uncertainty as it originates from degradation instead of a classification process. Consequently, the subpixel mapping error solely reflects the performance of the proposed method. Validation is facilitated, as the original hard classification can be used as reference material [14].

ZHONG AND ZHANG: REMOTE SENSING IMAGE SUBPIXEL MAPPING BASED ON ADAPTIVE DE

1315

Fig. 9. Subpixel mapping results in experiment 1 (scale = 4). (a) Original artificial image (400 × 400). (b) Original hard classification result (400 × 400). (c) Fraction images (100 × 100). (d) DNSM. (e) SASM. (f) BPSM. (g) GASM. (h) ADESM. (i) Zoom images for all the subpixel mapping results.

In this experiment, synthetic artificial imagery was used to enable the refinement of the technique. The synthetic image was created by the authors to aid in the design and development of the algorithm. A visual checking of the algorithm’s performance is easy with the synthetic image as it represents simple geometric figures. To better simulate the real remote sensing scene in the artificial image, we designed four landcover classes—water, tree, agricultural field, and impervious layer (e.g., building and road)—to illustrate the functionality of the technique, as shown in Fig. 9(a). A maximum-likelihood hard classification algorithm was applied to the synthetic artificial images to provide real imagery to be used for degradation and accuracy testing, which is shown in Fig. 9(b). Fig. 9(c) shows the degraded fraction image with a scale of four from a hard classification image. Subpixel mapping using the DNSM, SASM, BPSM, GASM, and ADESM algorithms is shown in Fig. 9(d)–(h). For a clear comparison of ADESM with the other algorithms, two small images of the S1 and S2 areas are zoomed and are shown in Fig. 9(i).

In this experiment, ADESM could adaptively choose the most user-defined running parameters, as described in Section IV-E. The number of iterations was equal to 100. The configuration for BPSM (using one hidden layer) was as follows: The number of hidden layers, the learning rate, the momentum rate, and the number of iterations were 1, 0.25, 0.9, and 1000, respectively. The crossover rate and mutation rate for GASM were set to 0.5 and 0.05. From the comparison between Fig. 9(b) and Fig. 9(d)–(h), it can be observed that ADESM shows better subpixel mapping results than DNSM, BPSM, SASM, and GASM do. Serrated edges exist between different classes in the subpixel mapping images obtained by DNSM and BPSM, with many subpixels of the image allocated to incorrect positions, as shown in Fig. 9(i), S1. It is worth noting that BPSM obtains better subpixel mapping for two ridges of the field on the middle of the results than the other three subpixel mapping algorithms from Fig. 9(i), S2. However, BPSM does not obtain the subpixel result of the edge, which is noted by the black pixel in Fig. 9(f). In general, SASM,

1316

IEEE TRANSACTIONS ON SYSTEMS, MAN, AND CYBERNETICS—PART B: CYBERNETICS, VOL. 42, NO. 5, OCTOBER 2012

TABLE I A DJUSTED C ONFUSION M ATRICES OF THE S UBPIXEL M APPING R ESULTS D ERIVED F ROM THE S YNTHETIC A RTIFICIAL I MAGE IN E XPERIMENT 1

TABLE II C OMPARISON OF F IVE S UBPIXEL M APPING A LGORITHM P ERFORMANCES FOR THE S YNTHETIC A RTIFICIAL I MAGE IN E XPERIMENT 1

GASM, and ADESM can obtain smoother visual results for all classes. For a more detailed verification of the results, we compared hard classification with the subpixel mapping results and assessed the accuracy of each algorithm quantitatively, using the overall accuracy (OA) measure and Kappa coefficient based on the confusion matrices [47], as shown in Tables I and II. In addition, an adjusted confusion matrix was utilized. The adjusted confusion matrix is identical to the original confusion matrix except that it is calculated only for mixed pixels, which means that it ignores the subpixels that have a pure pixel as

parent. These subpixels all belong to the same class and may raise the Kappa coefficient without providing information about the algorithm’s prediction abilities [9]. Based on the adjusted confusion matrix, we calculated the adjusted OA (OA∗ ), the adjusted Kappa coefficient (Kappa∗ ), the average of the producer’s accuracy, the average of the user’s accuracy, and the average of Short’s mapping accuracy index (Short’s index) [48]. Among them, the average of the producer’s accuracy, the average of the user’s accuracy, OA, and Kappa value are widely used in the validation of land-use/land-cover classification [4], [48]. The average of Short’s mapping accuracy index is the

ZHONG AND ZHANG: REMOTE SENSING IMAGE SUBPIXEL MAPPING BASED ON ADAPTIVE DE

arithmetic mean of this index [49], [50] (a monotonic function of the harmonic mean of the user’s and producer’s accuracies) that explicitly combines both user’s and producer’s accuracies in one measure [4]. This measure is supposed to be supplementary to the average of the user’s accuracy and the average of the producer’s accuracy. Table I lists the confusion matrices of comparisons between the hard reference classification image and the subpixel mapping results produced by five subpixel mapping algorithms: DNSM, SASM, BPSM, GASM, and ADESM. It is worth noting that the number of total subpixels for BPSM is less than that for the other subpixel mapping algorithms because BPSM did not place, on the subpixel map, the edge of the lowresolution image, as shown in Fig. 9(f). For a more detailed verification of the results, we assessed the accuracy of each method using the producer’s and user’s accuracy measures (see Table I). As an example, for C1 simulated water class, the producer’s accuracies of DNSM, SASM, BPSM, GASM, and ADESM are 73.38%, 96.69%, 82.62%, 94.88%, and 96.91%, respectively. The producer’s accuracies of these methods show that some of the C1 pixels were wrongly identified as C2 class. A similar situation occurred for C2 and C4 classes. A careful evaluation of the confusion matrix reveals that there is confusion when discriminating among C1, C2, and C4 classes in DNSM, SASM, GASM, and BPSM. It can be observed that there is also significant confusion when discriminating the C4 class from the other classes, particularly when using DNSM, SASM, GASM, and BPSM. ADESM recognizes the C4 class better than the other methods. It also exhibits the highest producer’s accuracy for C1, C2, and C4 classes among the five methods considered. It is worth noting that BPSM has obtained the highest producer’s and user’s accuracies for C3 class, i.e., 89.45%. However, it is difficult for BPSM to recognize the other three classes. Although the producer’s and user’s accuracies of SASM are higher than those of DNSM and BPSM, the producer’s and user’s accuracies provided by ADESM are further increased, to some degree, than SASM. Table II shows the global performance in terms of the subpixel mapping accuracy yielded by DNSM, SASM, BPSM, GASM, and the proposed ADESM. From the table, we can observe that ADESM exhibits the highest OA, Kappa value, adjusted OA, adjusted Kappa, and average of Short’s mapping accuracy index. DNSM has the lowest subpixel mapping accuracy because it did not consider the distance of the neighboring pixels. BPSM did not obtain a satisfactory result, as there is a sawtooth on the edges. Moreover, BPSM still needs an additional training data set. SASM utilized the spatialattraction models, improving the results of DNSM by considering neighboring information. However, employing SASM did not guarantee obtaining the optimal subpixel mapping result. To obtain the optimal result in GASM and ADESM, the subpixel mapping problem is transformed into an optimization problem of the SDI. ADESM, as a new subpixel mapping algorithm based on a DE algorithm, like GASM, belongs to the family of evolutionary computation and is a population-based stochastic global optimizer. However, there are some major philosophical differences between GASM and ADESM. Although ADESM also uses selection, crossover, and mutation operators, the op-

1317

erations need to be redefined. In ADESM, candidate subpixel mapping solutions are represented by chromosomes based on discrete numbers. In the mutation process of an ADESM algorithm, the weighted difference between two randomly selected population members is added to a third member to generate a mutated solution. Then, a crossover operator follows to combine the mutated solution with the target subpixel mapping solution, to generate a trial solution. Thereafter, a selection operator is applied to compare the fitness function values of both competing solutions, namely, target and trial solutions, to determine which can survive for the next generation. Following the aforementioned process, ADESM obtains the optimal subpixel mapping results for the synthetic artificial image in the experiment. B. Experiment 2: Synthetic Images—TM and HJ-1A Images Synthetic imagery is also used for real remote sensing images. This helps in avoiding the uncertainty inherent in real imagery that is caused by the sensor point spread function, atmospheric and geometric effects, and spectral unmixing or classification errors. Work in progress is aimed at applying the technique to real imagery and reducing such uncertainty. Two different types of remote sensing images, namely, Landsat TM images and Chinese HJ-1A images, were chosen to demonstrate the same spatial dependence present among features in both images at different resolutions. For ADESM algorithms, we carried out the experiments with the number of iterations equal to 100, and the other three parameters, i.e., the number of individuals in the population N P , the scaling factor F , and the crossover rate CR, were adaptively obtained. The number of hidden layers, the learning rate, the momentum rate, and the number of iterations were fixed to 1, 0.25, 0.9, and 1000, respectively, in the BPSM method. The crossover rate and mutation rate for GASM were set to 0.5 and 0.05. First, we tested the subpixel mapping algorithm proposed in experiment 2 using the 30-m resolution multispectral Landsat TM image shown in Fig. 10(a) [51]. The image (400 × 400 pixels) was acquired in Wuhan City, Hubei Province, China, on October 26, 1998. The survey area is part of the city, and the primary objective image was expected to fall into four classes: water, vegetation, road/bare land, and building. The classification image shown in Fig. 10(b) was obtained by an unsupervised artificial immune classifier [51]. Fig. 10(c) shows the degraded fraction images with a scale of four from a hard classification image. The subpixel mapping results, using DNSM, SASM, BPSM, GASM, and ADESM algorithms, are shown in Fig. 10(d)–(h). To ensure clarity in comparing ADESM with the other algorithms, two small images of S1 and S2 areas are zoomed and are shown in Fig. 10(i). To evaluate the subpixel mapping, the original hard classification result shown in Fig. 10(b) is considered as the real subpixel mapping result. As shown in Fig. 10 [specifically, Fig. 10(i) zoom images for all subpixel mapping results], SASM, GASM, and ADESM obtained satisfactory visual subpixel mapping results, which means that their results were smoother and more continuous than those of DNSM and BPSM. When comparing ADESM with SASM and GASM in

1318

IEEE TRANSACTIONS ON SYSTEMS, MAN, AND CYBERNETICS—PART B: CYBERNETICS, VOL. 42, NO. 5, OCTOBER 2012

Fig. 10. Subpixel mapping results for the Wuhan TM image in experiment 2 (scale = 4). (a) Wuhan TM image red-green-blue (RGB) (4, 3, 2) (400 × 400). (b) Original hard classification result (400 × 400). (c) Fraction images (100 × 100). (d) DNSM. (e) SASM. (f) BPSM. (g) GASM. (h) ADESM. (i) Zoom images for all subpixel mapping results.

Fig. 10(h), S2, one can see that ADESM can retain the linear features from a coarse image, e.g., a road (white) on the right of the S2 images. The adjusted confusion matrices and global performance indices obtained by applying DNSM, SASM, BPSM, GASM, and ADESM to this data set are shown in Tables III and IV. As an example, one can see from Table III that, for the road/bareland class, the producer’s accuracies of DNSM and BPSM are 49.42% and 54.28%, respectively. They are characterized by a significant amount of confusion between the road/bare-land class and the other classes. Many road/bare-land subpixels are misclassified as building or vegetation classes, as can be seen in Fig. 10. In addition, the subpixel mapping results obtained by DNSM and BPSM are fuzzy and not smooth, as a whole. By contrast, SASM and ADESM have better visual subpixel mapping results, as shown in Fig. 10(e) and (h), compared with Fig. 10(b). According to Table III, SASM and ADESM

also exhibited a higher accuracy (i.e., 76.76% and 77.38% for the water class, respectively) than the other subpixel mapping techniques. On the other hand, there is confusion when discriminating vegetation from building, and many vegetation subpixels are incorrectly allocated to the building class by SASM and GASM. It can also be observed that SASM and GASM misclassified many building subpixels to the vegetation class. In Table IV, we can observe that the proposed ADESM obtained the highest OA and Kappa value, i.e., 77.54% and 0.6911, respectively. ADESM also has the highest adjusted OA and Kappa value, compared with the other subpixel mapping results. These results are also confirmed by the other quality indices considered, including the average of Short’s mapping accuracy index. The reason for this may be attributed to the fact that ADESM can obtain the maximum of the SDI by using evolution operators, mutation, crossover, repair, etc. In

ZHONG AND ZHANG: REMOTE SENSING IMAGE SUBPIXEL MAPPING BASED ON ADAPTIVE DE

1319

TABLE III A DJUSTED C ONFUSION M ATRICES OF THE S UBPIXEL M APPING R ESULTS D ERIVED F ROM THE W UHAN TM I MAGE IN E XPERIMENT 2

TABLE IV C OMPARISON OF F IVE S UBPIXEL M APPING A LGORITHM P ERFORMANCES FOR THE W UHAN TM I MAGE IN E XPERIMENT 2

addition, ADESM uses the local search strategy to improve the connectivity property. It may also add the probability to find the optimal subpixel mapping solution. Hence, with this data set, the ADESM technique is superior to the other three algorithms considered. Another synthetic image, based on the HJ-1A image, was used to test the performance of the proposed subpixel mapping algorithm by comparing it with the DNSM, BPSM, SASM, and GASM algorithms. The HJ-1A satellite is a small Chinese environmental satellite. It was launched on September 6, 2008,

from the Taiyuan Satellite Launch Center of Shanxi Province [52]. It has a sun-synchronous circular orbit with an orbital altitude of 649 km. The HJ-1A satellite has a hyperspectral sensor with 115 bands, including blue, green, red, and nearinfrared spectra. It has a spatial resolution of 100 m and a spectral range of 0.45–0.95 μm. At present, the satellite imaging area can cover large parts of Asia, specifically, more than ten countries, including China, India, Pakistan, Kazakhstan, Mongolia, South Korea, North Korea, Japan, the Philippines, and Thailand [53]. The HJ-1A image (256 × 256), shown in

1320

IEEE TRANSACTIONS ON SYSTEMS, MAN, AND CYBERNETICS—PART B: CYBERNETICS, VOL. 42, NO. 5, OCTOBER 2012

Fig. 11. Subpixel mapping results for the Hanchuan HJ-1A image in experiment 2 (scale = 4). (a) Hanchuan HJ-1A image RGB (95, 73, 43) (256 × 256). (b) Original hard classification result (256 × 256). (c) Fraction images (64 × 64). (d) DNSM. (e) SASM. (f) BPSM. (g) GASM (h) ADESM. (i) Zoom images for all subpixel mapping results.

Fig. 11(a), was acquired on August 19, 2009 (path: 01; raw: 79), and was used as the original image for the experimental results. The study site is located in Hanchuan City, Hubei Province, central China, and its surrounding area. Four land-cover classes, i.e., city, agricultural land, water, and vegetation, characterize this image. Fig. 11(b) shows the hard classification result for the Hanchuan HJ-1A image obtained by a support vector machine (SVM) in ENVI software [54]. The fraction images of the four classes were obtained, as shown in Fig. 11(c), by resampling with a scale of four. Fig. 11(d)–(h) shows the subpixel mapping results obtained by using DNSM, SASM, BPSM, GASM, and the proposed algorithm ADESM (100 iterations). To quantitatively evaluate the subpixel mapping accuracy, the original hard classification result, as the real subpixel mapping result, was used to test the performance of the five algorithms. To allow us to evaluate all the subpixel mapping algorithms from the visual

results, two small images of the S1 and S2 areas are zoomed and are shown in Fig. 11(i). As shown in Fig. 11(d)–(i), DNSM and BPSM did not provide satisfactory visual results, as there is blurring and many structural features are lost. In addition, the sawtooth phenomenon can be observed in the BPSM result. By contrast, SASM, GASM, and ADESM obtained better visual results, being smoother and with most classes’ structural information preserved. Quantitative comparisons of the five algorithms, based on the adjusted confusion matrices and global performance indices obtained by applying DNSM, SASM, BPSM, GASM, and ADESM to this data set, are shown in Tables V and VI. As with the visual results, the subpixel mapping accuracies of DNSM and BPSM are lower than those of the other two algorithms. For instance, one can see from Table V that the average producer’s accuracies of DNSM and BPSM

ZHONG AND ZHANG: REMOTE SENSING IMAGE SUBPIXEL MAPPING BASED ON ADAPTIVE DE

are equal to 61.34% and 68.72%, respectively. SASM and GASM improve the subpixel mapping accuracy and locate the subpixels’ position better than the aforementioned algorithms, resulting in producer’s accuracies of 72.89% and 72.61%. It was ADESM, however, that obtained the highest subpixel mapping result, having a producer’s accuracy of 74.12%. This is because the proposed algorithm directly takes into account the spatial dependence as the optimization problem, to guarantee the connectivity property of subpixels and the smoothness of the result. In Table V, it can be seen that ADESM obtains the highest OA, Kappa coefficient, adjusted OA, and adjusted Kappa coefficient, which are 79.83%, 0.7263, 74.04%, and 0.6479, respectively. These results are also confirmed by the other quality indices considered, e.g., the average of Short’s mapping accuracy index being 0.59, thus pointing out the superiority of ADESM over the other four algorithms for this data set. C. Experiment 3: Real Image The synthetic imagery used in experiments 1 and 2 lacks coregistration errors, which have been shown to be an important source of errors in classification or spectral unmixing. Neither do the fraction images contain any uncertainty, as they originate from degradation instead of a classification or spectral unmixing process. A comprehensive discussion on the production of synthetic imagery is provided by Mertens et al. [12]. Although synthetic images may avoid the reflection of other errors (i.e., a coregistration error), the subpixel mapping methods were not directly applied to subpixel mapping for real imagery. In the third experiment, we used the subpixel mapping methods to obtain a higher resolution classification image from a real lowresolution remote sensing image (i.e., 30-m ETM image). To evaluate the result, a real high-resolution classification image based on a 4-m IKONOS image of the same research area was used to test the performance of these subpixel mapping algorithms. Although there are other errors, such as coregistration errors and spectral unmixing errors, these errors are the same for all the subpixel mapping methods, and relative comparisons are fair. The Landsat ETM standard product (98 × 30) over the Shenzhen City zone of China, shown in Fig. 12(a), was acquired on December 31, 2001. This product was georeferenced with a spatial resolution of 30 m. The IKONOS data over the same region, with a 4-m spatial resolution, shown in Fig. 12(b), were acquired on December 20, 2001 [55], and were used to obtain a high-resolution classification image as the evaluated standard of the subpixel mapping results. Although an 11-day difference exists between the two images, both of them were acquired in the same season; hence, their reflectance spectra are assumed to be similar. In such a case, the spectral signature may be obtained through inverse spectral mixture analysis to derive endmember signatures or training pixels, given the fractional covers that can be obtained from the IKONOS image. As a result, we first performed a geometric correction on the image so that each corrected pixel of the output image has the same Universal Transverse Mercator coordinate as the ETM data. The ETM spatial resolution was resampled to 32 m, which

1321

is eight times the IKONOS data (4 m). Three steps were included to calculate the fractional land cover for each ETM pixel. First, the IKONOS image was classified into four classes (city, vegetation, bare soil, and water) using an SVM with field surveying, assisted by interpretation, as the shade reflection intertwines with water. The classification image is shown in Fig. 12(c). Recent studies have found that the typical scale of urban reflectance is between 10 and 20 m [56]. Therefore, most 4 × 4 m IKONOS pixels are spectrally homogeneous and can be reasonably assumed to be pure materials. Second, the classified IKONOS image was registered with the ETM images using the ground control points with a total root mean square error of less than 0.03. Third, the fraction images of the ETM image were obtained by the least squares linear spectral mixture analysis method [3], using ENVI remote sensing image processing software [54]. The fraction images are shown in Fig. 12(d). Fourth, the subpixel mapping results using different methods—DNSM, SASM, BPSM, GASM, and ADESM—were acquired, based on the fraction image with scale = 8. The subpixel mapping results are shown in Fig. 12(e)–(i). Finally, the classification map of the IKONOS image was overlaid on the subpixel mapping results to estimate the performance of the different subpixel mapping algorithms at the IKONOS data scale. To ensure a clear comparison of ADESM with the other algorithms, two small images of S1 and S2 areas are zoomed and are shown in Fig. 12(j). In the subpixel mapping process, the number of iterations for ADESM and GASM was set to 100. The main parameters of BPSM were set as follows: number of hidden layers = 1, learning rate = 0.25, momentum rate = 0.9, and number of iterations = 1000. In addition, BPSM needs training samples to train the BP neural network. In this experiment, approximately 50 samples for each class were selected randomly as the training samples for BPSM. The hard classification result of the IKONOS image, as shown in Fig. 12(c), as the true result, was used to test the performance of the different subpixel mapping algorithms. As shown in Fig. 12, when comparing the subpixel mapping results of DNSM, SASM, BPSM, GASM, ADESM, and the highresolution IKONOS classification image, DNSM and BPSM did not provide satisfactory subpixel mapping results. For example, the structural information of the water class has been lost on the right of the images, and the structure of the playground is very fuzzy, as shown in Fig. 12(j), S1. The whole visual result of GASM is fuzzy, particularly for the vegetation and bare soil classes. Although SASM showed better visual results than DNSM, GASM, and BPSM and provided a more integrated subpixel mapping result, the continuity of the vegetation class was inadequate on the left of the results, as shown in Fig. 12(j), S2, in which many vegetation subpixels were cut. By contrast, ADESM had better visual results and provided more integrated and continuous subpixel mapping results, e.g., for the vegetation and water classes. To quantitatively evaluate the subpixel mapping accuracy, the adjusted confusion matrices and global performance indices were obtained by applying DNSM, SASM, BPSM, GASM, and ADESM to this data set (see Tables VII and VIII). It is worth noting that BPSM did not map the boundary of the image

1322

IEEE TRANSACTIONS ON SYSTEMS, MAN, AND CYBERNETICS—PART B: CYBERNETICS, VOL. 42, NO. 5, OCTOBER 2012

TABLE V A DJUSTED C ONFUSION M ATRICES OF THE S UBPIXEL M APPING R ESULTS D ERIVED F ROM THE H ANCHUAN HJ-1A I MAGE IN E XPERIMENT 2

TABLE VI C OMPARISON OF F IVE S UBPIXEL M APPING A LGORITHM P ERFORMANCES FOR THE H ANCHUAN HJ-1A I MAGE IN E XPERIMENT 2

in the subpixel layer, for example, the black lines shown in Fig. 12(g). The subpixels in the boundary were not considered, to test the performance of the BPSM subpixel mapping result. The number of real subpixels for different classes in the BPSM image is less than that in the subpixel mapping images of the other subpixel mapping algorithms. BPSM obtained the highest subpixel mapping producer’s accuracies of the vegetation and water classes, i.e., 22.28% and 21.57%, respectively. However, it was difficult for BPSM to differentiate the city class from the other classes. In addition, one can see from Table VII that, for the vegetation class, the producer’s accuracies of DNSM

and SASM are equal to 21.66% and 21.85%, respectively, and better than that of ADESM. However, this class is significantly confused with city and water, i.e., only 26.99% and 27.17%, respectively, of the vegetation subpixels are subpixel mapped correctly. By contrast, ADESM recognized the vegetation class better than the aforementioned algorithms, resulting in the user’s accuracy of 33.83%. Moreover, ADESM had the highest subpixel accuracy for the city class, with the producer’s and user’s accuracies of 76.71% and 64.19%, respectively. In Table VIII, we can observe that, in general, ADESM obtained the highest adjusted OA (OA∗ ) and adjusted Kappa

ZHONG AND ZHANG: REMOTE SENSING IMAGE SUBPIXEL MAPPING BASED ON ADAPTIVE DE

1323

Fig. 12. Subpixel mapping results for the Shenzhen TM image in experiment 3 (scale = 8). (a) Shenzhen ETM image RGB (4, 3, 2) (98 × 30). (b) Shenzhen IKONOS image (784 × 240). (c) Shenzhen IKONOS classification image. (d) Fraction images of ETM image (98 × 30). (e) DNSM. (f) SASM. (g) BPSM. (h)GASM. (i) ADESM. (j) Zoom images for all the subpixel mapping results.

(Kappa∗ ) value, with gains of 2.97%, 2.81%, and 2.08% over DNSM, SASM, and BPSM for OA∗ , respectively. It is worth noting that BPSM provides a higher subpixel mapping result than DNSM and SASM do for difficult real remote sensing images. Comparing ADESM with GASM, the two algorithms have similar adjusted OA, i.e., 53.89% and 54.21%, but the visual result of GASM is poorer than ADESM and the subpixel mapping result of GASM is fuzzy, as shown in Fig. 12(h) and (j). In addition, we can find that the subpixel accuracies using the real data are lower than those for the synthetic data for all algorithms. The possible reasons for this are as follows: 1) Georeferenced and geometric corrected errors may exist; 2) the spectral unmixing error of the ETM image and the classification error of the IKONOS image may influence the

evaluated accuracies; 3) there are radiometric correction and reflectance retrieval between IKONOS and ETM; and 4) the difference among homogeneous spectra is large. It can therefore be stated that the quality of the input soft classification or spectral unmixing result is a key factor to the success of the subpixel mapping. These results are also confirmed by the other quality indices considered, thus pointing to the superiority of ADESM over the other four algorithms for this data set, particularly the data set of the real image. VI. S ENSITIVITY A NALYSIS OF ADESM The scale of subpixel mapping is an important parameter for ADESM. In Section V, we only analyzed the subpixel mapping

1324

IEEE TRANSACTIONS ON SYSTEMS, MAN, AND CYBERNETICS—PART B: CYBERNETICS, VOL. 42, NO. 5, OCTOBER 2012

TABLE VII A DJUSTED C ONFUSION M ATRICES OF THE S UBPIXEL M APPING R ESULTS FOR THE S HENZHEN ETM I MAGE IN E XPERIMENT 3

TABLE VIII C OMPARISON OF F IVE S UBPIXEL M APPING A LGORITHM P ERFORMANCES FOR THE S HENZHEN ETM I MAGE IN E XPERIMENT 3

result when the scale was equal to four. In this section, the sensitivity analysis between scale and five subpixel mapping algorithms (DNSM, SASM, BPSM, GASM, and ADESM) was carried out to evaluate the performance of these algorithms. In ADESM, we added exchange, insertion, and local search operators to improve the subpixel mapping algorithm. To investigate the effectiveness of these operators, the sensitivity analysis between the improved operators and ADESM was undertaken. In addition, we analyzed the effectiveness of the self-adaptive strategy for ADESM. The proof of convergence of EAs with self-adaptation is difficult because control parameters

are changed randomly and the selection does not affect their evolution directly [40], [57], [58]. Since DE is a particular instance of EA, it is interesting to investigate how self-adaptivity can be applied to it [59]. In this section, to test the validity of the self-adaptive strategy in ADESM, we analyzed the sensitivity of three control parameters for the original DE subpixel mapping (DESM) algorithm: amplification factor of the difference vector F , crossover control parameter CR, and population size N P . In this section, the synthetic artificial image and two synthetic real images were tested using different parameter values.

ZHONG AND ZHANG: REMOTE SENSING IMAGE SUBPIXEL MAPPING BASED ON ADAPTIVE DE

1325

Fig. 13. Sensitivity of subpixel mapping algorithms in relation to s. (a) Artificial image. (b) Wuhan TM image. (c) Hanchuan HJ-1A image.

Fig. 14. Sensitivity of DESM in relation to improved operators. (a) Artificial image. (b) Wuhan TM image. (c) Hanchuan HJ-1A image.

A. Sensitivity in Relation to Parameter s To study the ADESM sensitivity in relation to scale s, the other parameters were kept the same as those in experiments 1 and 2. The values of the scale were assumed as follows: s = {2, 4, 8, 16}. As shown in Fig. 13, the higher the s value, the lower the adjusted OA (OA∗ ). Overall, the values of OA∗ , for all the subpixel mapping algorithms proposed in this paper, take on the tendency of descension. This is because the distribution of land-cover classes is more complex when the scale s increases, resulting in an increase in the difficulty of subpixel mapping. Although the OA∗ of the artificial image and two synthetic images decreases as the scale increases, those of the subpixel mapping images, overall, are superior than those of the degraded images and provide more subpixel information. Compared with the other subpixel mapping algorithms, ADESM is a more stable subpixel mapping algorithm, which can obtain better subpixel mapping results than DNSM, GASM, and BPSM. When the scale s is equal to two and four, ADESM has a higher accuracy than SASM. For scales 8 and 16, ADESM and SASM have similar subpixel mapping accuracy.

B. Sensitivity in Relation to Improved Operators To investigate the effectiveness of the improved operators, including exchange, insertion, and local search operators, we removed these operators of ADESM one by one. The other parameters were kept the same as those in experiments 1 and 2.

To denote these algorithms, the following notations are used. ADESM-Exchange represents the ADESM algorithm without exchange operator, ADESM-Insertion is the ADESM algorithm without insertion operator, ADESM-Exchange + Insertion represents the ADESM without exchange and insertion operators, and ADESM-Localsearch describes the ADESM algorithm without local search operators. In addition, to better evaluate the performance of ADESM, the subpixel mapping algorithms based on other adaptive-DE algorithms, for example, JADE (an adaptive differential evolution proposed by J. Zhang and A. C. Sanderson [38]), namely, the subpixel mapping algorithm based on JADE (JADESM), were used to compare with ADESM. As shown in Fig. 14, JADESM obtained the best OA∗ values, i.e., 90.72% and 67.23%, for experiments 1 and 2, respectively, and ADESM-Localsearch obtained the best OA∗ , i.e., 74.14%, for experiment 3. It was noted that the original JADE also cannot be used in the discrete space. Thus, the discrete-encoding and reversible-conversion strategies of ADESM were utilized in JADE for subpixel mapping. Although JADESM and ADESMLocalsearch slightly increased the subpixel mapping accuracy, their computational costs are too much, being approximately five to six times greater than ADESM. That is, the local search operator of ADESM can decrease the computational cost and find a better feasible solution in the local neighborhood. Compared with ADESM-Exchange, ADESM-Insertion, and ADESM-Exchange + Insertion, ADESM is a more stable subpixel mapping algorithm, which gave appropriate subpixel mapping results, with satisfactory computation times of 80.32, 386.85, and 163.85 s for the three images.

1326

IEEE TRANSACTIONS ON SYSTEMS, MAN, AND CYBERNETICS—PART B: CYBERNETICS, VOL. 42, NO. 5, OCTOBER 2012

Fig. 15. Sensitivity of DESM in relation to N P . (a) Artificial image. (b) Wuhan TM image. (c) Hanchuan HJ-1A image.

Fig. 16. Sensitivity of DESM in relation to CR. (a) Artificial image. (b) Wuhan TM image. (c) Hanchuan HJ-1A image.

C. Sensitivity in Relation to Parameter N P To compare our self-adaptive version of the DESM algorithm, i.e., ADESM, with the original DESM algorithm, the best control parameter setting for DESM may be needed. The DESM algorithm does not change the control parameter values during the subpixel mapping process, except for the analyzed parameter. The experiential parameters of the original DESM algorithm were set as follows: CR = 0.9 and F = 0.5, in [24] and [45]. The number of the initial population N P is very important in maintaining the diversity of the population and extending the search range in the feature space. To analyze the sensitivity in relation to parameter N P , the other parameters did not change, and N P assumed the following values for three experimental images: N P = {10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160}. Fig. 15 shows the sensitivity of DESM in relation to the parameter N P by comparing it with the adjusted OA provided by the DESM algorithms. The value of N P is of importance to the DESM algorithm in maintaining the diversity of the population, because N P determines the population size used to extend the search space and gives the mutation a good chance to proceed toward the solution. As shown in Fig. 15, there was an upward trend in the adjusted OA (OA∗ ) of the subpixel mapping algorithm based on the original DE algorithm (DESM) when the value of N P was changed from 10 to 160. When N P was equal to 120, 140, and 160, for the three experimental images, the highest OA∗ values of DESM are 90.73%, 66.93%, and 73.74%, respectively. It was noted that, the higher the N P value, the greater the computational cost of DESM, because

the population size is increased. When the N P is equal to 160, DESM has computational costs of 92.37, 444.88, and 188.42 s for the artificial image, Wuhan TM image, and Hanchuan HJ1A image, respectively. Compared with DESM, ADESM, with N P = 160, has computational costs of 80.32, 386.85, and 163.85 s for the three experiments. In a comparison of the OA∗ value with ADESM, ADESM obtained higher OA∗ values, i.e., 67.06%, and 74.04%, for the two synthetic real images.

D. Sensitivity in Relation to Parameter CR For all three experimental images, the DESM algorithm was performed with CR taken from [0.1, 1.0] by steps of 0.1, while the other parameters were set as follows: N P = 100, F = 0.5, and s = 4. First, we set the control parameter CR = 0.1 and kept it fixed during 30 independent runs. Then, we set CR = 0.2 for the next 30 runs, etc. The results were averaged over 30 independent runs. The experimental results for different images are shown in Fig. 16. The best adjusted OA (OA∗ ) values of DESM for the artificial image and Wuhan TM image, i.e., 90.74% and 67.06%, respectively, were obtained by CR = 0.4 (N P = 100, F = 0.5, and s = 4). The values are slightly higher than or equal to those of ADESM, i.e., 90.60% and 67.06%. For the Hanchuan HJ-1A image, ADESM obtained a better OA∗ value, i.e., 74.04%, than the best OA∗ value of DESM, i.e., 73.97%, when CR was equal to 0.5. Although DESM can obtain satisfactory results by adjusting the value of parameter CR, ADESM may adaptively provide similar or better subpixel mapping results, without prior knowledge or experience.

ZHONG AND ZHANG: REMOTE SENSING IMAGE SUBPIXEL MAPPING BASED ON ADAPTIVE DE

1327

Fig. 17. Sensitivity of DESM in relation to F . (a) Artificial image. (b) Wuhan TM image. (c) Hanchuan HJ-1A image.

E. Sensitivity in Relation to Parameter F According to (7), F is a key parameter for the DESM algorithm. Hence, we further investigated the effect of F on solution quality. In accordance with the conclusion from Storn and Price’s results, where F is in the range of 0.4–0.9 and CR is in the range of 0.5–0.9 [24], we set F from 0.1 to 1.0 with steps equal to 0.1 and fixed the other parameters with N P = 100 and CR = 0.9. Similar computational experiments were then conducted. The experimental results are presented for three different images in Fig. 17. From Fig. 17, it follows that, as F increases, the adjusted OA (OA∗ ) produced by the DESM algorithm decreases. The best adjusted overall accuracies of DESM, i.e., 90.72%, 66.88%, and 73.93%, for the three different images, were obtained when F was equal to 0.1. As shown in Fig. 17; although DESM may obtain higher OA∗ than ADESM for a synthetic artificial image, ADESM obtained better results than the DESM algorithm and did not need any other prior knowledge for the two synthetic real images. Based on the aforementioned sensitivity analysis, there are two disadvantages in the original DESM. First, parameter tuning requires multiple runs, and it is usually not a feasible solution for problems that are very time consuming. Second, the best control parameter settings of DESM are problem dependent. The proposed ADESM overcomes these disadvantages, so there is no need for multiple runs to adjust the control parameters; moreover, self-adaptive DESM is much more problem independent than the original DESM is. Therefore, we conclude that ADESM is a robust and effective subpixel mapping algorithm.

VII. C ONCLUSION Based on the DE theory, this paper has proposed a new mapping strategy for the subpixel mapping of remote sensing images. In line with this strategy, the subpixel mapping problem is transformed to an optimization problem in the discrete subpixels’ space by maximizing the SDI. The traditional DE algorithm only obtains the optimal solution in the continuous space, it does not process the subpixel mapping problem in the discrete space. In addition, DE needs to choose the proper control parameters, employing prior experience of the users, population size N P , crossover rate CR, and scaling factor F .

This is quite a difficult task because the best settings for the control parameters are not easy to determine for complex problems. In this paper, to apply DE to the subpixel mapping problem effectively, an improved DE approach, the ADESM algorithm, has been proposed. The proposed method is a combination of a discrete binary DE approach with an adaptive strategy, which is enhanced by advanced evolution operators and local search. It utilizes real transformation and integer transformation to encode candidate subpixel mapping solutions from the discrete space. During the subpixel mapping process, the enhanced DE operators—repair, exchange, and insertion—are added to improve the subpixel mapping results, along with the traditional evolution operators—mutation and crossover. The local search strategy, as an improvement strategy, is also used to find more feasible solutions in the local neighborhood. In addition, the proposed self-adaptive method is an attempt to determine the optimum values of control parameters F and CR. Our self-adaptive DESM algorithm (ADESM) has been implemented and tested on subpixel mapping problems using different remote sensing images taken from the literature. The experimental results in this paper consistently show that the ADESM algorithm, with its self-adaptive control parameter settings, provides better results than the traditional DNSM and SASM algorithms, as well as the latest BP neural network subpixel mapping algorithm from the literature. In the three experiments using different images (a synthetic artificial image, two synthetic remote sensing images, and a real remote sensing image), ADESM was able to achieve better results, with a higher subpixel mapping accuracy. This evinces that ADESM is appropriate for the subpixel mapping of remote sensing images and is an improved alternative subpixel mapping algorithm. In the comparison of the subpixel mapping accuracies between synthetic images and real images, ADESM accomplishes satisfactory subpixel mapping results for three synthetic images, with an accuracy of more than 67%. With the real data of ETM and IKONOS, the proposed method also achieves a higher accuracy than the other algorithms, but the accuracy is lower than that for the synthetic images because there are geometric correction, classification, or spectral unmixing errors. The sensitivity analysis of the scale demonstrates that the subpixel mapping accuracy decreases when the value of the scale is increased from 2 to 16. To evaluate the validity of the self-adaptive strategy in ADESM, we analyzed the sensitivity in relation to three parameters, namely, N P ,

1328

IEEE TRANSACTIONS ON SYSTEMS, MAN, AND CYBERNETICS—PART B: CYBERNETICS, VOL. 42, NO. 5, OCTOBER 2012

CR, and F , by using different user-defined values. The results show that our algorithm, with its self-adaptive control parameter settings, is better than, or at least comparable to, the standard DESM algorithms. However, our self-adaptive method decreases the user-defined process and can be simply used to solve the subpixel mapping problem with a satisfactory level of accuracy. Although we tried to apply ADESM to a real remote sensing image, and obtained better subpixel mapping results than traditional algorithms do, it is still difficult to process a real remote sensing image with many errors, e.g., spectral unmixing errors. Our future work will explore further methodologies and aims to find a significantly better data set to improve the subpixel mapping accuracy for real applications, thereby investigating the adaptability of the SDI for a complex real image (e.g., hyperspectral image [60]).

[13]

[14] [15]

[16] [17]

[18]

ACKNOWLEDGMENT The authors would like to thank the Editor-in-Chief, the Associate Editor, and the three anonymous reviewers for the helpful comments and suggestions that improved this paper. Dr. Y. Zhong would also like to thank Prof. K. Tang of the University of Science and Technology of China and Dr. A. M. Zhou of East China of Normal University for their helpful discussions and constructive suggestions.

[19] [20] [21] [22] [23]

R EFERENCES [1] A. Boucher and P. C. Kyriakidis, “Super-resolution land cover mapping with indicator geostatistics,” Remote Sens. Environ., vol. 104, no. 3, pp. 264–282, Oct. 2006. [2] C.-I. Chang, Hyperspectral Imaging: Spectral Detection and Classification.. New York: Plenum, 2003. [3] D. C. Heinz and C.-I. Chang, “Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery,” IEEE Trans. Geosci. Remote Sens., vol. 39, no. 3, pp. 529–545, Mar. 2001. [4] C. Yang, L. Bruzzone, F. Sun, L. Lu, R. Guan, and Y. Liang, “A fuzzystatistics-based affinity propagation technique for clustering in multispectral images,” IEEE Trans. Geosci. Remote Sens., vol. 48, no. 6, pp. 2647– 2659, Jun. 2010. [5] J. Verhoeye and R. De Wulf, “Land cover mapping at sub-pixel scales using linear optimization techniques,” Remote Sens. Environ., vol. 79, no. 1, pp. 96–104, Jan. 2002. [6] P. M. Atkinson, “Mapping sub-pixel boundaries from remotely sensed images,” in Proc. Innov. GIS, Z. Kemp, Ed., 1997, vol. 4, pp. 166–180. [7] P. M. Atkinson, E. Pardo-lguzquiza, and M. Chica-Olmo, “Downscaling cokriging for super-resolution mapping of continua in remotely sensed images,” IEEE Trans. Geosci. Remote Sens., vol. 46, no. 2, pp. 573–580, Jan. 2008. [8] M. W. Thornton, P. M. Atkinson, and D. A. Holland, “Sub-pixel mapping of rural land cover objects from fine spatial resolution satellite sensor imagery using super-resolution pixel swapping,” Int. J. Remote Sens., vol. 27, no. 3, pp. 473–491, 2006. [9] K. C. Mertens, L. P. C. Verbeke, T. Westra, and R. R. D. Wulf, “Using genetic algorithms in sub-pixel mapping,” Int. J. Remote Sens., vol. 24, no. 21, pp. 4241–4247, 2003. [10] P. M. Atkinson, “Sub-pixel target mapping from soft-classified, remotely sensed imagery,” Photogramm. Eng. Remote Sens., vol. 71, no. 7, pp. 839– 846, Jul. 2005. [11] M. W. Thornton, P. M. Atkinson, and D. A. Holland, “A linearised pixelswapping method for mapping rural linear land cover features from spatial resolution remotely sensed imagery,” Comput. Geosci., vol. 33, no. 10, pp. 1261–1272, Oct. 2007. [12] A. J. Tatem, H. G. Lewis, P. M. Atkinson, and M. S. Nixon, “Superresolution target identification from remotely sensed images using a Hop-

[24] [25] [26] [27] [28] [29] [30] [31]

[32] [33] [34] [35]

field neural network,” IEEE Trans. Geosci. Remote Sens., vol. 39, no. 4, pp. 781–796, Apr. 2001. A. J. Tatem, H. G. Lewis, P. M. Atkinson, and M. S. Nixon, “Increasing the spatial resolution of agricultural land cover maps using a Hopfield neural network,” Int. J. Geograph. Inf. Sci., vol. 17, no. 7, pp. 647–672, Oct./Nov. 2003. A. J. Tatem, H. G. Lewis, P. M. Atkinson, and M. S. Nixon, “Superresolution land cover pattern prediction using a Hopfield neural network,” Remote Sens. Environ., vol. 79, no. 1, pp. 1–14, 2002. K. C. Mertens, L. P. C. Verbeke, and R. R. D. Wulf, “Sub-pixel mapping with neural networks: Real-world spatial configurations learned from artificial shapes,” in Proc. 4th Int. Symp. Remote Sens. Urban Areas, ISPRS Int. Archives Photogramm., Remote Sens. Spatial Inf. Sci., Regensburg, Germany, Jun. 27–29, 2003, vol. 34-7/W9, pp. 117–121. L. Zhang, K. Wu, Y. Zhong, and P. Li, “A new sub-pixel mapping algorithm based on a BP neural network with an observation model,” Neurocomputing, vol. 71, no. 10–12, pp. 2046–2054, Jun. 2008. K. C. Mertens, L. P. C. Verbeke, T. Westra, and R. R. D. Wulf, “Subpixel mapping and sub-pixel sharpening using neural network predicted wavelet coefficients,” Remote Sens. Environ., vol. 91, no. 2, pp. 225–236, May 2004. A. Boucher, P. C. Kyriakidis, and C. Cronkite-Ratcliff, “Geostatistical solutions for super-resolution land cover mapping,” IEEE Trans. Geosci. Remote Sens., vol. 46, no. 1, pp. 272–283, Jan. 2008. K. C. Mertens, B. D. Baets, L. P. C. Verbeke, and R. R. D. Wulf, “Direct sub-pixel mapping exploiting spatial dependence,” in Proc. IEEE Int. Geosci. Remote Sens. Symp., 2004, pp. 3046–3049. K. C. Mertens, B. D. Baets, L. P. C. Verbeke, and R. R. D. Wulf, “A sub-pixel mapping algorithm based on sub-pixel/pixel spatial attraction models,” Int. J. Remote Sens., vol. 27, no. 15, pp. 3293–3310, Aug. 2006. T. Kasetkasem, M. K. Arora, and P. K. Varshney, “Super-resolution land cover mapping using a Markov random field based approach,” Remote Sens. Environ., vol. 96, no. 3/4, pp. 302–314, Jun. 2005. Y. Ge, S. Li, and V. C. Lakhan, “Development and testing of a subpixel mapping algorithm,” IEEE Trans. Geosci. Remote Sens., vol. 47, no. 7, pp. 2155–2164, Jul. 2009. Z. Q. Shen, J. G. Qi, and K. Wang, “Modification of pixel-swapping algorithm with initialization from a sub-pixel/pixel spatial attraction model,” Photogramm. Eng. Remote Sens., vol. 75, no. 5, pp. 557–568, May 2009. R. Storn and K. V. Price, “Differential Evolution—A simple and efficient heuristic for global optimization over continuous spaces,” J. Global Optim., vol. 11, no. 4, pp. 341–359, Dec. 1997. A. K. Qin, V. L. Huang, and P. N. Suganthan, “Differential evolution algorithm with strategy adaptation for global numerical optimization,” IEEE Trans. Evol. Comput., vol. 13, no. 2, pp. 398–417, Apr. 2009. L. Wang, Q. Pan, P. N. Suganthan, W. Wang, and Y. Wang, “A novel hybrid discrete differential evolution algorithm for blocking flow shop scheduling problems,” Comput. Oper. Res., vol. 37, no. 3, pp. 509–520, Mar. 2010. W. Gong, Z. Cai, C. X. King, and H. Li, “Enhanced differential evolution with adaptive strategies for numerical optimization,” IEEE Trans. Syst. Man, Cybern. B, Cybern., vol. 41, no. 2, pp. 397–413, Apr. 2011. S. Das, A. Konar, and U. K. Chakraborty, “Two improved differential evolution schemes for faster global search,” in Proc. GECCO, 2005, vol. 5, pp. 991–998. S. Das, A. Abraham, U. K. Chakraborty, and A. Konar, “Differential evolution using a neighborhood-based mutation operator,” IEEE Trans. Evol. Comput., vol. 13, no. 3, pp. 526–553, Jun. 2009. R. Joshi and A. C. Sanderson, “Minimal representation multisensor fusion using differential evolutions,” IEEE Trans. Syst. Man Cybern. A, Syst. Humans, vol. 29, no. 1, pp. 63–76, Jan. 1999. C. H. Chen, C. J. Lin, and C. T. Lin, “Nonlinear system control using adaptive neural fuzzy networks based on a modified differential evolution,” IEEE Trans. Syst. Man Cybern. C, Appl. Rev., vol. 39, no. 4, pp. 459–473, Jul. 2009. R. Storn, “Designing digital filters with differential evolution,” in New Ideas in Optimization, D. Corne, M. Dorigo, and F. Glover, Eds. London, U.K.: McGraw-Hill, 1999. S. Das and S. Sil, “Kernel-induced fuzzy clustering of image pixels with an improved differential evolution algorithm,” Inf. Sci., vol. 180, no. 8, pp. 1237–1256, Apr. 2010. S. Das and A. Konar, “Automatic image pixel clustering with an improved differential evolution,” Appl. Soft Comput., vol. 9, no. 1, pp. 226–236, Jan. 2009. S. Das, A. Abraham, and A. Konar, “Automatic clustering using an improved differential evolution algorithm,” IEEE Trans. Syst. Man Cybern. A, Syst. Humans, vol. 38, no. 1, pp. 218–237, Jan. 2008.

ZHONG AND ZHANG: REMOTE SENSING IMAGE SUBPIXEL MAPPING BASED ON ADAPTIVE DE

[36] U. Maulik and I. Saha, “Automatic fuzzy clustering using modified differential evolution for image classification,” IEEE Trans. Geosci. Remote Sens., vol. 48, no. 9, pp. 3503–3510, Sep. 2010. [37] X. Yuan, A. Su, H. Nie, Y. Yuan, and L. Wang, “Application of enhanced discrete differential evolution approach to unit commitment problem,” Energ. Convers. Manage., vol. 50, no. 9, pp. 2449–2456, Sep. 2009. [38] J. Zhang and A. C. Sanderson, “JADE: Adaptive differential evolution with optional external archive,” IEEE Trans. Evol. Comput., vol. 13, no. 5, pp. 945–958, Oct. 2009. [39] A. Ghosh, S. Das, A. Chowdhury, and R. Giri, “An improved differential evolution algorithm with fitness-based adaptation of the control parameters,” Inf. Sci., vol. 181, no. 18, pp. 3749–3765, Sep. 2011. [40] J. Brest, S. Greiner, B. Boškovi, M. Mernik, and V. umer, “Self-adapting control parameters in differential evolution: A comparative study on numerical benchmark problems,” IEEE Trans. Evol. Comput., vol. 10, no. 6, pp. 646–657, Dec. 2006. [41] D. Davendra and G. Onwubolu, “Forward backward transformation,” in Differential Evolution: A Handbook for Global Permutation-Based Combinatorial Optimization, G. Onwubolu and D. Davendra, Eds. Berlin, Germany: Springer-Verlag, 2009, pp. 35–80. [42] N. Noman and H. Iba, “Accelerating differential evolution using an adaptive local search,” IEEE Trans. Evol. Comput., vol. 12, no. 1, pp. 107–125, Feb. 2008. [43] S. Das and P. N. Suganthan, “Differential evolution: A survey of the state-of-the-art,” IEEE Trans. Evol. Comput., vol. 15, no. 1, pp. 4–31, Feb. 2011. [44] K. Price and R. Storn, Differential Evolution: A Practical Approach to Global Optimization. Berlin, Germany: Springer-Verlag, 2005. [45] S. Rahnamayan, H. R. Tizhoosh, and M. A. Salama, “Opposition-based differential evolution,” IEEE Trans. Evol. Comput., vol. 12, no. 1, pp. 64– 79, Feb. 2008. [46] L. Zhang, Y. Zhong, B. Huang, and P. Li, “Dimensionality reduction based on clonal selection for hyperspectral imagery,” IEEE Trans. Geosci. Remote Sens., vol. 45, no. 12, pp. 4172–4185, Dec. 2007. [47] J. A. Richards and X. P. Jia, Remote Sensing Digital Image Analysis, An Introduction., 4th ed. Berlin, Germany: Springer-Verlag, 2006. [48] C. Liu, P. Frazier, and L. Kumar, “Comparative assessment of the measures of thematic classification accuracy,” Photogramm. Eng. Remote Sens., vol. 107, no. 4, pp. 606–616, Apr. 2007. [49] N. M. Short, The Landsat Tutorial Workbook Basics of Satellite Remote Sensing. Greenbelt, MD: Goddard Space Flight Center, 1982. [50] G. H. Rosenfield and K. Fitzpatrick-Lins, “A coefficient of agreement as a measure of thematic classification accuracy,” Photogramm. Eng. Remote Sens., vol. 52, no. 2, pp. 223–227, 1986. [51] Y. Zhong, L. Zhang, B. Huang, and P. Li, “An unsupervised artificial immune classifier for multi/hyper-spectral remote sensing image,” IEEE Trans. Geosci. Remote Sens., vol. 44, no. 2, pp. 420–431, Feb. 2006. [52] G. X. Bai, “China environmental and disaster monitoring and forecasting small satellite—HJ-1A/B,” Aerosp. China, vol. 5, pp. 10–15, 2009. [53] S. Lu, B. Wu, N. Yan, and H. Wang, “Water body mapping method with HJ-1A/B satellite imagery,” Int. J. Appl. Earth Obs. Geoinf., vol. 13, no. 3, pp. 428–434, Jun. 2011, DOI: 10.1016/j.jag.2010.09.006. [54] American ITT Visual Information Solutions Company, ENVI Online Tutorials [EB/OL]. [Online]. Available: http://www.ittvis.com/ ProductServices/ENVI.aspx [55] L. Zhang, B. Wu, B. Huang, and P. Li, “Nonlinear estimation of subpixel proportion via kernel least square regression,” Int. J. Remote Sens., vol. 28, no. 18, pp. 4157–4172, 2007. [56] C. Small, “High spatial resolution spectral mixture analysis of urban reflectance,” Remote Sens. Environ., vol. 88, no. 1/2, pp. 170–186, 2003. [57] M. A. Semenov and D. A. Terkel, “Analysis of convergence of an evolutionary algorithm with self-adaption using a stochastic Lyapunov function,” Evol. Comput., vol. 11, no. 4, pp. 363–379, Winter 2003. [58] J. He and X. Yao, “Toward an analytic framework for analysing the computation time of evolutionary algorithms,” Artif. Intell., vol. 145, no. 1/2, pp. 59–97, Apr. 2003.

1329

[59] J. Liu and J. Lampinen, “A fuzzy adaptive differential evolution algorithm,” Soft Comput.—Fusion Found., Methodol. Appl., vol. 9, no. 6, pp. 448–462, Jun. 2005. [60] A. Plaza, J. A. Benediktsson, J. W. Boardman, J. Brazile, L. Bruzzone, G. Camps-Valls, J. Chanussot, M. Fauvel, P. Gamba, A. Gualtieri, M. Marconcini, J. C. Tilton, and G. Trianni, “Recent advances in techniques for hyperspectral image processing,” Remote Sens. Environ., vol. 113, no. Supplement 1, pp. S110–S122, Sep. 2009.

Yanfei Zhong (M’11) received the B.S. degree in information engineering and the Ph.D. degree in photogrammetry and remote sensing from Wuhan University, Wuhan, China, in 2002 and 2007, respectively. Since 2007, he has been with the State Key Laboratory of Information Engineering in Surveying, Mapping and Remote Sensing, Wuhan University, where he is currently a Professor. He was a Referee of Pattern Recognition. He has published more than ten peer-reviewed articles in international journals, such as IEEE T RANSACTIONS ON G EOSCIENCE AND R EMOTE S ENSING and International Journal of Remote Sensing. His research interests include multiand hyperspectral remote sensing image processing, artificial intelligence, and pattern recognition. Dr. Zhong was the recipient of the National Excellent Doctoral Dissertation Award of China (in 2009) and the New Century Excellent Talents in University of China (in 2009). He was a Referee of IEEE T RANSACTIONS ON S YSTEMS , M AN AND C YBERNETICS—PART B and IEEE J OURNAL OF S ELECTED T OPICS IN A PPLIED E ARTH O BSERVATIONS AND R EMOTE S ENSING.

Liangpei Zhang (M’06–SM’08) received the B.S. degree in physics from Hunan Normal University, Changsha, China, in 1982, the M.S. degree in optics from Xi’an Institute of Optics and Precision Mechanics, Chinese Academy of Sciences, Xi’an, China, in 1988, and the Ph.D. degree in photogrammetry and remote sensing from Wuhan University, Wuhan, China, in 1998. He is currently with the State Key Laboratory of Information Engineering in Surveying, Mapping and Remote Sensing, Wuhan University, as the Head of the Remote Sensing Division. He is also a “Chang-Jiang Scholar” Chair Professor appointed by the Ministry of Education, China. He is currently the Principal Scientist for the China State Key Basic Research Project (2011–2016) appointed by the Ministry of National Science and Technology of China to lead the remote sensing program in China. He is an Executive Member (Board of Governor) of the Chinese National Committee of International Geosphere–Biosphere Programme. He also serves as an Associate Editor of International Journal of Ambient Computing and Intelligence, International Journal of Image and Graphics, International Journal of Digital Multimedia Broadcasting, Journal of Geo-spatial Information Science, and Journal of Remote Sensing. He has more than 230 research papers and is the holder of 5 patents. His research interests include hyperspectral remote sensing, highresolution remote sensing, image processing, and artificial intelligence. Dr. Zhang is a fellow of the Institution of Electrical Engineers, an executive member of China Society of Image and Graphics, and others. He regularly serves as a Cochair of the series SPIE Conferences on Multispectral Image Processing and Pattern Recognition, Conference on Asia Remote Sensing, and many other conferences. He edits several conference proceedings, issues, and the Geoinformatics Symposiums.