A remote sensing assessment index for urban ecological livability and its application

ABSTRACT Remote sensing provides us with an approach for the rapid identification and monitoring of spatiotemporal changes in the urban ecological environment at different scales. This study aimed to construct a remote sensing assessment index for urban ecological livability with continuous fine spatiotemporal resolution data from Landsat and MODIS to overcome the dilemma of single image-based, single-factor analysis, due to the limitations of atmospheric conditions or the revisit period of satellite platforms. The proposed Ecological Livability Index (ELI) covers five primary ecological indicators – greenness, temperature, dryness, water-wetness, and atmospheric turbidity – which are geometrically aggregated by non-equal weights based on an entropy method. Considering multisource time-series data of each indicator, the ELI can quickly and comprehensively reflect the characteristics of the Ecological Livability Quality (ELQ) and is also comparable at different time scales. Based on the proposed ELI, the urban ecological livability in the central urban area of Wuhan, China, from 2002 to 2017, in the different seasons was analyzed every 5 years. The ELQ of Wuhan was found to be generally at the medium level (ELI ≈0.6) and showed an initial trend of degradation but then improved. Moreover, the ecological livability in spring and autumn and near rivers and lakes was found to be better, whereas urban expansion has led to the outward ecological degradation of Wuhan, but urban afforestation has enhanced the environment. In general, this paper demonstrates that the ELI has an exemplary embodiment in urban ecological research, which will support urban ecological protection planning and construction.


Introduction
In recent years, rapid economic growth has enabled a significant improvement in people's quality of life (Tang et al. 2017).However, the accompanying urbanization has also caused problems, such as social pressure (Graafland and Smid 2017), environmental degradation (Adams and Klobodu 2017), and declining health (Liu et al. 2017b).As a global phenomenon, urban sustainability also seeks a balance between local socioeconomic benefits and eco-environmental benefits (Yan et al. 2018).From this perspective, urban livability could be regarded as a local subset of sustainability (Alijani et al. 2020).Moreover, some environmental, economic, and social aspects are usually contained in urban livability (Kashef 2016).In this case, some sustainability concepts, such as living quality and living environment, have also emerged.Nevertheless, urban livability differs from the wellbeing of residents.The concept of well-being often refers to many quantitative indicators, such as income, health, public service, and environmental quality, as well as qualitative indicators, such as pleasure and joy (Cummins 1995;Wei and Gao 2017).
The concept of urban ecological livability focuses more on the surrounding environmental quality of local residents, which is related to the characteristics of places, environmental comfort, and the ecology of communities (Aulia 2016;Liu et al. 2017a).It is obvious that local dwellers may embrace a wonderful life with a healthy natural environment.In summary, one pivotal emphasis of urban ecological livability is the livable urban ecological environment (Fu, Yu, and Zhang 2019;Valcárcel-Aguiar and Murias 2019).Timely and comprehensive evaluation of the livable urban ecological environment and its change trends is the basis of urban planning and human habitation, and the key to sustainable urban development.To better understand the status and change trend of the urban ecologically livable environment, various urban ecological evaluation systems and monitoring technologies have been applied to evaluate the quality and function of the ecological system to formulate ecological protection countermeasures (Chang et al. 2017;Cobbinah, Poku-Boansi, and Peprah 2017;Grubert 2018).
Due to a lack of flexibility in traditional monitoring systems, which also have a very high cost in human, material, and financial resources, the data provided by traditional environmental monitoring systems insufficiently reflect the change trends of the urban ecologically livable environment.Remote sensing not only helps regulators promote the monitoring of environmental changes (Patino and Duque 2013;Bhardwaj et al. 2016;Fassnacht et al. 2016;Mushore et al. 2017) but also facilitates urban environmental research.Furthermore, remote sensing can help us formulate feasible plans to provide ecological baseline maps for the early development of cities (Willis 2015;Chen et al. 2018c;Dlamini et al. 2019).Remote sensing data are now widely utilized in the field of ecology, due to the advantages of large-area, real-time, rapid, and periodic repeated observations (de Araujo Barbosa, Atkinson, and Dearing 2015;Rocchini et al. 2017).
Researchers have applied different types of remotesensing data to analyze various aspects of the urban ecologically livable environment (Vogt et al. 2015;Wang and Zhao 2016;Li et al. 2017).However, researchers have tended to focus on a single aspect of the urban ecological environment.Vegetation indices and impermeable surface coverage have been employed to evaluate urban ecology (Jiang et al. 2006;Gupta et al. 2012;Yu et al. 2016;Zhang, Zhu, and Fan 2016), land surface temperature (LST) has been utilized to evaluate the urban thermal environment (Buyantuyev and Wu 2010;Shen et al. 2016), and air quality indices have been applied to represent urban air pollution (Zhan et al. 2018;Zhou et al. 2018).These studies have tended to analyze environmental changes from the following perspectives: regionalscale range from global to national, to watersheds, and even to metropolitan areas; data scale from statistical data to remote sensing data (Gu et al. 2015;Li, Fang, and Wang 2016;Huang et al. 2017;Nguyen and Liou 2019;Shao et al. 2021a).
Although the effectiveness of remote sensing data for multi-scale environmental monitoring is constantly improving, current studies are still limited by the inconsistency of the imaging period and the spatial resolution of the remote sensing images (Shen et al. 2015;Shao, Wu, and Li 2021b).For example, some urban indices based on remote sensing data are usually extracted with only a single scene from each year, due to the imaging time, clouds, missing pixels, and other degradation factors, disregarding the differences of each remote sensing index in the different seasons (Castrence et al. 2014;Li, Gong, and Liang 2015;Alqurashi and Kumar 2016).
Each indicator reflects the urban ecosystem in different ways, and they cannot be viewed separately (Xu et al. 2018).Therefore, different scholars have begun to extract a variety of complex urban environmental factors and have developed different comprehensive indices for efficient and real-time evaluation of the urban ecologically livable environment.These indices or frameworks include the Ecological Footprint (EF) (Lin et al. 2018), the Environmental Performance Index (EPI) (Hsu, Lloyd, and Emerson 2013), the Ecological Index (EI) (Rsei 2013), the Remote Sensing based Ecological Index (RSEI) (Hu and Xu 2018), the Comprehensive Evaluation Index (CEI) (He et al. 2017), the Urban Environmental Quality Index (UEQI) (Musse, Barona, and Santana Rodriguez 2018), the theme-based sustainability framework developed by the United Nations Commission on Sustainable Development (UNCSD) (Huang, Wu, and Yan 2015), the Procedure for Ecological Tiered Assessment of Risk (PETAR) (Xu et al. 2016), and multi-dimensional indicators within a Drivers, Pressures, States, Impacts, and Responses (DPSIR) framework (Liu, Song, and Deng 2019).These studies are mainly based on natural resources, environmental performance, or ecological factors.Although the comprehensive indices and frameworks of these studies are relatively mature, they still face some problems.For example, the indicators and standards are diverse, and different types of data are combined in various evaluation studies, resulting in relatively complex data acquisition, processing, and calculation standards.On the other hand, due to the difficulty of selecting indicators for calculating urban remote sensing ecology, some studies have either set up indices with less coverage or set up complex indices that are not suitable for small areas, such as inner cities.For the RSEI index, the effects of the atmosphere are not taken into account; factors such as urban dryness and humidity are not covered by the CEI; and the data employed in the evaluation frameworks are difficult to obtain, resulting in a real-time decline in the acquisition of urban ecological conditions.Thus, the construction of a comprehensive, remote sensing-based, ecological index could become popular because of its advantages of being a rapid and comprehensive reflection of the urban ecologically livable environment.
In addition to the establishment of an ecological index system, another crucial issue for the evaluation of urban ecological livability is to weigh each factor according to its relative effects on ecological livability.Various methods have previously been utilized and developed for urban eco-environmental assessment.Equal weighting, a simple assigning method, tends to disregard the validity and transparency of the indicators (Rowley et al. 2012).The Analytic Hierarchy Process (AHP) is a quantitative evaluation method for a multi-index system that can decompose a complex problem into several layers or factors (Song et al. 2010).However, the weight results are easily influenced by the cognition level of experts, as well as the subjective consciousness.Although the Principal Components Analysis (PCA) takes into account both qualitative and quantitative analysis, it is also affected by ambiguous principal components or target quantity (Mikulić, Kožić, and Krešić 2015).There are other patterns, such as gray or fuzzy evaluation (Phillis, Kouikoglou, and Verdugo 2017), artificial neural-network evaluation method (Park et al. 2004), and the benefit of the doubt approach (Van Puyenbroeck and Rogge 2017).Nevertheless, the complex parameters in these models are not always easy to obtain and calculate.The objective weighting method which is widely practiced in many research fields, often produces satisfactory results in decision analysis (Lan, Wu, and Lee 2012;Alemi-Ardakani et al. 2016;Sahoo et al. 2016), while the entropy method is an important allocation method (Zhao et al. 2018).
Hence, to better evaluate the regional, ecologically livable environment, it is necessary to develop an objective and comprehensive urban ecological livability quality index, which is referred to as the Ecological Livability Index (ELI).In addition, a spatiotemporal continuous basic dataset framework to monitor the urban regional environment is presented.Furthermore, the status of obtaining evaluation results derived from only one snapshot from remote sensing images of in situ reality may be broken via ELI.To address these issues, the main objectives of this paper are: 1) to establish a regional evaluation index of Ecological Livability Quality (ELQ) from five ecological aspects, i.e. greenness (Normalized Difference Vegetation Index, NDVI), temperature (Land Surface Temperature, LST), dryness (Normalized Difference Built-up and bare Soil Index, NDBSI), water-wetness (Near-Water Distance, NWD), and atmospheric turbidity (Aerosol Optical Depth, AOD), integrating multi-source, temporally continuous, and spatially fine remote sensing data and 2) to verify the universality of the ELI and explore the spatial and temporal changes in different seasons and years in Wuhan, China.
The rest of the paper is organized as follows: In Section 2, the study area and data are described.In Section 3, the generation of the basic data sets for remote sensing ecological evaluation based on spatiotemporal fusion are introduced, and the evaluation method for urban ecological livability quality and the construction of the ELI are presented.Section 4 describes the spatiotemporal results of the ELI in Wuhan.Some trend reasons, feasibility, and directions for the ELI are discussed in Section 5, while conclusions are shown in Section 6.

