extfieldsofstudy:
Computer Science
Mathematics

text:
Some Improvements in Fuzzy Turing Machines

In this paper, we modify some previous definitions of fuzzy Turing machines to define the notions of accepting and rejecting degrees of inputs, computationally. We use a BFS-based search method and obtain an upper level bound to guarantee the existence of accepting and rejecting degrees. We show that fuzzy, generalized fuzzy and classical Turing machines have the same computational power. Next, we introduce the class of Extended Fuzzy Turing Machines equipped with indeterminacy states. These machines are used to catch some types of loops of the classical Turing machines. Moreover, to each r.e. or co-r.e language, we correspond a fuzzy language which is indeterminable by an extended fuzzy Turing machine.

Introduction
In [14], Lotfi Zadeh defines the notion of fuzzy algorithm. His definition is based on a fuzzification of Turing machines. However, that work was not deep enough in the recursion theoretical aspects of the mentioned model. That work is followed by the same setting in [6]. The equivalency of previous fuzzy models is shown in [8,9]. Afterwards, the research in this field is revisited in [2,3,5]. In [2], Biacino and Gerla generalize the definition of recursive enumerability introduced in [5]. Next, Wiedermann proposes a formal fuzzy computing model based on Turing machine model [12,13]. He claims that it is possible to accept r.e. sets and co-r.e. sets by fuzzy Turing machines and so these machines can solve the halting problem. So he claims that these fuzzy Turing machines have more computational power than the classical Turing machines. In [1], Bedregal and Figueira analyse Wiedermann's statement about the computational power of fuzzy Turing machines and show that Wiedermann's statement is not completely correct. They give a characterization of the class of n-r.e. sets in terms of associated fuzzy languages recognized by fuzzy Turing machines. They also show that there is no universal fuzzy Turing machine which can simulate each machine in the class of all fuzzy Turing machines. More recently, Moniri defines the class of Generalized Fuzzy Turing Machines [7]. His machines are equipped with both accepting and rejecting states. He studies some basic computability aspects of his machines and proves that a fuzzy language L is decidable if and only if L and L c are acceptable.
In Section 3, we make some essential modifications in previous definitions of fuzzy Turing machines in [1,7,13]. In these works, the notion of accepting or rejecting degree of an input is not computationally well defined and there are some cases which these degrees can not be computed. We modify these notions by (1) applying a BFS-based search in the computational tree of a given fuzzy Turing machine on an input and (2) obtaining an upper bound on the number of levels in our BFS-based search. If we reach an accepting or a rejecting configuration in a level, then the upper bound indicates the number of next levels which are needed to be traversed (at most) to determine the existence of another accepting or rejecting configuration. We also prove that computational power of fuzzy, generalized fuzzy Turing machines, and classical machines are the same.
In Section 4, we establish a new class of fuzzy Turing machines that we call Extended Fuzzy Turing Machines. These machines are equipped with indeterminacy states as a new type of states. There are silent transitions between indeterminacy states and all other non-accept and non-reject states with degree 1. By silent transition we mean a transition which does not change the position of head and the tape's content of the machine. Indeterminacy states are applied to catch some types of loops. The loop on an input is identified when its indeterminacy degree equals 0. Although extended fuzzy Turing machines are used to catch some types of loops, however they are not strong enough to solve the halting problem. We study some basic computability properties of these machines. Moreover, to each r.e. or co-r.e language, we correspond a fuzzy language which is indeterminable by an appropriate extended fuzzy Turing machine.

Preliminaries
In this section, we review some preliminaries from the literature.
where for each x ∈ X, µ A (x) represents the grade of membership of the element x ∈ X to A.

Definition 2.2. [4]
A t-norm is a binary operation * on [0, 1] satisfying commutativity, associativity, non-decreasing in both arguments, with the properties 0 * x = 0, 1 * x = x, for all x. Definition 2.3. [7] Let L 1 and L 2 be two fuzzy languages. Also, let * be a t-norm and * ′ be its dual t-conorm. The languages L 1 * L 2 and L 1 * ′ L 2 are defined as follows: In Turing machine T , Q is a set of states, Σ is a set of input symbols, Γ is a set of tape symbols, transition relation ∆ is a subset of Q × Γ × Q × Γ × {R, L}, q s ∈ Q is the starting state and F ⊆ Q is the set of final states. , η), where T = (Q, Σ, Γ, ∆, q s , F ) is a non-deterministic Turing machine such that Q is a set of states, Σ is a set of input symbols, Γ is a set of tape symbols, ∆ is a fuzzy subset of Q × Γ × Q × Γ × {R, L}, q s ∈ Q is the starting state, F ⊆ Q is a set of accepting and rejecting states, * is a t-norm, *

Corrections in Some Previous Definitions of FTMs
In this section we explain some computational objections about the definitions of fuzzy Turing machines in [1,7,13]. The Wiedermann's idea in defining a Fuzzy Turing Machine (FTM) is to establish an uncertainty degree for the acceptance of a given string or the membership degree of a string to the language of a FTM. In order to compute this degree, he applies a composition on a t-norm evaluation. He considers the class of FTMs as a fuzzy extension of non-deterministic Turing machines (NTMs), where each transition has a membership degree. In order to define the acceptance degree of an input w, he defines the degree d(C t ) of an arbitrary configuration C t , as the maximum degree of all computational path C 0 ⪯ C 1 ⪯ ... ⪯ C t leading from the initial configuration C 0 = q 0 w to C t : He also defines the accepting degree of a string w in a FTM F = (T , * , η) as follows: where, T = (Q, Σ, Γ, ∆, q s , F ) is a non-deterministic Turing machine and F contains only accept states. Note that q 0 w ⪯ * uq f v denotes that there is a computational path from q 0 w to uq f v. In [1,7], the same setting is followed. Bedregal and Figueira [1] change the maximum to supremum in definitions of d(C t ) and deg F (w). Moniri [7] assumes that F may contain both accept and reject states. Analogously, he defines accepting and rejecting degrees of a string w using a composition on a t-conorm * ′ instead of maximum, where * ′ is the dual of the t-norm * . For simplicity, we only continue discussion in the case of accepting degrees and our discussion can be extended to the case of rejecting degrees naturally. Here we use the notation deg F (w) for accepting degree of w.
The important point is that none of the definitions above for accepting degree of an input is not computationally well defined. There are some cases in which these degrees can not be computed. The machine T is non-deterministic and when the machine is running on an input w there are different computational branches. If there is a loop branch, then the machine continues searching in computational tree forever to find all accepting configurations. This prevents the process to halt. For example, the machine illustrated in Figure 1, has infinitely many computational paths leading from q s 01 to 0q I 1 and does not stop searching in the computational tree on input 01. Thus, degree of configuration 0q I 1 can not be computationally obtained using definitions in [1,7,13].
In order to resolve this problem, we make some changes in previous definitions. We take Moniri's definition of deg F (w) as the base and make/use the following changes/facts: • Any non-deterministic Turing machine can be simulated by a deterministic Turing machine using a BFS search method (See Theorem 3.16 [10]). By applying and using some modifications in the introduced BFS search method, we simulate the non-deterministic Turing machine T by a deterministic Turing machine. • We find an upper bound for the number of levels (in our BFS-based search) which are needed to be traversed (at most) to guarantee the existence of deg F (w), • We define path independent degree just for accepting and rejecting configurations.
gives values 0 and 1 in the following special cases: (i) a direct transition from initial state to an accept or a reject state is allowed to take degree 0, (2) transitions into or out from an indeterminacy state are allowed to take degree 1.
is a nondeterministic Turing machine. Let D be the deterministic Turing machine which simulates T using the BFS search method described in Theorem 3.16 [10]. We make the following modifications in the introduced BFS search method: (i) If in BFS search an accepting configuration in a node is reached, then we don't traverse any of the following configurations initiated from that node. This modification is important because for example if we have a silent move from accept state to itself, i.e. for some a ∈ Γ we have (q accept , a, q accept , a, S) ∈ ∆, then we have a branch consisting infinitely many accepting configurations which are not effective in computing deg F (w). (ii) Assume that we are at level n when D is simulating T on an input w. If we visit an accepting configuration in this level, then we can find an upper bound for the number of levels enough to be traversed (at most) to find another accepting configuration. If we could not find another accepting configuration in this level bound, then it is guaranteed that there is no other accepting configuration and we can halt and compute deg F (w).
In the sequel, we find the upper's level bound. Suppose that we reach an accepting configuration C with fuzzy degree d in level n. Consider our modified BFS in traversing the computational tree on input w. If we visit no new accepting configuration in next t = max d ′ ⌈log k d d ′ ⌉ levels, then we are sure that there is no other accepting configuration and we can compute deg F (w). In this upper bound, maximum is taken over all degrees d ′ of another configurations in level n and k is the maximum value in the range of η. Note that ∆ is a finite set and so the image of η over ∆ is a finite set and has a maximum value k.
Now we explain the idea to construct the given upper bound. In the definition of the degree of a given path, the t-norm * is applied between degrees of transitions along the path. So, if the degree of transitions is strictly less than 1, then the path's degree strictly decreases when we traverse from the current level to next one in the computational tree. Assume that the current configuration in level n is C ′ with degree d ′ and i is the minimum natural number which satisfies the condition d ′ * k i ≤ d. Thus, for configurations in an arbitrary path initiated from C ′ , the only ones effective in computing deg F (w) are (at most) in the next "i" levels. If we apply the same process for each arbitrary configuration in level n, then the upper bound t = max d ′ ⌈log k d d ′ ⌉ can be defined as the maximum of these numbers "i". Consequently, by the discussion above, we can modify previous definitions of FTMs or GFTMs to compute the accepting or rejecting degrees of inputs in general. In Section 4, we apply the modifications so far in defining the class of Extended Fuzzy Turing Machines.
Remark 3.1. Let E = (T , * , η) be an FTM and the computational tree of T on an input w has no loop branch, then the accepting degree of each arbitrary configuration can be computationally defined and also there is no restriction on degrees of transitions. A similar discussion also holds for GFTMs.
Remark 3.2. If in an FTM E = (T , * , η) (or similarly in a GFTM), the machine T be a deterministic Turing machine, then all accepting or rejecting degrees in [1,7,13] are computationally well-defined. Now we compare the computational power of fuzzy and generalized fuzzy Turing machines to classical Turing machines. Concerning Church-Turing thesis, this comparison is an important issue in computability theory. It is known that FTMs and GFTMs are extensions of NTMs. In the following it is shown that there is a classical Turing machine which can simulate each FTM (also, in the same setting for GFTMs) and so, FTMs, GFTMs and classical Turing machines have the same computational power. Proof. We give the sketch of the proof. Let F = (T , * , η) be an FTM, where, T = (Q, Σ, Γ, ∆, q s , F ) is a non-deterministic Turing machine. First note that the language of F includes pairs of the form (w, deg F (w)). Assume that D is the deterministic Turing machine which simulates T by applying our modified BFS search. We propose a classical 3-tapes deterministic Turing machine M which can simulate the machine F such that for each input w: Intuitively the machine M works as follows: • Its first tape is the input tape which remains unchanged during the computation, • Its second tape is the work tape and the simulation of T on w is executed on this tape, • Its third tape holds both the degree of the current configuration and also the degrees of visited accepting configurations.
If after visiting all accepting configurations (considering the upper level bound) the machine D halts on w, then the machine M computes the degree of w using the contents of its third tape and outputs this degree.
Likewise, it can be shown that there exists a classical Turing machine which can simulate each GFTM and give the accepting or rejecting degree of a given input.