Study area
Wuhan, the capital of Hubei Province in central China, is the only subprovincial city and megacity in the six central provinces and is the core city of the Yangtze River Economic Belt.The city is located between 11341'−11505' E and 2958'−3122' N. By the end of 2018, the total area of Wuhan was 8569.15 km 2 , accounting for approximately 4.6% of Hubei Province (the seven central urban areas covered approximately 955.15 km 2 , and the built-up areas covered approximately 812.39 km 2 ).
In this paper, the study area covers the seven central urban areas of the city of Wuhan (abbreviated as WHC), including Wuchang (Wuchang District, Qingshan District, and Hongshan District), Hankou (Jianghan District, Jiang'an District, and Qiaokou District), and Hanyang (Hanyang District).The three areas in the Yangtze River and Hanjiang River interchange form the "three towns of Wuhan" (Figure 1).

Data
In this study, two kinds of remote sensing datasets were used to evaluate the ecological livability quality of the WHC area.Landsat surface reflectance products in 2003, 2006, 2007, 2008, and 2011(Landsat 5 TM) and in 2013, 2016, 2017, 2018 (Landsat 8 OLI/TIRS) were downloaded from the United States Geological Survey (USGS, https://earthexplorer.usgs.gov/).The NDVI, LST, NDBSI, Normalized Difference Water Index (NDWI), and Modified Normalized Difference Water Index (MNDWI) were calculated from the Landsat images.The standard MODIS land product (MOD09A1, MOD11A2, MOD13A3, and MAIAC AOD) provided by the National Aeronautics and Space Administration (NASA) was adopted, and all MODIS data corresponding to the time of Landsat data, as well as data for a five-year interval from 2002 to 2017, can be downloaded from the NASA website (https://ladsweb.modaps.eosdis.nasa.gov/search/).The specific date of Landsat and MODIS data pairs are listed in Table 1.

Methods
The study followed the steps shown in Figure 2. To better understand the changes in the urban ecologically livable environment in different years and seasons, first, it was necessary to construct a temporally continuous and spatially fine basic dataset.Thus, the processing of the MODIS images involved the construction of a Terra-Aqua satellite data composite, harmonic analysis, missing information reconstruction, and quarterly data synthesis.The NDVI, LST, NDBSI, NDWI, and MNDWI were obtained through Landsat images.The corresponding bands and parameters of the MODIS and Landsat data were then taken as references for spatiotemporal fusion and band calculation to obtain the basic remote sensing evaluation dataset, i.e.NDVI, LST, NDBSI, NWD, and AOD.Furthermore, five remote sensing ecological indicators were selected to reflect the quality of the local ecological environment, which roughly cover the greenness, temperature, dryness, water-wetness, and atmospheric turbidity of urban ecology.Second, in importance, we constructed an appropriate evaluation index -the ELI -which was geometrically aggregated with the abovementioned indicators by non-equal weights based on an entropy method, to reflect the urban ecological livability quality.More details are provided in the following sections.

Remote sensing dataset generation for ecological assessment
In this study, we processed remote sensing images of the WHC area with an interval of every 5 years from 2002 to 2017 and generated remote sensing ecological evaluation data for each year with a quarterly timescale and spatial resolution of 30 m.

Reconstruction and fusion of missing information
Due to the missing information caused by cloud cover (Li et al. 2019), sensor faults, and the discontinuity of the spatial and temporal distribution of single-site data, we mainly addressed the surface reflectance and meteorological data of the MODIS and Landsat images.For the meteorological data, seamless AOD daily data were obtained by mean synthesis and twosatellite fusion of Terra and Aqua, and then mean data of each year and season were generated by averaging the AOD data.In addition, the harmonic analysis method was utilized to remove cloud interference and fill the data gaps (Yang et al. 2015) of the MODIS data.It represents the time-series data of the reflectance or LST as a combination of multiple sines and cosines in Equation ( 1), and uses several important harmonics smoothing data to obtain the continuous results.
where Y T ð Þ represents the result after reconstruction; A 0 is the average of time-series data; A i and B i are parameters of sines and cosines, respectively;F i represents the frequency;T is the number of points in time series while N is the length of time series.
Moreover, according to the reconstructed data, the quarterly mean values of each year were obtained for the subsequent spatiotemporal fusion.All the satellitederived mean data results were considered and utilized for cloud-free conditions.Contrary to each index from only one snapshot or a single period of the in-situ reality, the average value avoids instantaneity during the evaluation and reduces randomness in the short term.

Spatiotemporal fusion of Landsat and MODIS
The remote sensing spatiotemporal fusion method based on nonlocal filtering (Cheng et al. 2017) was used to fuse the Landsat, MODIS surface reflectance, LST, NDVI, and NDBSI for the basic remote sensing ecological evaluation dataset.
Considering the nonlocal characteristics of similar pixels, this fusion method optimizes the selection of similar pixels and the construction of a weight function.For high-and low-resolution images, it is assumed that the change in reflectance from these two data during the same period (t k to t P ) is linear, disregarding the difference in atmospheric correction and anthropogenic factors.During the spatiotemporal fusion, it is necessary to ensure that the coordinate system and spatial scope of the two resolution images are consistent.Thus, all the regression coefficients are assumed to be unchanged.The details are from Equation (2) to (3) .
where L x i ; y j ; t P À � and L x i ; y j ; t k À � are the surface reflectance of Landsat image pixels whose coordinates are (x i ,y j ) at time t P and t k , respectively; M x i ; y j ; t P À � and M x i ; y j ; t k À � are the surface reflectance of MODIS image pixels whose coordinates are (x i ,y j ) at time t P and t k , respectively; A and B represent the linear regression coefficients of the pixel from time t P to t k , respectively.
In the real experiment, the categories of urban objects and each pixel of its remote sensing images may change significantly over time.Therefore, auxiliary information needs to be introduced to meet the accuracy requirements of fusion.Here, similar pixels in each fusion reference image are obtained in a local window range, and the weight function is calculated for the weighted average in Equation ( 4).
where w is the window size; L x w=2 ; y w=2 ; t P À � is the surface reflectance of the high-resolution Landsat data obtained from the center of the window at the predicted time and w ijk represents the weight of each similar pixel, consisting of individual and whole weights.
To ensure the possibility of ecological environmental change in the different years and seasons and the accuracy of the fusion, we attempted to choose images that were close to the year and consistent with the quarter (Table 1).The corresponding remote sensing parameters for the urban ecologically livable environment were selected for each ecological aspect.The calculation and normalization methods for each index are expressed as follows.

Greenness
Green vegetation interacts with urban ecology in many aspects and is a necessary condition for promoting the urban environment and maintaining ecological balance.The most direct manifestation of greenness in remote sensing data is a vegetation index.The NDVI, which can reflect the background influence of the plant canopy, is the most common vegetation index for the analysis of plant growth state and spatial distribution density.For the monitoring of vegetation growth, the NDVI is commonly used to reflect the growth pattern of vegetation.When the value of the NDVI is positive, the larger the value of the NDVI, the greater the vegetation coverage.Moreover, the NDVI has been successfully applied to examine the temporal and spatial variation trends and dynamic distribution of vegetation (Xu and Zhang 2013).The NDVI calculation method and normalization method are defined in Equations ( 5) and ( 6) : where ρ NIR and ρ Red are the surface reflectance values of the near-infrared and red bands, respectively; NDVI i is the value of the NDVI in pixel i; NDVI max and NDVI min are the maximum value and minimum value, respectively, of the NDVI in the study year; and NDVI i represents the normalized value of the NDVI in pixel i.

Temperature
Air temperature and LST, two strongly correlated variables, are usually used as complementary parameters to measure urban thermal environment (Jin and Dickinson 2010;Gallo et al. 2011;Mildrexler, Zhao, and Running 2011;Benali et al. 2012;Marzban, Sodoudi, and Preusker 2018).Meteorological stations can provide fine air temperature data with a height of 2 m, but it is difficult to completely reveal the spatial heterogeneity of urban thermal environment (Lin et al. 2012).Meanwhile, LST, a parameter that is now often used to estimate air temperature, may complement that (Yuan et al. 2020).In view of the characteristics of spatial continuity, surface physical processes at regional or global scales and lower boundary conditions to ecological livability assessment, LST is a better factor for the surface energy balance and greenhouse effect than air temperature, especially in urban areas (Mutiibwa, Strachan, and Albright 2015;Shen et al. 2020; Venter Zander, Chakraborty, and Lee 2021).However, there are many LST retrieval algorithms.In this study, the radiative transfer equation algorithm was selected to obtain the LST from the Landsat images (Sobrino, Jiménez-Muñoz, and Paolini 2004).
The LST shows obvious differences in different seasons.Considering the living environment, temperature is particularly important for urban dwellers.When the temperature of a certain area is similar to the appropriate temperature, it is determined to be more livable.Therefore, we set a comfort temperature as the reference temperature (LST R ) and calculated the difference between LST and LST R so that we could obtain the positive or negative effect of this index by judging the deviation degree between LST and LST R .We utilized the following temperature normalization method in Equation ( 7): where LST i represents the normalized LST in pixel i; LST i is the LST in pixel i; and LST max and LST min are the maximum LST and minimum LST, respectively, in the study year after comparison with the reference temperature.LST R is the reference temperature, which may vary from city to city.Some studies have suggested that the "human thermal sensation" is comfortable and without thermal pressure between 23 and 27°C (Omonijo 2017).Generally, LST is slightly higher than the local temperature, and 25°C was selected as the appropriate temperature basis for Wuhan in this study.

Dryness
The boom in urban infrastructure, such as buildings, squares, roads, and parking lots, is a typical characteristic of urbanization.Population growth, impervious surface increase have been demonstrated to affect urban landscapes and ecological livable in a certain extent in previous work (Ellis et al. 2006).The original natural landscapes dominated by vegetation are gradually substituted for numerous artificial impervious surfaces.Moreover, a considerable part of bare soil, which causes a "dry" surface.The negative effects of human interventions, project construction, and impervious surface reflect in those poor ecoenvironment quality areas (Zhang, Zhao, and Gu 2014;Zhang, Zhu, and Fan 2016;Shan et al. 2019).Therefore, understanding impervious surfaces and bare soil is also very important for urban ecological planning.The NDBSI integrates the Index-based Built-up Index (IBI) and Soil Index (SI) in Equations ( 8) and ( 9) to better combine their respective advantages and produce a more complete urban ecological dryness index (Xu 2008).The NDBSI calculation method and normalization method are expressed in Equations ( 10) and ( 11) (Hu and Xu 2018): where ρ Red , ρ Green ,ρ Blue ,ρ NIR , andρ SWIR1 represent the surface reflectance values of the red, green, blue, nearinfrared, and short-wave infrared bands, respectively, from the MODIS and Landsat satellites;NDBSI i is the NDBSI in pixel i; NDBSI max and NDBSI min are the maximum NDBSI and minimum NDBSI, respectively in the study period; and NDBSI i represents the normalized NDBSI in pixel i.

Water-Wetness
For urban ecological planning, the protection of wetlands, lakes, and other water resources can not only meet the needs of public and social development but also help to regulate the ecological balance of urban areas while improving the quality of the ecological environment and promoting sustainable urban development (Liu et al. 2013;Feyisa et al. 2014).In addition, people may have some real and substantial psychological ripple effects from viewing or touching surrounding environment (Yu et al. 2016).Therefore, it is indeed important to understand and measure open surface water accessibility because of its entertainment or other function after urban dwellers' visual and physical experience.Thus, a water-wetness index is utilized to represent water-wetness in the living environment.
Before calculating the water-wetness, it is necessary to extract the water range of the urban area through a waterbody index, such as the NDWI or MNDWI in Equation ( 12) and ( 13).After the segmentation and extraction of the waterbody range by these two methods whose thresholds are both greater than zero (Xu 2006), the Euclidean distance from each pixel of a non-water body area to the boundary of the waterbody range in the urban area is obtained (Salonen et al. 2012).
whereρ Green , ρ NIR and ρ MIR are the surface reflectance values of the green, near-infrared, and midinfrared bands, respectively.
For brevity, the shortest distance from any pixel of a non-water body area to the boundary of the waterbody vector is referred to as the NWD (Figure 3).The WHC region is rich in water resources and has a wide range of waterbodies.
The areas closer to the water appear light blue, while the areas further from the water appear orange.In the WHC area, these areas further from open water, namely, orange areas, are usually builtup areas or have large amounts of bare soil.The low ecological quality in these areas is mainly due to the human high-stress intensity activities and large-scale urban expansion (Zhu et al. 2021).Regions with water relatively nearby, such as lakes and rivers, may be considered more livable and have a relatively high ecological livability quality (Soranno et al. 2015;Zachariasz and Porada 2019).Hence, we select a Reference Water Distance (NWD R ) and calculate the Difference Distance between the NWD and the NWD R (DNWD) so that we can obtain the positive or negative effects of this index by judging the deviation degree between the two distances.As the distance between the urban area and the water reaches a certain range, the impact on it will not substantially change in different cities, so we select one threshold distance (NWD T ).When the DNWD exceeds this threshold, we can set it as a constant distance (NWD C ). NWD i is calculated in Equation ( 14).
where NWD i is the NWD in pixel i; DNWD max and DNWD min are the maximum DNWD and minimum DNWD, respectively; NWD i represents the normalized NWD in pixel i; and NWD R is the reference distance, which may vary from city to city.In this study, we selected 100 m as the NWD R basis for Wuhan and set 1000 m as the NWD T to reasonably show the influence of NWD on the ELI.

Atmospheric turbidity
As an important part of the environment, the atmosphere plays a key role in environmental health.The relation of the AOD with regional air quality and public health is widely discussed as one of the main stressors of global climate change (Chin et al. 2014;Huang et al. 2014;Wang et al. 2017;Setti et al. 2020;Madadi et al. 2021).Moreover, the AOD is a reference parameter of aerosol with related to other air pollutants for the comprehensive characterization of atmospheric turbidity in urban ecological livability, and one of the most influential parameters for satellite-based PM modeling (Gupta et al. 2013;Chen et al. 2018aChen et al. , 2018b)).Air pollutants from large areas and long time series remote sensing inversion methods have been estimated or obtained in many studies, including AOD inversion (Ghotbi, Sotoudeheian, and Arhami 2016;Zang et al. 2017;Soni, Payra, and Verma 2018).In terms of the data applicability for intra-urban (Hammer et al. 2022), the correlation and modeling with other pollutants (Gupta et al. 2006;Lu, Zhang, and Streets 2011;Zang et al. 2018), the AOD is utilized as a reference parameter for air pollution in urban ecological livability.Generally, the higher the value of the AOD, the greater the altitudinal accumulation of aerosols, and the corresponding reduced visibility of the atmosphere.This finding reflects the turbidity of the atmosphere and usually has a negative impact on urban ecological livability.The following normalization method in Equation ( 15) is employed to obtain AOD i : where AOD max and AOD min are the maximum AOD and minimum AOD in the study period, respectively; and AOD i represents the normalized AOD in pixel i.

ELI construction
From the perspective of the five main environmental themes (freshwater; land; atmosphere; biodiversity; and oceans, seas and coasts) in the UNCSD framework (Huang, Wu, and Yan 2015), the ELI can also be composed of five indicators derived from remote sensing datasets at the pixel level.Unfortunately, the last two elements of the UNCSD framework are not easy to apply in urban areas.With regard to a more comprehensive and reasonable index, another benchmark, referred to as the RSEI based on the Pressure-State-Response (PSR) framework for the ELI, was employed (Hu and Xu 2018).First, the main pressures on landscape ecological security during urbanization are primarily derived from urban expansion and population growth.Therefore, some remote sensing building indices can be employed to explain the intensity of pressure on the landscape from external social activities.Second, vegetation coverage and its growth activity calculated by vegetation indices inversely reveal the status quo and trend of ecosystems.Last, temperature, air quality, and open surface water indicate responses to urban changes in the thermal, atmospheric, and water environments.
Based on this index evaluation framework, the NDVI, LST, NDBSI, NWD, and AOD were selected to generate the spatiotemporally continuous evaluation data of ELQ with a spatial scale of 30 m and a seasonal timescale.These indicators represent five ecological aspects of urban areas: greenness, temperature, dryness, water-wetness, and atmospheric turbidity.Ecological problems such as noise and other waste pollution were not considered in this study as it is difficult to obtain effective indicators from remote sensing data at the urban scale.The ELI is regarded as being affected by changes in these five ecological aspects.Therefore, we established the following method to calculate the ELI in Equation ( 16 The change in ELQ in urban areas can be evaluated by the ELI.When constructing a remote sensing evaluation index for the urban ecologically livable environment, normalization of each index can ensure that the ELI is an increasing function, that is, the higher the ELI value is, the better the ecological livability quality.In addition, this approach can ensure that the weight imbalance of regions with no water affects the universal applicability of the ELI.The methods used to determine the weights of this comprehensive evaluation index system include both subjective and objective weighting methods.The entropy method is an objective weighting method with high reliability and accuracy, which can reflect the utility value of the information entropy of each index.Thus, in this study, the entropy method was selected to determine the index weights to eliminate subjective preference (Chen, Lu, and Zha 2010).We calculated the ELI according to the normalization principle and the geometric mean with the weights of these five indicators (Gan et al. 2017).Therefore, the calculation of the ELI was conducted by Equation ( 17): where ELI i is the value of ecological livability quality at pixel i.The larger the value of ELI i is, the better the ecological livability quality and the higher the ecological livability.NDVI i , LST i , NDBSI i , NWD i , and AOD i represent the normalized values in pixel i; w 1 ; w 2 ; w 3 ; w 4 ; w 5 are the corresponding weights of each index, and the sum of their values is one.The weight is calculated by the pixel number of each indicator.The main steps are listed from Equations ( 18) to ( 21): Calculate the proportion of the ith pixel under the jth indicator: Calculate the information entropy of the indicator: Calculate the redundancy of the information entropy: Calculate the weight of each indicator: where x ij is the result of the homogenization of the heterogeneous indicators, and the normalization method is calculated with reference to each previous indicator.In addition, m is the number of indicators, and its value is five.n is the number of pixels for each indicator.
According to Equation ( 17), we can determine that the ELQ result calculated by the five indices will fall between one and two.To better measure and compare the changes in ELQ during the different study periods, the results should be further normalized.The maximum value of normalization is two, and the minimum value is one.The calculation is conducted in Equation ( 22): where ELI n i is the final remote sensing assessment index of ecological livability quality, ranging from zero to one.When the resulting ELI values in a certain region are nearly one, the ELQ of this area is superior and the ecological livability is high and vice versa.

Spatiotemporal characteristics analysis of the ELI
Studying the distribution characteristics and changing trends of urban ecological livability is a key prerequisite for urban construction or planning, and comprehensive and coordinated development.In this study, we quantitatively illustrated the characteristics of the ELI in urban areas from the spatiotemporal perspective with statistical and spatial heterogeneity analyses.The results were applied to investigate whether the ELI is practical for areas, such as the WHC area.The analyses paid more attention to the two aspects of temporal variation trend and spatial clustering, i.e. the overall distribution and variation trend of the ELI at annual and seasonal scales and the local pattern of the regional urban ecological livability quality.
Most trends over time can be intuitively obtained from the statistics of each grade of the ELI or its visual results.Moreover, spatial autocorrelation analysis can be performed to study the spatial correlation of a certain attribute at an adjacent location from the ELI, including global autocorrelation and local autocorrelation.To explore the spatial pattern or distribution characteristics of the ELQ in the WHC area, a spatial autocorrelation index was applied to measure its spatial correlation degree.Simultaneously, this paper gives more consideration to the local characteristic differences in the distribution of the ELI.The Local Indicators of Spatial Association (LISA) index (Anselin 1995) was utilized to further measure the degree of spatial association between one region and its adjacent spatial units of the ELQ in the WHC area, that is, the local Moran's I index in Equation ( 23), rather than Moran's I, was employed to determine whether the phenomenon of variable agglomeration existed and to indicate the hot and cold spots of the WHC area in this study.All the results were calculated using the GeoDa program (Anselin, Syabri, and Kho 2006).
where I is the local Moran's I index; i and j are any one of the pixels in the region; x i and x j are the ELI of pixels i and j, respectively; � x is the average ELI of a type of pixel; S 2 ¼ , which represents the ELI variance of pixel I; n is the total number of pixels; and w ij is the element of row i, column j in the spatial weight matrix.
In addition, the Standard Deviation Ellipse (SDE), which is a common tool for the geoscience analysis of directional characteristics (Ghasemi, Hamzenejad, and Meshkini 2018), was also employed to measure the centrality, spreading, direction, and spatial shape of the ELI.The center of the ellipse is the average center of the spatial distribution, and the standard deviation in the X and Y directions represents its long axis and short axis, representing the dispersion degree of the main direction and the secondary direction of the spatial element distribution, respectively, that is, the larger the length is, the more discrete it is.The azimuth of the ellipse reflects the main trend of its distribution.The area of the ellipse indicates the degree of concentration or dispersion of the spatial distribution of the geographical elements.
The larger the area, the more dispersed the spatial distribution of the good ecological livability quality.In contrast, the more concentrated the ecological livability quality is, the better its ecological environment is.The calculation formula for the SDE average center is presented from Equations ( 24) to (25): where (x i ,y i ) represents the spatial location of the research object, w i represents the weight, and (X w ,Y w ) represents the weighted average center.

Basic dataset generation for the urban ecologically livable environment
First, a spatiotemporal information fusion method based on a nonlocal means filter (Cheng et al. 2017) was utilized to provide basic data for the evaluation of the urban ecologically livable environment.Figure 4 shows the real and fused Landsat 8 surface reflectance results for the near-infrared band in the WHC area.The fused Landsat images have a highly similar spatial distribution to the real Landsat images, as shown in the water information and features of the Yangtze River and East Lake, indicating that this fusion method can be applied to the generation of regional ecological environment evaluation datasets.According to the quantitative results (Table 2.), the Root-Mean-Square Error (RMSE) of the green, red, and near-infrared bands are all below 0.05, while the correlation coefficient (R 2 ) of the green band reaches 0.7246 and that of the red and near-infrared bands exceeds 0.79.Furthermore, given the complexity and uncertainty of the various indicators for ELI, we conducted a preliminary comparison and verification of the fusion results of the three indices, namely, NDVI, LST, and NDBSI.The results of their comparison with real data are presented as follows: Figure 5 shows the real and fused difference results for these three indicators in the WHC area.Figure 5 3.), the RMSEs of the NDVI and NDBSI are 0.1147 and 0.0690, respectively, while that of the LST exceeds three.The correlation coefficient (R 2 ) of the NDBSI was approximately 0.6938, and the R 2 of the NDVI and LST exceed 0.7.

Weight of each index in the ELI of the WHC area
The weight of each index in the ELI for the WHC area was determined according to the entropy method.The specific results are shown in Table 4.A comparison of the weight values of each index reveals no obvious difference in the quantitative characteristics between two seasons, most of which float at approximately 0.2.Among them, the largest weight value in spring and autumn was the LST, while the minimum value in the former season was the NWD, and the weight value of the NDBSI was the smallest in winter.Furthermore, the maximum weight value in summer and winter is the NDVI, whereas the minimum of the former season is the NDBSI, while the latter is the LST.The changes in time (different seasons or years) and space of the ELQ in the WHC area were then obtained according to the established evaluation system and the results of the ELI.

Temporal change in the ecological livability quality in the WHC area
Figure 6 represents the temporal change trend of the ELQ results for the WHC area from 2002 to 2017 based on the normalized ELI (green and red represent the excellent and poor ecological livability quality levels, respectively).First, in terms of the seasonal variation, the ELQ in spring and autumn is distinctly better than that in summer and winter, with spring being the best and winter being the worst, that is to say, the results of the ELI in summer and winter are mainly occupied by yellow and orange areas, and for the winter results shown in Figure 6, the red and orange areas dominate (fourth row of Figure 6).These findings may be attributed to the weight results for each season: while the weights of the NDVI and LST are relatively high, the AOD in winter is also high and even exceeds the weight of the LST.However, it is difficult for the NDVI to dominate, even with a high weight, due to the low vegetation coverage in winter.This finding also indicates that the results of the ELI in winter are mainly affected by the AOD and NWD (both of their weights exceed 0.2).
The results of the ELI in the other seasons are mainly based on the weights of the NDVI and LST.Even if the temperature is higher in the spring and summer, the high vegetation cover offsets the corresponding influence.Furthermore, although the weight of the NDVI is lower than that of the LST in autumn, the relative suitability of the temperature adequately reflects the habitability.
In addition, the quantitative ratios of each quality grade based on the ELI were calculated.For brevity, the ELI was classified into three grades (poor, medium, and good), which refer to the 5-leveled RSEI map, according to its value range (Xu et al. 2018).Among them, ELI values of less than 0.4 and greater than 0.6 belong to the poor grades and good grades, respectively, and the medium grade for the ELI ranges from 0.4 to 0.6.The statistical result for Figure 6 is shown in Figure 7.In spring, the ELI from 2002 to 2017 in Figure 7 (a) was mainly approximately 0.6, indicating that the overall ecological livability quality level of the WHC area was acceptable.Specifically, the good grade of the ELI in the spring of 2002,2007,2012, and 2017 accounted for 85.419%, 49.750%, 42.461%, and 86.774%, respectively.This change indicates that the ELQ of the WHC area showed a deteriorating trend from 2002 to 2012, and then improved from 2012 to 2017.The ELQ in 2017 then recovered to basically the same condition as that of 2002.
The ELI of the WHC area in the other three seasons showed a similar pattern to that of spring, with all three seasons being above the medium grade, as shown in Figure 7 (b) to Figure 7 (d).The percentage of the good grade in summer decreased from 56.025% to 47.805% and then to 39.004% but increased to 66.711% from 2002 to 2017.The corresponding percentage in autumn decreased from 66.128% to 60.500% and then to 56.596% but then increased to 78.821%.The good grade in winter occupied a very low ratio, and the percentage above the general grade changed from 80.886% to 90.669%, but it decreased significantly to 65.561%, and finally rose to 88.695% in the next year.

Spatial distribution of the ecological livability quality in the WHC area
The results of spatial autocorrelation (Figure 8) and SDE (Figure 9) are utilized for the spatial characteristics analysis of the ELI for a better  formulation of urban ecological planning.The local Moran's I values for each year and season were approximately 0.9, indicating that the ELQ in the WHC area has a high local autocorrelation.In combination with Moran's I scatter map, the spatial pattern of the ELQ in the WHC area was visualized and its distribution characteristics were established when the significance level was set to 0.05.The local indicators of spatial association were employed to detect the spatiotemporal clusters of the ELI: High-High (HH), Low-Low (LL), Low-High (LH), High-Low (HL), and not significant (Figure 8).In general, the purpose of the spatial clusters was to identify the HH (hot-spot) or LL (cold-spot) patterns in the WHC area within a certain time period.In Figure 8, HH/LL means that the value of the ELQ of a spatial unit in the WHC area is high/low, and the surrounding spatial units are also high/low, reflecting the codirectional influence relationship.LH shows that the ELQ of a spatial unit is relatively low, but it is surrounded by nearby spatial units with highquality scores, while HL is just the opposite.
The first line in Figure 8 shows clustering results of the local indicators of spatial association in spring from 2002 to 2017.The HH regions are mainly distributed east of the WHC area, near East Lake, the southwest corner of Hongshan District, the southeast area of Hanyang District, and part of the riverside area.In contrast, the LL regions are mainly distributed in the center of the WHC area, including Jianghan District and Jiang'an District.Moreover, the LL regions show a tendency of spreading to the surrounding areas over time, especially from 2007 to 2017, and HH aggregation regions gradually appear in the central area.On the whole, the other three seasons show roughly the same pattern as the spring distribution.
Double SDE is applied in this paper.The second standard deviation range includes approximately 95% of the elements of the good grade ecological livability quality.Based on the shape of the SDE and the movement of the centers of gravity, the spatial concentration degree and the direction change trend of the good grade ELQ in the WHC area were analyzed.Figure 9 shows the SDE distribution and center of gravity positions of the four seasons from 2002 to 2017.In terms of the overall spatial distribution (Figure 9, left), the spatial pattern of the SDE in each season is roughly the same, and can be roughly divided into northeast and east by north, according to the direction of the long axis of the ellipse.This finding shows that the good grade ELQ in this direction is higher than that in the short axis direction and concentrated in the middle and east as the eastern part of Wuchang District and central part of Hongshan District have many lakes and forest parks, and the government has greatly invested in their protection and development.
Furthermore, the SDE direction generated each year is basically consistent with the Yangtze River, which means that the ecological and environmental impact in the WHC area is closely related to the Yangtze River.Although the SDEs in the winter of 2002 and 2012 were slightly off, the direction of the long axis, approximately east by north, was similar to the distribution direction of the Yangtze River and all the major lakes in the WHC area.Therefore, a close relationship between the ecological environmental quality of WHC and the distribution of open water sources, such as the Yangtze River in this region, may be observed.For the change in the center of gravity of the SDEs (Figure 9, right), the center of gravity of the good grade ELQ during the study period mainly fluctuated Wuchang District, and an overall trend of small amplitude fluctuation was observed in the east of the WHC area.Although the winter results for 2002 and 2012 show some deviation, the centers of gravity of these two SDEs are still east of the WHC area, which is consistent with the previous overall SDE distribution.Meanwhile, Figure 6 shows that the ELQ in winter is generally poor, and the area with an ELI greater than 0.6 is reduced.Therefore, the centers of gravity of the SDEs migrate to the areas with lakes and forest, i.e. east of the WHC area.

Discussions
The factors of human activities and the urban environment itself were discussed in relation to these spatiotemporal changes to ELI in the WHC area.The feasibility and universality of the ELI were investigated.

What accounts for the spatiotemporal trends of the ELI in the WHC area?
With the impact of human activities and climate change, the ecological environment and livability of urban areas are receiving more attention (Huang and Wang 2020;Wu 2020;Li et al. 2021).This paper shows the change trend and spatial difference of the ELQ in the WHC area in the last 15 years via the ELI.
In terms of all four seasons in each year, as shown in Figure 6, the ELQ in the northwest is the worst, and the worst area is located in Hankou, which is the economic and financial center of Wuhan.The main reason for these findings is the large number of compact buildings, a large population flow, and fewer water resources, such as lakes in these areas, leading to the prominent influence of the negative indicators in the ELI, such as the NDBSI, LST, and NWD.
On the other hand, from an annual perspective, the green areas in 2007 and 2012 were significantly smaller than those in 2002 and 2017.This phenomenon indicates that the ELQ of the WHC area first decreased and then increased.In particular, the results for winter in 2012 can be partly attributed to the extreme air pollution events in China in 2012, where haze blanketed eastern China for several days, including Wuhan (Tao et al. 2013).Meanwhile, during the Spring Festival, fireworks caused serious particulate matter pollution, and the air pollution index exceeded the standard.The AOD in 2012 was relatively high, contributing to the overall poor ecological livability quality (red and orange results of the ELI in Figure 6) in the prosperous areas (Jianghan District and the surrounding Qiaokou District).
Moreover, the ELQ of the WHC area showed a trend of degradation from 2002 to 2012, due to city expansion and development.Unlike this situation, under the policies of green development and ecological protection measures from the Twelfth Five-Year Plan of China, the good grade of the ELQ in the WHC area gradually increased from 2012 to 2017, and the ELI also increased accordingly.These results comprise response of human environmental protection actions and policy implementation to the ELQ of the WHC area.
Overall, an excellent ELQ was obtained for the vegetation and lake areas from 2002 to 2017, especially near East Lake.Due to the protection of the East Lake Scenic Area, particularly the opening of the East Lake Greenway in 2016, the ELQ near East Lake has gradually improved over the past 15 years.A poor ELQ was obtained around the second ring road of the WHC area, which was mainly employed for the construction of Wuhan's central business circle and residential buildings.For example, Qingshan District, which has many heavy industries, including chemical and manufacturing factories, often exhibited a poor level of ecological livability quality and showed a certain tendency for deterioration, due to the long-term impact of industrial pollution and the lack of lakes and vegetation (Trinder and Liu 2020).
In addition, there are many reasons for the phenomenon of the spatial distribution of the ELQ in Figure 8 from 2002 to 2017.On the one hand, the HH regions, such as those near East Lake, are well protected by green development and are less affected by urbanization, so the ELQ in these areas is high.On the other hand, many impermeable surfaces in the WHC area appeared as the city gradually expanded, leading to the outward diffusion of the LL regions with the rapid development of Wuhan.This situation is evident between the second ring road and third ring road in the WHC.However, after the implementation of a series of ecological environmental protection measures by Wuhan's government (Luo et al. 2019), the greening of the central area along the river was improved, and the ELQ in these areas was also improved.

The ELI can quickly and comprehensively evaluate the urban ecologically livable environment
Through the real-time feedback of the high-efficiency and fast ELI, it is possible to realize the management and analysis of dynamic remote sensing data and the real-time monitoring or identification of the urban ecologically livable environment from a multi-scale and high-precision perspective.Compared with LSTand-EVI-regulated-NTL-city index (LERNCI) for inner-city patterns (Liu et al. 2017c), the ELI pays more attention to the ecological livability, taking the ecological dimensions felt by citizens as the priority.Unlike neighborhood-level environmental performance in global cities from Urban Environment and Social Inclusion index (UESI) (Hsu et al. 2020), the ELI based on remote sensing data primarily reveals the quantitative characteristics of urban ecological livability in each landscape pixel grid.On the other hand, considering the comprehensiveness of the urban ecologically livable environment, the ELI constructed in this study adequately utilizes the five ecological indices of greenness, temperature, dryness, water-wetness, and atmospheric turbidity, to meet the evaluation needs of urban ecological livability quality.However, there are still some limitations to this study.First, in the process of remote sensing evaluation dataset generation or evaluation, it is difficult to avoid errors that may affect the accuracy of the characteristics or spatial distribution of the results in order to ensure the lower boundary of the model, for example, due to the limitation of meteorological station data, the overestimated value of temperature from the LST is selected (Venter Zander, Chakraborty, and Lee 2021).Second, there is a lack of investigation of other living indicators and the relationship between urban ecologically livable environment change in continuous years and potential driving factors, such as the social economy and policies.

Some future directions for the ELI
The ELI is primarily intended to show the intra-urban ecological livability or get its variability.In the future, some ecological indicators based on remote sensing analysis that are highly dependent on the spatial variability in urban environment need to be employed more rationally rather than overestimated, such as the biases between the LST and air temperature in different seasons.Some indicators, such as the planar PM 2.5 , may also be introduced in the ELI for comparative spatial variation.
In addition, research with long time series in continuous years that covers large areas could be carried out to explore the change rules of the urban ecologically livable environment in different cities, especially environmental events, through historical data and different image processing methods to provide a reference for urban ecological planning.More datasets combined with multisource remote sensing data, statistical data, and big data with higher spatial and temporal resolutions could be reasonably generated.Moreover, reasonable and objective comprehensive evaluation indices for the ecological environment should be established to provide decisionmaking support for ecological protection and environmental management.

Conclusions
pIn this study, the five urban ecological indices of greenness, temperature, dryness, water-wetness, and atmospheric turbidity were used to construct the ELI, which is a comprehensive evaluation index for remote sensing ecological livability.The ELI allows an efficient and real-time comprehensive assessment of the urban ecologically livable environments, and will be of great significance for the in-depth analysis of urban development.The advantages of remote sensings, such as large-area dynamic observation, timeliness, periodicity, and data comparability, have been considered by the ELI.Compared with site or statistical data and other evaluation frameworks, the ELI takes into account five aspects of urban ecology, and its results can be clearly visualized at national, urban, or even regional levels.Moreover, the spatial and temporal distribution pattern of urban ecological livability quality and its dynamic change characteristics can be analyzed by the ELI at different scales based on multi-source data and information processing technologies.
Specifically, the ELQ in the WHC area has been at a medium level (ELI ≈0.6), and its change showed a trend of stability or improvement from 2002 to 2017.From 2002 to 2012, the ELQ showed a stable or declining trend in different seasons but gradually turned to an improving trend from 2012 to 2017.In terms of the spatial distribution, the surrounding areas, such as East Lake and the riverside green parks, have a good grade of ecological livability quality.Meanwhile, Qingshan District, which has heavy and chemical industries, and the southeastern Hongshan District, which has limited water resources, have a poor grade of ecological livability quality.Although urban expansion led to ecological degradation and a poor grade of ELQ in the WHC area, urban afforestation and the implementation of government policies improved the urban environment.
In general, a remote sensing assessment index framework has been constructed to rapidly reflect the urban ecological livability quality with multi-source data, but there are still some to be further explored.More ecological factors could be considered to further optimize the evaluation model in future studies, such as air temperature, PM 2.5 , noise and other issues.In the future, more effective urban ecological evaluation systems and environmental protection policies should be implemented to provide decision-making support for ecological protection and environmental management, and a further reference for urban construction with harmonious development between humans and nature.

Figure 2 .
Figure 2. The Flow chart for establishing the ELI.

Figure 3 .
Figure 3. NWD results for the WHC area.
): where NDVI; LST; NDBSI; NWD; andAOD are the values of the five ecological indicators after normalization.
(a 1 ), Figure 5 (b 1 ), and Figure 5 (c 1 ) show the numerical differences of the three indicators, which are within the tolerable fluctuation range.The results in the right column of Figure 9 are the spatial distribution differences of the three indicators.As shown in Figure 5 (a 2 ), Figure 5 (b 2 ), and Figure 5 (c 2 ), the three fused indicators have a relatively high similar spatial distribution to the real distribution.Based on the quantitative results (Table

Figure 4 .
Figure 4. Real and fusion surface reflectance results for the near-infrared band of Landsat 8 for 15 September 2018.

Figure 5 .
Figure 5. Real and fused difference results for NDVI, LST, and NDBSI.

Figure 6 .
Figure 6.Temporal variation of the ELI in the WHC area.

Figure 7 .
Figure 7. Area percentage results for the ecological quality grades in the different seasons from 2002 to 2017.

Figure 8 .
Figure 8. LISA index clustering results for the ELI from 2002 to 2017 in every season.

Figure 9 .
Figure 9.Standard deviation ellipse distribution in every season from 2002 to 2017.

Table 1 .
Spatiotemporal fusion reference data applied in this study.

Table 4 .
Weight results for each index of the ELI in the WHC area.