Extended Fuzzy Turing Machines and Some Computability Results
In this section, we extend FTMs and GFTMs to machines with some new type of states that we call indeterminacy states. We also study some computability properties of these extended machines.

Extended Fuzzy Turing Machines.
Considering the modifications mentioned in Section 3, we define the class of Extended Fuzzy Turing Machines as an extension of the class of Generalized Fuzzy Turing Machines (GFTMs) defined by Moniri [7]. In defining extended fuzzy Turing machines, we consider some modifications with respect to Moniri's definition: (1) we specify indeterminacy states as a new type of states; (2) unlike Moniri, we assume that the transition relation ∆ is a crisp set and it is not fuzzy. where, q a and q r are accept and reject states, respectively.
Remark 4.2. Note that although the transition relation ∆ in a GFTM is considered as a fuzzy subset, but we define ∆ as a crisp subset in our model. In this way, η has the desired meaning.
Instantaneous description (ID) which is the unique description of a machine's tape, is defined as usual (see [12]). If α, α ′ are two IDs, then α ⪯ r α ′ means that there is a move in ∆ with truth degree r leading from α to α ′ in one step. Let α t be reachable from α 0 in t steps through the computational path α 0 ⪯ r1 α 1 ⪯ r2 ... ⪯ rt α t , then the truth degree of this path is defined as r 1 * ... * r t , if t ≥ 1 and 1, if t = 0. Due to non-determinism, an arbitrary configuration such as α t can be reached from α 0 through different computational paths, and so the degree of a configuration should be defined path independently. We define the path independent degree of a configuration just for accepting and rejecting configurations. Let α h be an accepting or a rejecting configuration, then we define its degree as follows: where, {a j ∶ j ∈ J} is the set of all truth degrees of all computational paths leading from the initial ID α 0 to α h such that the computational paths are traversed considering the obtained upper level bound for the given modified BFS search.
We call the computational path α 0 ⪯ r1 α 1 ⪯ r2 ... ⪯ rt α t , an accepting, a rejecting or an indeterminacy path, if α t is an accepting, a rejecting or an indeterminacy ID, i.e. in the form uq a v, or uq r v or uq i v, where u, v are two strings in Γ * and q a , q r , q i are states in A, R, I, respectively. If there exists at least one accepting (or rejecting) path on input w, then the accepting (or rejecting) degree of w, denoted by e E (w) (or e ′ E (w)), is defined as follows: where, α a (or α r ) is an accepting (or a rejecting) ID. Otherwise, if there is no desired path we define e E (w) = 0(= e ′ E (w)). Remark 4.4. In the definition of an EFTM, the set of states Q and the set of symbols Γ are finite and so the sets of accepting and rejecting configurations on a given input are finite. Therefore, unlike Moniri's, in Definition 4.3 we use maximum instead of supremum in the definitions of accepting and rejecting degrees.
Above, we computationally defined path independent truth degree for an accepting or a rejecting configuration only. Now, we define the path independent truth degree of an indeterminacy configuration α i such that it is not computationally definable in general. Assume that E = (T , * , * ′ , η) is an EFTM and w ∈ Σ * T . Let α i be an indeterminacy configuration on w, then we define: where, for each j, b j is the truth degree of a computational path leading from α 0 to α i . Note that d(α i ) = Inf j∈I a j is mathematically definable. If there exists at least one indeterminacy path on an input w, then we define the indeterminacy degree of w as follows: Consider the EFTM shown in Figure 1. It can be verified that the accepting degree of w 1 = 0 equals 0.14, and the indeterminacy degree of w 2 = 01 equals 0.
In the sequel, we give a proposition which characterizes the loops of classical Turing machines using EFTMs. We define the notion of an acceptable or a decidable fuzzy language similar to [7], but we define the new notion of an indeterminable fuzzy language in the following definition.  i.e. for each w ∈ L we have e E (w) < e ′′ E (w). Note that L c is defined as L c (x) = 1 − L(x). In the following we restate Proposition 3.5, Corollary 3.6 and Proposition 3.8 of [7], which also hold here. Corollary 4.10. If a fuzzy language L is decidable, then L c is decidable, Proposition 4.11. Let L 1 and L 2 be two fuzzy languages. Assume that L 1 and L 2 are accepted by EFTMs E 1 and E 2 equipped with the same t-norm * and let * ′ be the dual t-conorm of * . Then L 1 * ′ L 2 is accepted by an EFTM equipped with the same t-norm and t-conorm.

4.2.
Classical Languages and Extended Fuzzy Turing Machines. In this section, we propose a correspondence between the class of extended fuzzy Turing machines and the class of all classical r.e. or co-r.e. languages. Proposition 4.12. Let L be a classical r.e. language. There is an EFTM E such that for each arbitrary input w, we have w ∈ L if and only if E accepts w with a non-zero degree b, rejects it with degree 0 and indeterminates it with a degree strictly greater than b.
Proof. Suppose that M is a classical Turing machine which recognizes L. Without loss of generality we assume that M has only accepting states as its final states and eventually accepts an input w if w is in L and never halts on other inputs. Let t be a real number in (0,1). Construct the EFTM E as follows: • change the starting state of M to a non-starting state q, • correspond degree t to all transitions of M , • consider a starting state q s , • consider two nondeterministic transitions from q s to: (i) state q of the new modified version of M with degree t, (ii) a new rejecting state with degree 0 (note that here the degree 0 is allowed here by Definition 4.1), • consider transitions from all non-accept and non-reject states to an indeterminacy state with degree 1 and vice versa. It can be shown that for each input w if w ∈ L, then there is a number b ∈ (0, 1) such that e E (w) = b ≤ t and e ′′ E (w) > b, also e ′ E (w) = 0. Therefore, E is the desired machine.
Remind that in above proposition, the language L is indeterminable by E. Proposition 4.13. Let L be a (classical) co-r.e language. There is an EFTM E such that for each input w, we have w ∉ L if and only if E rejects w with a non-zero degree r, accepts it with degree 0 and indeterminates it with a degree strictly greater that r.
Proof. Since L is a co-r.e language then there exists a Turing machine M which recognizes L c . Let s be a real number in (0,1). Construct the EFTM E as follows: • replace the starting state of M with a new non-starting state q, • correspond degree s to each transition of M , • change the accepting states of M to rejecting states, • consider a starting state q s , • add two non-deterministic transitions from the start state q s : (1) to the state q with degree s, (2) to an accepting state with degree 0 (note that here the degree 0 is allowed here by Definition 4.1), • consider transitions from all non-accept and non-reject states to an indeterminacy state with degree 1 and vice versa. It can be shown that for each input w if w ∈ L c , then there is a number r ∈ (0, 1) such that e ′ E (w) = r ≤ s and e ′′ E (w) > r, also e E (w) = 0. Thus, E is the desired machine.

Final remarks
We gave some modifications in previous definitions of fuzzy Turing machines proposed in [1,7,13]. We applied a BFS-based search method, also obtained an upper level bound to define the notions of accepting and rejecting degrees of a given input, computationally. We introduced the class of Extended Fuzzy Turing Machines which are equipped with indeterminacy states. Finally, we used indeterminacy states to catch the loops of classical Turing machines.


extfieldsofstudy: Geography Medicine

text:
Urban Ecological Security Simulation and Prediction Using an Improved Cellular Automata (CA) Approach—A Case Study for the City of Wuhan in China

Ecological security is an important research topic, especially urban ecological security. As highly populated eco-systems, cities always have more fragile ecological environments. However, most of the research on urban ecological security in literature has focused on evaluating current or past status of the ecological environment. Very little literature has carried out simulation or prediction of future ecological security. In addition, there is even less literature exploring the urban ecological environment at a fine scale. To fill-in the literature gap, in this study we simulated and predicted urban ecological security at a fine scale (district level) using an improved Cellular Automata (CA) approach. First we used the pressure-state-response (PSR) method based on grid-scale data to evaluate urban ecological security. Then, based on the evaluation results, we imported the geographically weighted regression (GWR) concept into the CA model to simulate and predict urban ecological security. We applied the improved CA approach in a case study—simulating and predicting urban ecological security for the city of Wuhan in Central China. By comparing the simulated ecological security values from 2010 using the improved CA model to the actual ecological security values of 2010, we got a relatively high value of the kappa coefficient, which indicates that this CA model can simulate or predict well future development of ecological security in Wuhan. Based on the prediction results for 2020, we made some policy recommendations for each district in Wuhan.

Introduction
Ecological security refers to a non-threatened state [1][2][3], which may include the state of human life, health, well-being, basic rights, livelihood resources, necessary resources, social order, and human ability to adapt to environmental changes. Ecological security reflects the degree of ecosystem integrity and health [1], and it may provide guarantee for sustainability of an eco-economic system [4]. Ecological security research began in the last century [5][6][7]. Important literature on this concept includes the "Ecological risk assessment outline" published by the United States Environmental Protection [8] and the "Global vulnerability assessment framework for environmental risk" proposed by Clark et al. [9]. Ecological security is a concept opposite to ecological risk. Ecological security measures the safeguard of ecological systems while ecological risk assesses the danger or destruction of ecological systems. However ecological security and ecological risk are related. In fact, parts of the concept of ecological security come from ecological risk [10][11][12]. Because many ecological disasters or degradation are caused by human activities, measurement of ecological security should represent the capability of human activities. Ecological security may be assessed using an indicator system for environmental pressures or degradation.
Urban ecological security is an important research topic [13]. The urban area, as the center of socio-economic development, is a highly populated ecological system and is more fragile in terms of ecological security [2,4,14]. It is, therefore, imperative to study urban ecological security. In addition, studying urban ecological security may benefit for carrying out efficient environmental management [15]. With increasing concerns over urban ecological security, many methods have been proposed in the literature for assessing urban ecological security. These methods mainly include mathematical models, system dynamics methods, ecological footprint methods [16], and landscape ecology methods [15,17,18]. The commonly used mathematical model methods include the pressure-state-response (PSR) model, the comprehensive index method [19], the analytic hierarchy process method [20,21], the gray correlation degree method [22], the matter element evaluation method, and the principal component projection method. However, some of these methods, such as the comprehensive index methods and ecological footprint methods have limitations in the study of urban ecological security. These methods may be good at evaluating the status of an administrative region, but they cannot reflect the dynamic process of the ecological security change, especially at a fine scale. Among these mathematical methods, the PSR model for human and environmental systems [23][24][25][26], which was jointly developed by the United Nations Economic Cooperation and Development (OECD) and the United Nations Environment Program (UNEP), and the comprehensive evaluation index system established by U.S. Environmental Protection Agency (USEPA) [27] are particularly important.
However, most of the research on urban ecological security using the aforementioned methods focused on evaluating current or past status of the ecological environment. Very few studies in the literature have been carried out to simulate or predict future changes in urban ecological security. Nevertheless, predicting future changes in urban ecological security is important [28,29]. The predicted change of urban ecological security could be used to make reasonable environmental policies for urban sustainability development [30][31][32]. The predicted change of urban ecological security could also be useful for urban planners to develop more realistic urban development plans. However, there is a lack of literature in simulating or predicting future changes of urban ecological security. Based on our best knowledge, we only found a couple of publications in literature for forecasting urban ecological security. For example, Gong et al. [28] and Yang et al. [29] forecasted urban ecological security using the classical cellular automata (CA) models. However, the classical CA models consider only short-range spatial dependence and do not consider long-range spatial dependence (e.g., not immediately connected neighbor cells but located far away cells) [33]. Therefore, it is difficult to make an accurate simulation or prediction of the dynamic urban ecological security changes among the neighborhoods.
In addition, there is little literature exploring the urban ecological environment at a fine scale. In most existing literature, urban ecological security was studied at coarse scales (global, regional, city, or county levels) [34,35]. In those studies, a whole city or a county was given only one security value [21,36]. It is, therefore, impossible to figure out how things are going on inside the city [33]. Nevertheless, urban areas are heterogeneous [37]. Urban ecological security may be different at different locations in the fine scales such as district level. It is, therefore, improper to treat the city as homogeneous and give the whole city one urban ecological security value. One security value cannot represent the uneven status of ecological security inside a city. Different districts or neighborhoods in a city may have different ecological security values. Therefore, we need to simulate urban ecological security values at a fine scale to give us more information about the heterogeneity of a city. This may be used for better ecological understanding of integrated human-natural urban systems and undertake wise urban planning. Compared to the coarse scales such as city scale, a fine scale could help to find out the uneven status of ecological security inside a city and could give more information about the heterogeneous ecological security inside a city.
To overcome these limitations and fill-in the literature gaps, in this study, we simulated and predicted urban ecological security at a fine scale (district level) using an improved CA approach to better understand spatial heterogeneity of urban ecosystems. There are two major new contributions of the improved CA approach. First, we used the pressure-state-response (PSR) method based on grid-scale data to evaluate urban ecological security. We considered urbanization as the most significant driving force causing ecological and environmental problems. We carefully quantified the ecological status of the urbanization areas, because the biggest impetus to the ecological security problem is the increased urbanization processes. To better quantify urban ecological security, we used three dimensions to characterize urban growth. The three dimensions include urban growth degree (UGD), urban growth capacity (UGC), and urban growth pattern (UGP). We evaluated ecological security states from the view of environmental structure, ecological function value, and resilience of the ecosystem.
Second, in simulating and predicting the urban ecological security process, we extended the classical CA model by introducing the geographic weighting regression (GWR) concept into the transformation rules of the classical CA model. GWR can take into account spatial variation of the sampling unit [38], therefore, it may make simulation of the dynamic urban ecological security changes among different neighborhoods more accurate.
We applied the improved CA approach in a case study-simulating and predicting urban ecological security for the city of Wuhan of Central China. Through the case study, we aimed to verify the potential of the improved CA method for more accurately simulating and predicting urban ecological security.

Study Area
The city of Wuhan is one of the fastest growing cities in central China, located in the middle reaches of the Yangtze River (see Figure 1). It is the capital of Hubei Province. Wuhan city constitutes seven districts and has a total area of 8549 km 2 , of which 39.27% consists of plains and 18.17% is hilly and mountainous regions. Wuhan has experienced rapid urbanization and a significant economic growth since last century. From 2005 to 2010, the urban area expanded to 552.61 km 2 (Wuhan Municipal Bureau of Statistics 2005, 2010). Its population grew from 8.58 million to 10.61 million and its economy expanded from 226.12 billion to 1090.56 billion yuan. The city has been witnessed an obvious and rapid urbanization.

Data Sources
Land use data were mainly derived from remotely sensed images. The remotely sensed images include two multi-spectral Landsat MSS/TM/ETM+ images in the years 2005 and 2010. The resolution of the Landsat MSS/TM/ETM+ images is 30 m. The images were downloaded from the U.S. Geological Survey (USGS). To ensure accuracy of the image classification results, all the selected images are cloudless and with thick vegetation coverage. Before classification, all the images and other Geographic Information System (GIS) data were geometrically rectified to the World Geodetic System-1984 Coordinate System (WGS-84) and the Universal Transverse Mercator coordinate system. In total, 90 ground control points were selected as supplementary data. We did the classification using the maximum likelihood method. According to natural ecological and land use attributes, land uses/covers from remotely sensed images were always classified into six types. However, grassland has a very small proportion in Wuhan. The images were, therefore, classified into other five types of land use/cover, namely construction land, forest, farmland, water, and bare land. Meanwhile, we randomly selected 200 training samples with sufficient information on the spectral characteristics for each land use type. In addition, we selected 100 field points for checking classification accuracy. The overall classification accuracy values for the images are from 81.41% to 88.00%, and the kappa coefficient values are from 0.77 to 0.85, which indicate a high classification accuracy [39,40]. Based on the classification results, the proportions of construction land, forest, farmland, water, and bare land in Wuhan in 2010 are 7.78%, 8.34%, 61.99%, 21.40%, and 0.49%, respectively. Population data were collected from the fifth census and the sixth census in China at the street level. Socio-economic data mainly include national and regional planning documents and policies, which were downloaded from the government website [41].

Data Sources
Land use data were mainly derived from remotely sensed images. The remotely sensed images include two multi-spectral Landsat MSS/TM/ETM+ images in the years 2005 and 2010. The resolution of the Landsat MSS/TM/ETM+ images is 30 m. The images were downloaded from the U.S. Geological Survey (USGS). To ensure accuracy of the image classification results, all the selected images are cloudless and with thick vegetation coverage. Before classification, all the images and other Geographic Information System (GIS) data were geometrically rectified to the World Geodetic System-1984 Coordinate System (WGS-84) and the Universal Transverse Mercator coordinate system. In total, 90 ground control points were selected as supplementary data. We did the classification using the maximum likelihood method. According to natural ecological and land use attributes, land uses/covers from remotely sensed images were always classified into six types. However, grassland has a very small proportion in Wuhan. The images were, therefore, classified into other five types of land use/cover, namely construction land, forest, farmland, water, and bare land. Meanwhile, we randomly selected 200 training samples with sufficient information on the spectral characteristics for each land use type. In addition, we selected 100 field points for checking classification accuracy. The overall classification accuracy values for the images are from 81.41% to 88.00%, and the kappa coefficient values are from 0.77 to 0.85, which indicate a high classification accuracy [39,40]. Based on the classification results, the proportions of construction land, forest, farmland, water, and bare land in Wuhan in 2010 are 7.78%, 8.34%, 61.99%, 21.40%, and 0.49%, respectively. Population data were collected from the fifth census and the sixth census in China at the street level. Socio-economic data mainly include national and regional planning documents and policies, which were downloaded from the government website [41].

The Flow Chart of the Methods
The methods include two main components: evaluation of urban ecological security using the PSR model and simulation and prediction of ecological security using the CA model with the GWR concept. Figure 2 illustrates the flow chart of the methods. The procedure of the methods is introduced as follows: Step 1: We evaluated the ecological security classes of 2005 and 2010 using the PSR model.
Step 2: We used the data of 2005 to conduct a five-year simulation based on the CA rules. We obtained the simulated result of 2010 based on the five-year simulation.
Step 3: We compared the simulated result of 2010 and the evaluated result of 2010 to verify the simulation accuracy in the second step.
Step 4: We repeated step 2 and step 3 until a satisfactory level of accuracy was obtained in step 3. Then we recorded the CA rules of the simulation. The total number of the iterations of the 5-year simulation was recorded as N.
Step 5: We then obtained the predicted result for 2020, which is a 15-year simulation, by using the rules obtained in step 4. We took 3N as the total number of iterations in this 15-year simulation.

The Flow Chart of the Methods
The methods include two main components: evaluation of urban ecological security using the PSR model and simulation and prediction of ecological security using the CA model with the GWR concept. Figure 2 illustrates the flow chart of the methods. The procedure of the methods is introduced as follows: Step 1: We evaluated the ecological security classes of 2005 and 2010 using the PSR model.
Step 2: We used the data of 2005 to conduct a five-year simulation based on the CA rules. We obtained the simulated result of 2010 based on the five-year simulation.
Step 3: We compared the simulated result of 2010 and the evaluated result of 2010 to verify the simulation accuracy in the second step.
Step 4: We repeated step 2 and step 3 until a satisfactory level of accuracy was obtained in step 3. Then we recorded the CA rules of the simulation. The total number of the iterations of the 5-year simulation was recorded as N.
Step 5: We then obtained the predicted result for 2020, which is a 15-year simulation, by using the rules obtained in step 4. We took 3N as the total number of iterations in this 15-year simulation.

PSR Evaluation Model
PSR model is a comprehensive evaluation model for environmental problems proposed by the United Nations Economic Cooperation and Development (OECD). It is composed of the pressure index for ecological environment problems, state index for eco-environment status, and response index for social and economic measures that quantify human's response to ecological problems. The pressure index consists of three indicators: urban growth degree (UGD) [29], urban growth capacity (UGC) [42][43][44][45][46][47][48], and urban growth pattern (UGP) [49][50][51]. The ecological security status of an ecosystem was considered from three aspects: ecological function, ecological structure, and ecological resilience [52]. In order to make the research scale more refined and better reflect ecological state of the administrative units, the indicators of these three aspects were chosen at landscape scale [53]. The landscape ecological function (LEF) was assessed mainly based on the ecological service values. The most frequently used method proposed by Costanza et al. [54] and Xie et al. [55] was used to measure

PSR Evaluation Model
PSR model is a comprehensive evaluation model for environmental problems proposed by the United Nations Economic Cooperation and Development (OECD). It is composed of the pressure index for ecological environment problems, state index for eco-environment status, and response index for social and economic measures that quantify human's response to ecological problems. The pressure index consists of three indicators: urban growth degree (UGD) [29], urban growth capacity (UGC) [42][43][44][45][46][47][48], and urban growth pattern (UGP) [49][50][51]. The ecological security status of an ecosystem was considered from three aspects: ecological function, ecological structure, and ecological resilience [52]. In order to make the research scale more refined and better reflect ecological state of the administrative units, the indicators of these three aspects were chosen at landscape scale [53]. The landscape ecological function (LEF) was assessed mainly based on the ecological service values. The most frequently used method proposed by Costanza et al. [54] and Xie et al. [55] was used to measure the ecological service values for different land use types in this study. The landscape ecological structure (LES) was measured using (1) Shannon's diversity index (SHDI), which measures the number of different land cover classes present; (2) the area-weighted plots fractal dimension index that reflects shape complexity across a range of spatial scales; and (3) modified Simpson's evenness index, which ensures an even distribution of areas among patch types [51].
Landscape ecological resilience (LER) is an important factor for sustainability of regional ecology. LER measures the disturbance that an ecosystem can bear until its characteristics change [56]. In this study, we used the fragility index to measure this disturbance. It is calculated by the inverse of number of patches per unit area (PD). The 8-way neighbors method was used to define patches that can well represent the ecological resilience of a region [52].
The ecological security responses were interpreted by policy measures such as ecological zoning and environmental plans. From 2005 to 2010, the specific ecological responses of Wuhan City included "ecological bottom line", "ecological development area", "prospective ecological area", "tricyclic rectification area", and "urban concentrated construction area".
All of the aforementioned landscape indices were calculated using the software Fragstats [57,58]. Delphi method was used to assign the weight for each index. Compared to the objective weight allocation methods, Delphi method can better reflect the regional differences, and is more in line with the objective of this study so this method was employed. Considering the mechanism of the PSR model and the richness of the data, we assigned the state a weight of 0.5, which can better display the ecological security of Wuhan city. Since both the UGD and UGC in the pressure part and LEF and LES in the state part played important roles in the evaluation of urban ecological security, they were assigned higher weights. Table 1 lists the specific PSR evaluation indices and their weights.
Before performing the evaluation, the original data of different dimensions were standardized. The extremum standardization method was used for this purpose.

GWR-Based Cellular Automata (CA) Model
Cellular automata (CA) was a commonly used model for simulating and predicting urban growth [42,58,59]. CA models can take into account dynamic changes of a city caused by complex factors such as nature and socio-economy, and use simple rules to simulate the complicated temporal and spatial dynamic processes [60][61][62]. In addition, it has the advantage of easily integrating itself with high-resolution remotely sensed images and geographic information systems [63]. This makes it popular for simulating urban dynamic changes. Like general urban dynamic change processes, urban ecological security could be described as a synthesis of components of coupled human and natural ecosystems. Urban ecological security exhibits complex spatial and temporal processes. Therefore, urban ecological security could be simulated using CA models.
The CA simulation model used in this paper is based on the concept of GWR. The cellular automata is composed of five basic elements: cell, cell space, neighbor, transformation rules, and time sets. The neighborhood function is the core component of CA's conversion rules, which embody the essence of CA's "bottom-up" self-organization evolution. The classical CA algorithms usually employ the traditional global regression statistics such as a logistics regression model to calculate the conversion rules. There is, however, a limitation in using the global regression: the estimated parameters using the global logistics regression model are just mean values of the explanatory variables. The global parameters cannot reflect spatial heterogeneity of the relationships between the dependent variable and the explanatory spatial variables. Based on the first law of geography, near objects are more related than distant objects. Therefore, the explanatory variables may have different degrees of influence on conversion rules for different geographical locations. The conversion rules may not be a trend distribution with a single spatial structure. They may have spatial instability.
To reflect spatial instability in conversion among different land use types, we employed a local regression model-GWR-instead of the global logistics regression model to calculate the conversion rules. GWR is an extension of the traditional multiple linear regression toward a local regression. The regression coefficients of GWR are specific to a location rather than being global estimates [64,65]. The local estimation of the regression coefficients of GWR is expressed by the following equation: where (m i , n i ) is the spatial location of the ith observation and β j (m i , n i ) is the value of the jth parameter at point location i. The regression parameters of this equation are estimated at each location i(m i , n i ). GWR was specifically designed to deal with the spatial non-stationarity of regression coefficients between the dependent variable and explanatory variables by measuring those coefficients locally using local data [66][67][68]. GWR has capability of incorporating various auxiliary variables with spatially varied correlation coefficients. The GWR concept, therefore, may more accurately reflect the characteristics of spatial heterogeneity of ecological dynamic processes. Hence, in the processes of deciding conversion rules, we used GWR instead of global logistics regression to obtain cell conversion probability. We pplied the GWR in the CA rules mainly for calculating the probability from one assigned class to other classes. We first used the GWR models to find out relationships between a cell state value and each explanatory variable over space. Based on the relationships we then calculated the conversion probability of each cell to one specified type at each spatial position using the following formula: where p i represents the transition probability to one cell state value at position i, and y i is the regression coefficient equation for position i that was obtained by GWR.
In order to combine with the CA model, it is necessary to give one specific probability for each spatial location. We used the concept of roulette wheel to decide which class would be selected substantially based on the aforementioned calculated probability values. At the same time, we also considered whether or not the selected class values would satisfy the pre-set cell transition probability values.
We forecasted the ecological security status and patterns of 2020 based on the ecological security levels of 2005 and 2010 using the GWR tool in ArcGIS 10.2 software (Environmental Systems Research Institute, Inc., Redlands, CA, USA) [69]. We took the ecological security level as the dependent variable and the indicators calculated in the aforementioned PSR evaluation as the independent variables to obtain the regression coefficients of each variable. Taking simulation of the urban ecological security of 2005-2010 for example, to simulate ecological security status in 2010, we used the number of ecological evaluation classes in 2010 as the termination condition based on the ecological security patterns of 2005. In the process of cell conversion iterations, we first estimated the probability of the currently iterated cell that would turn to each ecological security class. We then determined a class by the method of roulette wheel and the judgement about whether or not the conversion probability of the class satisfied the pre-set cell transition probability. If it satisfied, then the cell would be turned to that class. We stopped the iteration after all of the ecological security classes were satisfied. The total number of iterations was recorded as N. For forecasting the urban ecological security of 2020, we took 3N as the total number of iterations from the ecological security patterns of 2005.

PSR Evaluation
The PSR evaluation results of ecological security values were taken into account of ecological pressure, ecological state, and ecological response. The evaluated ecological security values are distributed in the range of 0-0.7. For a better view, we divided the range into five levels, from Class 1 to Class 5. If a value is less than 0.15, the grid cell will be classified into Class 1. Class 1 means a location with the lowest security. If a value belongs to the range from 0.15 to 0.3, the location will be classified into Class 2. Similarly, Class 3 has values ranging from 0.3 to 0.4. A cell location with Class 3 indicates that it is located at the area with the mid-security level. Class 4 has values ranging from 0.4 to 0.5. The highest level of the urban ecological security is Class 5, which has values greater than 0.5.       To find out the reason for ecological security changes from 2005 to 2010 in Wuhan, we analyzed the urban ecological security from the perspectives of ecological pressure (P), ecological state (S), and ecological response (R), respectively. Figure 6 illustrates the P/S/R changes from 2005 to 2010. As Most areas in the 13 administrative districts also retain unchanged ecological conditions. For example, more than 80% of the areas in Wuchang and Jianghan Districts maintain their original ecological security levels. Except for Hannan District, the other 12 districts have enhanced conditions as the second largest proportion. In total 16.76% of the areas in the Hannan District have weakened ecological conditions while only 11.41% of areas in that district have enhanced conditions. Hannan District is the only district that has a larger proportion of areas with weakened conditions than areas with enhanced conditions. There are relatively large ecological security-enhanced areas in Jiangxia, Huangpi, Xinzhou, Caidian, Hongshan and Dongxihu Districts. The proportions of enhanced areas in Hongshan and Qingshan Districts are the two highest proportions in all of the districts, at 39.27% and 38.21%, even if their enhanced areas are not as large as those in the previously mentioned districts because of their small total district areas. Except for Hongshan and Qingshan Districts, Hanyang and Dongxihu Districts have enhanced proportions larger than 30%.
To find out the reason for ecological security changes from 2005 to 2010 in Wuhan, we analyzed the urban ecological security from the perspectives of ecological pressure (P), ecological state (S), and ecological response (R), respectively. Figure 6 illustrates the P/S/R changes from 2005 to 2010. As shown in Figure 6, the overall ecological environment pressure of Wuhan has significantly increased from 2005 to 2010, and only small areas around the main city have decreased ecological pressure. The ecological pressure in north Huangpi District and northeast Xinzhou District has little change. The areas of ecological improvement are mainly located in the southern part of Huangpi District and some areas in Qingshan, Hongshan, and Dongxihu Districts. The northeastern part of Wuhan city that belongs to the Xinzhou area and the southwest part that belongs to the Hannan District show a weakened ecological environment.
The ecological state basically remained unchanged from 2005 to 2010, and the ecological response was not significant either. There is only a slight decline of the ecological security in the central city. Nevertheless, the construction expansion from the city center to the surrounding areas of the city was significant, which caused most of the surrounding urban areas to have increased ecological environment pressure. The ecological responses in these areas were the most significant. However, the ecological security state of these areas remained basically unchanged or was even slightly enhanced. This is the main ecological security improvement in Wuhan. The ecological response in the suburb areas was small. If combined with urbanization, the ecological security situation remained basically unchanged or slightly decreased from 2005 to 2010.
We used expert knowledge and a literature review to find out the important drivers of ecological security changes in Wuhan. We found that the better ecological planning management and policies greatly accounted for ecological security changes of Wuhan City from 2005 to 2010. The changes in urban planning and land use policies were linked to modifications of landscapes and ecosystems in Wuhan. Good planning, management, and policies played an important role in driving the security changes. These planning management and policies were mainly applied to the outer ring areas around the main urban area, which include Huangpi, Dongxihu, Caidian, Jiangxia, Xinzhou, Hongshan Districts, and the eastern part of Qingshan District. Based on the expert knowledge and literature review we also found that the ecological security environment changes during 2005-2010 resulted from comprehensive multi-factors such as social, political and economic factors. For example, population growth and changing economic conditions were two causes of ecological security change. Urbanization and extended peri-urban settlement threatened various ecosystem processes and drove ecological security change from 2005 to 2010. Rapid economic development and construction efforts in Wuhan led to major ecological security changes from 2005 to 2010. In general, Wuhan has undergone rapid development period in recent years, especially in the central city, where development and construction efforts continued to increase in its main areas. This resulted in a widespread ecological pressure. However, even with the wide and fast development of the city, there was no sharp decline of the ecological state.

The Simulation by the Proposed CA Model
We used the proposed CA model with the GWR concept to simulate the ecological security situation in 2010 based on the ecological security situation in 2005. We then compared and verified the simulation results with the actual ecological security values in 2010. Kappa coefficient statistics were used to measure the simulation results. The calculated kappa coefficient is 0.756, which indicates that the simulation results are good and the proposed CA model can be used to predict future ecological security.
The evaluated ecological security results in 2010 and the simulated ecological security results in 2010 are shown in Table 3 and Figure 7. From Figure 7, it can be seen that the proportion of each class in the simulated results in 2010 is roughly the same as the actually evaluated results in 2010. From Table 2 it can be seen that the areas with a fair ecological security in both the evaluated results and the simulated results still occupy the largest proportion. However, Class 4 occupies 12.80% in the simulated results and 21.85% in the evaluated results. Class 4 has the largest difference among the five classes. Only the simulated areas of this class are smaller than the actual areas, and the simulated areas of the other four classes are larger than the actual areas.

The Simulation by the Proposed CA Model
We used the proposed CA model with the GWR concept to simulate the ecological security situation in 2010 based on the ecological security situation in 2005. We then compared and verified the simulation results with the actual ecological security values in 2010. Kappa coefficient statistics were used to measure the simulation results. The calculated kappa coefficient is 0.756, which indicates that the simulation results are good and the proposed CA model can be used to predict future ecological security.
The evaluated ecological security results in 2010 and the simulated ecological security results in 2010 are shown in Table 3 and Figure 7. From Figure 7, it can be seen that the proportion of each class in the simulated results in 2010 is roughly the same as the actually evaluated results in 2010. From Table 2 it can be seen that the areas with a fair ecological security in both the evaluated results and the simulated results still occupy the largest proportion. However, Class 4 occupies 12.80% in the simulated results and 21.85% in the evaluated results. Class 4 has the largest difference among the five classes. Only the simulated areas of this class are smaller than the actual areas, and the simulated areas of the other four classes are larger than the actual areas.

Prediction Results of 2020
We used the proposed CA model to forecast the ecological security situation of Wuhan in 2020. Figure 8 illustrates the predicted results of the ecological security in 2020. As shown in the figure, areas with comparative insecurity (Class 1 and Class 2) are projected to be more concentrated by 2020, especially in Jianghan, Jiangan and Qiaokou Districts. Both of Class 1 and Class 2 occupy more than 70% of the administrative areas in these three districts, and they are the major classes in these three districts. Hanyang District and Wuchang District belong to the areas with a low ecological security, because the proportions with ecological unsafe areas in these two districts are larger than those with secure areas. Caidian, Huangpi, Xinzhou, and Jiangxia Districts belong to the districts with good ecological conditions. These four districts have large areas with mid-security level. The total area with security level higher than the middle level is larger than that with security level lower than the middle level in these four districts. Hannan and Hongshan Districts are the districts with the best ecological conditions in Wuhan because of their large proportions of Class 4 and Class 5. Qingshan District has a big internal difference: the ecological safe areas and unsafe areas account for 42% and 44%, and both are widely distributed. On the contrary, the Dongxihu District has the smallest internal heterogeneity.

Prediction Results of 2020
We used the proposed CA model to forecast the ecological security situation of Wuhan in 2020. Figure 8 illustrates the predicted results of the ecological security in 2020. As shown in the figure, areas with comparative insecurity (Class 1 and Class 2) are projected to be more concentrated by 2020, especially in Jianghan, Jiangan and Qiaokou Districts. Both of Class 1 and Class 2 occupy more than 70% of the administrative areas in these three districts, and they are the major classes in these three districts. Hanyang District and Wuchang District belong to the areas with a low ecological security, because the proportions with ecological unsafe areas in these two districts are larger than those with secure areas. Caidian, Huangpi, Xinzhou, and Jiangxia Districts belong to the districts with good ecological conditions. These four districts have large areas with mid-security level. The total area with security level higher than the middle level is larger than that with security level lower than the middle level in these four districts. Hannan and Hongshan Districts are the districts with the best ecological conditions in Wuhan because of their large proportions of Class 4 and Class 5. Qingshan District has a big internal difference: the ecological safe areas and unsafe areas account for 42% and 44%, and both are widely distributed. On the contrary, the Dongxihu District has the smallest internal heterogeneity.

Conclusions
In this paper, we proposed an improved CA approach for simulate and predict urban ecological security at a fine scale (district level). We proposed to first evaluate urban ecological security using the PSR method, and then simulate and predict urban ecological security using the CA model with the GWR concept. We applied the proposed CA approach in a case study-simulating and predicting urban ecological security in Wuhan, China. We first evaluated the ecological security levels in 2005 and 2010 in Wuhan city using the PSR method. We then simulated and predicted the future urban ecological security situations in 2010 and 2020 using the improved CA model with the GWR concept. The results show that our improved CA model is effective and feasible for simulation or prediction of future ecological security.
Although the simulated results are good, there were some limitations in our case study. For example, the Landsat images were used in this case study. The Landsat imagery has a low resolution and many mixed pixels. The mixed pixels made it difficult to accurately classify the satellite images. As the input data for our models, the inaccurate classification results from the satellite images affected the accuracy of our final simulation results. For example, Class 4 generated the largest error in the simulation results of 2010-there is about 9% difference between the simulation results and the evaluated results of Class 4. This has caused the overall simulated ecological security situation lower than the overall actual situation. The other limitation is validation of the CA model. It is impossible to obtain a perfect simulation result using any model. Therefore, it is important to validate the improved CA model for practical applications. In this study, we only used the kappa coefficient indices for validation purpose. However, the Kappa indices have limitations, and sometimes they may be misleading for purposes of accuracy assessment and map comparison [70]. Therefore, we may consider adopting more measures, such as quantity disagreement and allocation disagreement from a cross-tabulation matrix, for validation purposes in the future.

Conclusions
In this paper, we proposed an improved CA approach for simulate and predict urban ecological security at a fine scale (district level). We proposed to first evaluate urban ecological security using the PSR method, and then simulate and predict urban ecological security using the CA model with the GWR concept. We applied the proposed CA approach in a case study-simulating and predicting urban ecological security in Wuhan, China. We first evaluated the ecological security levels in 2005 and 2010 in Wuhan city using the PSR method. We then simulated and predicted the future urban ecological security situations in 2010 and 2020 using the improved CA model with the GWR concept. The results show that our improved CA model is effective and feasible for simulation or prediction of future ecological security.
Although the simulated results are good, there were some limitations in our case study. For example, the Landsat images were used in this case study. The Landsat imagery has a low resolution and many mixed pixels. The mixed pixels made it difficult to accurately classify the satellite images. As the input data for our models, the inaccurate classification results from the satellite images affected the accuracy of our final simulation results. For example, Class 4 generated the largest error in the simulation results of 2010-there is about 9% difference between the simulation results and the evaluated results of Class 4. This has caused the overall simulated ecological security situation lower than the overall actual situation. The other limitation is validation of the CA model. It is impossible to obtain a perfect simulation result using any model. Therefore, it is important to validate the improved CA model for practical applications. In this study, we only used the kappa coefficient indices for validation purpose. However, the Kappa indices have limitations, and sometimes they may be misleading for purposes of accuracy assessment and map comparison [70]. Therefore, we may consider adopting more measures, such as quantity disagreement and allocation disagreement from a cross-tabulation matrix, for validation purposes in the future.
We discussed the reasons for the ecological security changes in Wuhan from 2005 to 2010. The ecological security environment changes between 2005 and 2010 resulted from multiple factors. The ecological security changes were mainly caused by the rapid economic development and construction efforts in Wuhan. Better ecological planning management practices and policies in 2010 also led to the ecological security improvement. These planning management and policies were mainly applied to the outer ring areas around the main urban areas.
Based on the prediction results of 2020, we would like to make the following policy recommendations for each district in Wuhan. For the districts that are continuing to have a comparably high ecological security and have an obvious improvement from 2005 to 2020, such as Hannan, Hongshan, and Jiangxia Districts, the government only needs to maintain the ecological environment. For the districts that will continue to be the concentration areas with the worst ecological security situations of Wuhan city in 2020, such as the Jianghan, Jiangan, and Qiaokou Districts, the government should pay more attention to a centralized development for them. The government should undertake more contiguous land use planning instead of fragmented land use planning in order to make full use of each piece of land. For some districts such as Qingshan District, diverse plans should be made inside the district considering its high internal spatial heterogeneity in 2020.


Experiments of two pupil lateral motion tracking algorithms using a Shack–Hartmann sensor

Abstract Pupil stability is one of the factors which limit the performance and operational stability of adaptive optics (AO) systems. This paper analyses two pupil-tracking methods to measure the lateral pupil shift: the first one utilizes the fluxes in all outer edge sub-apertures of a Shack–Hartmann sensor and the second one utilizes the real-time interaction matrix used in an AO system. Experiments with 9 × 9 Shack–Hartmann sensor are conducted to verify both pupil-tracking algorithms. The results show that both algorithms are effective, after two correction steps, the residual pupil shift is reduced to less than 5% of a Shack–Hartmann sub-aperture.

Introduction
Nowadays, adaptive optics (AO) has become one of the key technologies for large ground-based telescopes (1) and has been adopted in almost all 8-10 m class telescopes to improve their optical quality while observing through the turbulent atmosphere. To make full use of the AO system, the wavefront sensing channels and correcting device channels need to be registered precisely during closed-loop operation to guarantee optimal performance of the AO system (2).
GRAVITY is an AO-assisted beam combiner for the second generation very large telescope (VLT) interferometer (VLTI) instrumentation (3,4) within which there is a system named Coudé infrared adaptive optics (CIAO) (5,6) used to ensure the required wavefront quality for the interferometer. Figure 1 shows the schematic of a CIAO system, in total four CIAO units are available for GRAVITY to correct the wavefronts of the beams of the four telescopes of VLT. For each telescope, after being corrected by the deformable mirror (MACAO DM), the light is split into two parts by the Star Separator: the light of the wavefront reference (star) goes into the near-infrared (IR) wavefront Sensor, which composes the AO system together with the MACAO DM. The light of the object and the reference star goes into the Beam Combiner Instrument, where they interfere with the corresponding signals from other telescopes.
Misalignments or the motions of the optical elements between the deformable mirror and the Shack-Hartmann wavefront sensor (SHWFS) can induce both lateral and rotational pupil shifts during the telescope's operation. These shifts lead to the mis-registration between the optical pupils defined by the DM and the wavefront sensor subapertures which are registered to the actuators of the DM. As a result, pupil shift can lead to the loss of light as well as the degradation of the performance of the AO system (7). This paper focuses on the lateral pupil shift in the CIAO system which exists between the DM and SHWFS. A pupil actuator (field lens) located in CIAO cryostat can be used to stabilize the pupil position (8) (Figure 2). This requires a precise measurement of the actual pupil position. This paper analyses two pupil-tracking methods: Section 2 describes the first method which utilizes the fluxes in all outer edge sub-apertures of the Shack-Hartmann sensor and Section 3 describes another method which utilizes the real-time interaction matrix obtained during the closed-loop control of the AO system. Experiments are conducted to verify these two methods; the results are shown, respectively, in Sections 2 and 3. In Section 4, we make a comparison of these two algorithms. Section 5 summarizes our conclusions.
where θ i is the angle between the centre of the sub-aperture O i and X-axis ( Figure 3).

Result of the experiment using the intensities method
We use the CIAO SHWFS in the laboratory to perform the experiment. As shown in Figure 3, the fluxes in the In the intensities method, we choose the fluxes in the 24 outer edge sub-apertures to measure the pupil shift. The integrated fluxes of the sub-apertures are marked with O 1 , O 2 , … , O 24 (outer ring sub-apertures) and I 1 , I 2 , I 3 , I 4 (inner ring sub-apertures). When the lateral pupil shift is smaller than the length of a sub-aperture, the pupil's central position P X and P Y can be calculated via Equations (1) and (2):  . optical path inside the ciao Wfs cryostat. after entering the entrance window, the beam passes through a field lens (mounted on an X-Y piezo stage -n-111 neXline oeM linear actuator) which can shift the pupil on two orthogonal directions, an achromatic doublet, and a filter (mounted in a filter wheel), before imaging the pupil on the lenslet array and the detector.
24 outer edge sub-apertures of the 9 × 9 Shack-Hartmann sensor are used to measure the pupil shift. Due to the internal optical alignment of the wavefront sensor as well as the stroke of the pupil actuator, we could only introduce a lateral pupil shift between -0.4 and 0.4 in the length of a sub-aperture on the SHWFS, which corresponds to a shift of -0.36 to 0.36 m on the VLT's primary mirror.
To verify the performance of the algorithm, we conduct the following experiment: we move the pupil's centre along Y-axis from -0.4 to 0.4 sub-aperture while keeping the pupil's centre on X-axis fixed at different values. We measure 200 frames of fluxes in the SHWFS at each position, and then use the intensities of all the 24 outer edge sub-apertures of each pupil's position to calculate pupil's shift P Y . Finally, we compare P Y with the applied shift to evaluate the performance of this algorithm. To compensate the deviation in the measurement of the pupil shift caused by the non-uniformity of the illumination in the pupil, a flat-field calibration is implemented. The measurement of pupil shift on X-axis is similar to the measurement of pupil shift on Y-axis. The result of the experiment is shown in Figure 4. Figure 4 shows that the measurement of the pupil shift on Y-axis is quite precise when 'shiftx' = 0 (the dot-dashline), the maximal deviation is less than 0.05 sub-aperture. The deviation of the measurement gets larger as 'shiftx' departs from 0.
We adopt the maximal deviation between the measured pupil shift and the applied pupil shift to be the indicator to evaluate the performance of the algorithm. The result is shown in Figure 5. As the range of the applied pupil shift on Y-axis shrinks from ±0.4 (-0.4 to 0.4 sub-aperture) to 0, while keeping 'shiftx' at 0, the maximal deviation (dotted curve) reduces from 0.042 to 0, other curves which represent different 'shiftx' also have the same tendency. On the other hand, as 'shiftx' varies from -0.4 to 0, the average values of the curves also decline, which means the deviation   (3) and (4) is to normalize the value of the formula to 1 or -1 when the pupil shift is an entire sub-aperture in either positive or negative direction.

The factor 2 in Equations
We analyse the data with this algorithm to compare it with the former one which just utilizes the outer subapertures. The result is shown in Figure 6. Compared with Figure 4, although when 'shiftx' = 0, the precision of the measurement in Figure 6 is pretty high, it degrades when 'shiftx' departs from 0 -even approaches about 0.15 sub-aperture when 'shiftx' = -0.4 and 'shifty' = 0.4which means the precision of the algorithm degrades fast when the pupil shift departs from the axes. As a result, in general, the performance of the algorithm which only utilizes the fluxes in the 24 outer sub-apertures is better than the one which utilizes the fluxes in both inner and outer sub-apertures.

Integration time
Disturbances and noises such as the atmosphere turbulence (to simulate the atmosphere turbulence in the actual observation, an artificial turbulence with a seeing of 0.5 arcsec is added in the experiment), photon noise and detector noise can induce errors to the measurement of the pupil shift and these errors can be reduced The results of the measurements, when 'shiftx' varies from 0 to 0.4, are not shown here because they are basically the same with the results when 'shiftx' varies from 0 to -0.4. The maximum of the maximal deviations of the measured pupil shift in Figure 5 is 0.095 sub-aperture (the diamond block marked with '1st correction'). It means for an initial pupil shift within -0.4 to 0.4 sub-aperture on both axes, the residual pupil shift after one pupil shift correction will decrease to under 0.095 sub-aperture. Furthermore, this residual deviation decreases to less than 0.01 sub-aperture (the '2nd correction' diamond block) after performing the correction for the second time. In CIAO, the budget of the wavefront error caused by the misalignment between the DM and the SHWFS is 50 nm, which corresponds to a residual pupil shift after the correction of 10% of the length of a sub-aperture (9,10). As the smaller the pupil shift is, the better the CIAO system works, here we adopt a stricter demand: 5% of the length of a sub-aperture -this could be achieved within two correction steps with the intensities method when the initial pupil shift on both axes are within -0.4 to 0.4 sub-aperture.

Using additional sub-apertures
For completeness, we consider another similar algorithm in which the four inner sub-apertures are also utilized. The formulas of this algorithm are: We test the impact of the number of gathered frames in our experiment on the variance of the measured pupil shift for the intensities algorithms. The result is shown in Figure 7 (grey curve).
The grey curve in Figure 7 shows that the variance of the measurement for the intensities algorithms decreases as more frames are gathered. If we want to reduce 2 P X to by gathering more frames of data in the experiment (in other words, extending the integration time). Assume 2 P X is the variance of the measurement signal P X , to detect a pupil shift as small as 5% of a sub-aperture, 2 P X needs to be less than 10 −4 d 2 with d as the length of a sub-aperture (which means an average measurement error of about 1% sub-aperture).  experiment. As the length of a sub-aperture of the SHWFS in our experiment corresponds to a pupil actuator motion of 300 μm and we want the deviation caused by the pupil actuator to be smaller than 1% of a sub-aperture, the location precision of the actuator needs to be smaller than 3 μm. We test the error distribution for this actuator in the stroke -0.4 to 0.4 sub-aperture, the result is shown in Figure 9.
The error distribution of the actuator complies with a Gaussian distribution whose mean value is -0.08 μm and standard deviation is 0.21, thus the precision of the actuator is -0.50 to 0.34 μm with a confidence interval of less than 10 −4 d 2 , the least number of frames that should be gathered for the intensities method is 150. In fact, we set the number to 200 in our experiment. As the frequency of the AO closed-loop operation is 500 Hz, the integration time is 0.4 s.

Non-uniformity of the illumination of the pupil
The pupil in the experiment is not illuminated uniformly, this non-uniformity leads to deviation in the intensities method -it causes a deviation of 0.13 sub-aperture -when the pupil's centre is at the original point.
A flat-field calibration is implemented to eliminate the deviation and it records the intensities in the sub-apertures when the pupil's centre is at the original point. Every frame of intensity data has to be divided by this calibration frame before further process. Figure 8 shows the comparison of the measurements of the pupil shift before and after the flat-field calibration with 'shiftx' = 0. The deviation caused by the non-uniformity of the illumination in the pupil is significantly eliminated after implementing the flat-field calibration.

Precision of the pupil actuator
The pupil position located on the lenslet array of the SHWFS is moved via a pupil actuator based on an internal piezo walk mechanism. The precision of this actuator determines the precision of the applied pupil shift in the

Principle of the matrix method
The matrix method computes the pupil shift by comparing two interaction matrices (IM) of the DM which are measured, respectively, before and during the closed-loop AO operation. The former IM is also called the nominal IM (M0) which is obtained when the applied pupil shift is zero. Figure 10 shows the structures of matrix M0 in CIAO SHWFS, as each of the 68 sub-apertures has two tilts (X 95% -which is much smaller than 3 μm, as a result, the deviation caused by the pupil actuator can be neglected in our analysis.

Pupil rotation
Pupil rotation can lead to the change of the intensities' distribution in the SH sensor. But in our experiment, the piezo walk mechanism can only move the pupil along two orthogonal directions, therefore the pupil rotation is negligible.  where: M1 is the real-time IM (dimension: 136 × 60), Δ meas is the difference between two successive SHWFS measurements (dimension 136 × 1), Δ comm is the difference between two successive DM command vectors (dimension 60 × 1), Δ turb is the turbulence increment between the two successive frames We are able to measure a wealth of pairs (Δ meas , Δ comm ) by acquiring slopes and voltages data. As Δ turb is uncorrelated from Δ comm , it can be considered just like a 'noise' term. Thus, getting plenty of measurements allow us to solve Equation (9) for M1 by the least square estima- with N the number of frames acquired, M1 can be calculated: The experiment is conducted as follows: move the pupil's centre along Y-axis from -0.4 to 0.4 sub-aperture while keeping the pupil's centre on X-axis fixed at different values. We measure 2000 frames of the slopes in the 68 sub-apertures of the SH sensor and the voltages of the 60 actuators during the closed-loop control at each pupil position, then make use of the slopes and voltages to calculate M1 with Equation (10). Finally, we compute the pupil shift by comparing M1 with the nominal IM M0 using Equation (7).
The result of the experiment is shown in Figure 12: the matrix method works well when 'shiftx' is small, but its precision degrades fast as 'shiftx' departs from 0, the maximal deviation reaches more than 50% when 'shiftx = -0.4' . To evaluate the performance of the measurement, we adopt the maximal deviation between the measured pupil shift and the applied pupil shift to be the indicator. As the deviation of the measurement is too high, when the pupil centre gets far from the original point, we verify this algorithm's performance when the applied pupil shift is within only -0.2 to 0.2 sub-aperture on both axes. The result is shown in Figure 13. Figure 13 shows that the average values of the curves decline when 'shiftx' varies from -0.2 to 0, in addition the maximal deviations also decrease as the range of applied pupil shift on Y-axis ('shifty') shrinks from ±0.2 (-0.2 to 0.2 sub-aperture) to 0. The deviation of the pupil shift measured by this algorithm decreases fast after applying multiple pupil shift corrections. The two diamond blocks in Figure 13 show that an initial pupil shift within -0.2 to 0.2 sub-aperture on both axes can be reduced to 0.097 and Y directions), M0 can be decomposed into M0 X and M0 Y . Because there are 60 actuators on the back of the DM and each actuator has its own IM, both M0 X and M0 Y are 60 × 68 matrices. The DM's IM varies when the pupil shifts during the closed-loop AO operation. Let M1 be the real-time DM's IM obtained during the closed-loop AO operation, M1 has the same structure as M0 and can be decomposed into M1 X and M1 Y .
To apply the matrix algorithm, we need to differentiate the matrix M0 on X and Y directions to get two new matrices, dM0 X and dM0 Y , which represent M0's gradients on the two directions. Figure 11 shows the structures of dM0 X and dM0 Y . As with M0 X and M0 Y , both dM0 X and dM0 Y are 60 × 68 matrices.
The matrix method uses M0, M1 and dM0 to calculate the pupil shift. Assume a pupil shift (P X , P Y ) which is much smaller than the length of a sub-aperture and the corresponding real-time IM obtained in the AO closed-loop control is M1. Because the shift is small, we have: Equation (5) with least square estimation, P X and P Y can be written as: As P X or P Y is just a number instead of a matrix, Equation (6) can be written as:

Result of the experiment using the matrix method
We conduct the experiment of the matrix method on the CIAO bench in the lab. As shown in Figure 3, there are 68 sub-apertures in the SH sensor, each sub-aperture can measure the average wavefront slopes at X and Y directions in every frame. The deformable mirror (MACAO DM) has 60 actuators and we need to measure the voltages applied to the 60 actuators to calculate the real-time IM. Assume a DM command C comm and the SHWFS's measurement W meas after performing C comm , we have: where T turb is the atmosphere turbulence during the SHWFS's measurement and M1 is the real-time interaction matrix. Furthermore, we have:

Integration time
We test the impact of the number of gathered frames in our experiment on the variance of the measured pupil shift and 0.036 length of a sub-aperture after 1 and 2 pupil shift correction steps, which means the requirement on residual pupil shift -5% of the sub-aperture -can be achieved within two correction steps with the matrix method.  suitable for situation with faint object, low-frequency pupil correction and small pupil shift.

Conclusion
We analyse two pupil-tracking methods which can measure the lateral pupil shift with a SH sensor based on AO system. The first one utilizes the fluxes in all outer edge sub-apertures of a Shack-Hartmann sensor and the second one utilizes the real-time interaction matrix obtained during closed-loop operation of the AO system. We conduct two experiments with the CIAO system to verify these pupil-tracking algorithms. The results of the experiment show that both algorithms are effective, after two correction steps, the residual pupil shift is reduced to less than 5% of a sub-aperture. The intensities algorithm works better in situation with bright object, high-frequency pupil correction and large pupil shift, while the matrix algorithm is more suitable for situation with faint object, low-frequency pupil correction and small pupil shift.

Disclosure statement
No potential conflict of interest was reported by the authors.
for the matrix algorithm. The result is shown in Figure 7 (black curve).
The black curve in Figure 7 shows the variance of the measurements for the matrix algorithm decreases as more frames are gathered. To reduce 2 P X to less than 10 −4 d 2 , the least number of frames that should be gathered is 1100. In fact, we set the number to 2000 in our experiment. As the frequency of the AO closed-loop operation is 500 Hz, the integration time is 4 s.

Pupil rotation and the precision of the pupil actuator
Pupil rotation can lead to the mis-registration seen by the interaction matrix of the AO system. But in our experiment, the piezo walk mechanism can only move the pupil along two orthogonal directions, thus the pupil rotation is negligible.
The pupil actuator used in the experiment for the matrix method is identical with the one used for the intensities method. Thus, the result of the test shown in 2.4.3 is still valid here, which means the deviation caused by the pupil actuator can be neglected in the experiment of the matrix algorithm. Table 1 shows the comparison of the two algorithms. Firstly, in our experiment, the deviation of the matrix algorithm reaches more than 50% when the pupil shift's range exceeds -0.2 to 0.2 sub-aperture, while the intensities algorithm can keep the deviation lower than 25% within the range -0.4 to 0.4 sub-aperture, thus the intensities algorithm has a larger valid range. Secondly, compared with the matrix algorithm's 4 s integration time, the intensities algorithm only needs a 0.4 s integration time, which makes it more suitable for those situations which need high-frequency pupil corrections. Thirdly, unlike the intensities algorithm, the matrix algorithm does not directly utilize the fluxes in the sub-apertures, thus it has a lower requirement on the brightness of the observed object and is more suitable for observing faint science objects.

Comparison of the two algorithms
Generally, the intensities algorithm works better in situation with bright object, high-frequency pupil correction and large pupil shift, while the matrix algorithm is more