Studying the Eco-Environmental Quality Variations of Jing-Jin-Ji Urban Agglomeration and Its Driving Factors in Different Ecosystem Service Regions From 2001 to 2015

Exploring the regional eco-environmental quality (EEQ) and its driving factors is of great significance for regional management. Although existing studies have paid much attention to evaluate EEQ, few studies have been performed to investigate the spatiotemporal variations of EEQ and its driving factors in different ecosystem service regions (ESR) at an urban agglomeration scale. In this study, we selected Jing-Jin-Ji urban agglomeration (JJJ) as the study area to evaluate its EEQ, analyze its spatiotemporal variations, and investigate potential driving factors explanatory power based on the geographical detector methods in different ESR during 2001 ~ 2015. The main conclusions were as follows: (1) The EEQ of JJJ had improved from 2001 to 2015, with the average RSEI increased from 0.43 to 0.46; among them, Bashang Plateau and Western Hebei Ecosystem Service Region (BWHE) had the highest RSEI change rate (+26.19%) and the highest NTEDI value (0.13), while Central Hebei Plain Ecosystem Service Region (CHPE) had the lowest RSEI change rate (−5.41%) and the lowest NTEDI value (−0.02). (2) The EEQ of JJJ had strong spatial agglomeration effects, with the global Moran’s $I$ increased from 0.82 to 0.88. Spatially, the LL regions mainly changed into the HH regions in the northwestern part, while in the central and eastern areas, some isolated LL regions displayed an aggregated trend. (3) In terms of the driving factors, soil type and elevation were primary factors in explaining the variations of EEQ. Specifically, natural factors explained the highest variations in BWHE. The interaction of topographical and socio-economic factors had high explanatory power in Yanshan and Taihang Mountain Ecosystem Service Region (YTME) and CHPE; To Bohai and Coastal Ring Ecosystem Service Region (BCRE), the interaction of meteorological and socio-economic factors accounted for the high variations of EEQ. All these findings could provide more valuable advice for relevant policy-makers.


I. INTRODUCTION
Since the implementation of the reform and opening-up policy in 1978, tremendous changes have taken place in mainland China, especially in the following aspects, urbanization, population expansion, and industrialization [1]. Nevertheless, along with these processes and intensive human activities, the regional eco-environment quality (EEQ) has also The associate editor coordinating the review of this manuscript and approving it for publication was Lefei Zhang . degraded to some extent, which had caused various environmental problems [2]- [4]. Furthermore, these problems had impeded the achievement of the goal of regional coordinated development. Hence, it is exceptionally urgent to evaluate regional EEQ.
Eco-environment defines as ''the total quantity and quality of water resources, land resources, biological resources and climate resources that affect human survival and development.'' It is a social-economic-natural compound system and an essential element for human subsistence [5], [6].
EEQ assessment is a useful tool to help decision-makers fully understand the natural and anthropogenic effects on the regional eco-environment and provide them with valuable information to formulate corresponding measures [7]. Till now, numerous researches have been conducted to evaluate the regional eco-environment at various scales, such as farmland [8], province [6], [9] prefecture-level city [10], urban agglomeration [5], river basin [7], [11], island [12], open pits [13], and tableland [14]. However, based on existing studies, it can be found that few studies have performed at the urban agglomeration scale. An urban agglomeration is a highly developed spatial form of integrated cities [15]. In 2014, the China government issued a development roadmap of the ''National New Urbanization Plan,'' which pointed out to establish a coordinated mechanism for urban agglomerations development [5]. To date, China has proposed building a hierarchical urban agglomeration system with five national-level large urban agglomerations, nine regional-level mediumsized urban agglomerations, and six sub-regional-level smallsized urban agglomerations [15]. Overall, as the final urban spatial form, evaluating the EEQ at the urban agglomeration scale will provide support to relevant policy-makers.
A comprehensive evaluation of the regional EEQ should include two parts, which are the evaluation of the local EEQ and the exploration of its driving factors. Firstly, regional EEQ evaluation mainly divided into two types, which are qualitative and quantitative evaluation methods. Qualitative evaluation places more emphasis on the ''state''; however, it can not provide the degree of these ''state.'' For instance, region A and region B are both in a good state, but we cannot know which region is better? Compared with qualitative evaluation, quantitative evaluation can give objective evaluation value of the eco-environment. To date, numerous quantitative evaluation methods have been proposed, such as the ecological footprint method [16], artificial neural network method [17], matter element analysis method [18], and so forth. However, these traditional methods mostly evaluate the whole regional eco-environment with one quantized value, which fails to acquire the result at any location. Under such background, integrating with remote sensing datasets and geographical information system technology to evaluate the regional eco-environment has attracted the attention of scholars all over the world [19]. Currently, remote sensing with its advantages of rapid, real-time, and large-scale monitoring has made much progress in land use classification [20]- [23], ecoenvironment evaluation [12]. For example, the normalized difference vegetation index (NDVI) has been widely used to monitor the regional vegetation situation [24]; the land surface temperature (LST) retrieved from thermal images is applied to evaluate the urban heat island [25]. However, the eco-environment is one complicated system; one single index can hardly evaluate it comprehensively [12], [19]. Therefore, it is common practice to construct one index system to evaluate the regional eco-environment [19]. He et al. developed a comprehensive evaluation index (CEI) integrating fine particulate matter (PM2.5) concentration, land surface temperature (LST), and vegetation cover (VC) to evaluate the urban environmental change [26]. Chang et al. established one index system, including fourteen indicators, to assess the eco-environment of the upper Hanjiang River [27]. Musse et al. integrated remote sensing and census data to assess Cali city's environment based on one proposed urban environment quality index (UEQI) [28]. He et al. (2019) combined remote sensing and statistical data to explore the regional differences of ecosystem health based on VORS (vigor, organization, resilience, and ecosystem service) framework [29]. Liou et al. synthesized 12 variables to assess the eco-environment vulnerability [30].
To sum up, existing researches mostly more or less selected some indicators which are cumbersome or difficult to obtain. Besides, some studies included too much subjectivity in the EEQ assessment, especially in the indicator weight setting aspect [12]. To make up for these shortcomings, one comprehensive aggregated index, remote sensing ecological index (RSEI) proposed by Xu [31], has made progress. This index integrated four aspects (greenness, heat, wetness, and dryness), which could reflect the intuitive eco-environment [12]. What is more, the weight of each indicator was set objectively based on the spatial principal component analysis method (SPCA). Till now, various studies have utilized RSEI to assess regional EEQ successfully [8], [32], [33]. Due to the comprehensive evaluation of the regional eco-environment, RSEI was selected to evaluate the EEQ of an urban agglomeration.
Secondly, explorations of the driving factors that affect the local eco-environmental changes have drawn much attention from scholars globally. Currently, numerous methods have performed, such as regression analysis [34], correlation analysis [35], and principal component analysis [36]. However, these traditional methods can only acquire the influence degree of individual factors and fail to explore the interactions of different driving factors and their combining effects on regional eco-environment variation [29], [37]. Compared with these methods, one promising approach, named the Geodetector method, proposed by Prof. Wang provides a solution to make up for these existing shortcomings [38], [39]. Geodetector method can analyze the single and combined effects of driving factors by investigating the spatial heterogeneities of geographical phenomena and the changes of various driving factors [38], [39]. Till now, several studies have performed by applying this innovative method to investigate the driving factors in the ecosystem health field or the urban thermal environment [29], [40]. However, in the ecoenvironment evaluation field, existing studies seldom explore the driving factors of regional eco-environment at different years and different ESR. Spatial heterogeneity widely exists in large regions. Hence, it is of considerable significance to investigate the regional EEQ, and its driving factors influence degree from the perspective of different ESR.
In this study, Jing-Jin-Ji urban agglomeration (JJJ) was selected as the study area to evaluate the EEQ and explore their driving factors from the whole and ESR views.
The objectives of this study are to a) understand the spatiotemporal distribution and variation of JJJ's ecoenvironment; b) explore the driving factors affecting the JJJ's eco-environment at the whole and ESR level. Our work will provide useful suggestions for relevant policy-makers to better formulate eco-environmental protection measures. The rest of this paper is organized as follows. Section II introduces the materials and methods, including the study area, study framework, RSEI construction, NTDEI construction, spatial autocorrelation, and geographical detector method. Section III details our results, which includes the spatiotemporal changes of EEQ and NTDEI in JJJ and different ESR, the spatial autocorrelation analysis, the factor, and interaction detection analysis of JJJ and different ESR. Discussions are presented in Section IV to compare the result of RSEI and EI and analyze the implications, limitations, and further study. Finally, Section V provides the conclusions and suggestions for relevant policy-makers.  (Fig. 1). This region has a temperate semi-humid and semi-arid monsoon climate with the average temperature in July and annual precipitation of 18 ∼ 27 • C and 524.4 mm, respectively. Based on JJJ's special geographical features, in this study, it was divided into four ESR, which were Bashang Plateau and Western Hebei Ecosystem Service Region (BWHE), Yanshan and Taihang Mountain Ecosystem Service Region (YTME), Central Hebei Plain Ecosystem Service Region (CHPE), and Bohai and Coastal Ring Ecosystem Service Region (BCRE) [41]. Nowadays, JJJ is considered as ''the three engines of China's economic growth'' in the 21st century. In 2017, the total population and gross domestic product reached 95.74 million and 8058.04 billion yuan, which accounts for 6.89% and 9.77% of that of the whole country (http://www.stats.gov.cn). Overall, with the fast urbanization, it is a hot topic to evaluate JJJ's eco-environment and explore their driving factors in JJJ and each ESR.

B. DATA SOURCES AND PROCESSING
In our study, the data can be divided into five categories: remote sensing data, statistical data, meteorological data, topographic data, and socio-economic data. The remote sensing data included two datasets, namely MOD09A1 and MOD11A2, which were provided by the NASA LP DAAC at the USGS EROS Center (https://lpdaac.usgs.gov) with a spatial resolution of 500 m and 1000 m. Statistical data included industrial Sulphur dioxide emission (SF) and second industry proportion (SI), which were acquired from the China City Statistical Yearbook (http://www.stats.gov.cn). Meteorological data, involving the annual average temperature (AT) and annual precipitation (PR), was collected from the China Meteorological Data Service Center (https://data.cma.cn). Topographic data included the elevation (EV) and slope (SP) at the raster format with a spatial resolution of 90 m, among it, the elevation data was downloaded from SRTM 90 m DEM Digital Elevation Database (http://srtm.csi.cgiar.org). The slope data was derived from the elevation data. Socio-economic data included the gross regional product (GP) and population density (PD) with a spatial resolution of 1 km, providing by the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (http://www.resdc.cn). Besides, this website can also give the soil type (ST) and land use (LU) at the raster format with a resolution of 1 km. Furthermore, MOD09A1 and MOD11A2 datasets were collected and processed at the Google Earth Engine platform (https://code.earthengine.google.com). All other datasets were preprocessed to at the 1 km spatial resolution and the same projection with MOD09A1 and MOD11A2 by using the ArcGIS 10.6 software.

C. STUDY FRAMEWORK
In this study, we proposed a framework to evaluate the JJJ's EEQ during 2001 ∼ 2015. This framework consists of three parts, which were the assessment of JJJ's EEQ based on RSEI, the spatiotemporal variations and spatial aggregation characteristics of JJJ's EEQ based on NTEDI and spatial autocorrelation method, and the investigation of driving factors explanatory power in different ESR. Figure 2 displays the detailed information.

D. RSEI CONSTRUCTION
The regional EEQ assessment was represented by the comprehensive remote sensing ecological index (RSEI). It was composed of four sub-indexes, which were normalized difference vegetation index (NDVI), wetness (WET), normalized difference build-up and soil index (NDBSI), and land surface temperature (LST). Existing studies found that this index could be influenced by the crop phenology [42], and was mainly applied at the city level. To overcome these problems, we integrated it with MODIS datasets in the period of 1st June to 31st October. Specifically, LST was directly derived from the MOD11A2 datasets, the left three sub-indexes were calculated based on Eq. (1)- (3).
where PC1 is the first component result; RSEI origin_i represents the RSEI origin value at the ith pixel; RSEI min is the minimal value of RSEI origin value; RSEI max represents the maximum value of RSEI origin .

E. NTEDI CONSTRUCTION AND INTRODUCTION OF SPATIAL AUTOCORRELATION
In this study, to understand the change situation of JJJ's ecoenvironment, one new index, named normalized time ecoenvironmental difference index (NTEDI), was proposed to analyze the regional eco-environment variations in a different time. The value of this index belongs to −1 ∼ 1. The higher the value, the better the EEQ, and vice versa. If the value equals 0, it shows that the EEQ does not change. The formula is as follows.
NTEDI = EQ T end − EQ T start / EQ T end + EQ T start (6) where EQ denotes the EEQ, which was represented by RSEI, T start and T end represent the start and end time.
The concept of spatial autocorrelation describes the correlation between a variable at a specific position and its neighboring position. If one observed spatial value changes in the same direction with the surrounding values, it is called a positive spatial correlation. If these two values are inversely related, it is called a negative spatial relationship. If the observed values show randomness in spatial distribution, this indicates that the spatial correlation is not apparent, and it is called a random distribution. The spatial autocorrelation measurement method can be divided into global and local spatial autocorrelation. Global spatial autocorrelation describes the overall distribution of some phenomena and judges whether these phenomena have an aggregated characteristic in particular areas. However, it fails to precisely point out these areas [43], [44]. Local spatial autocorrelation calculates the local spatial aggregation and points out their locations [45]. Moran's I is the most commonly used indicator to detect the spatial autocorrelation situation [43]. Eq. (7) and Eq. (8) are the formulas of global Moran's I (I g ) and local Moran's I (I l ) VOLUME 8, 2020 respectively.
where x i and x j denote the variable value at pixel i and j; represents the average value; w ij is the spatial weight function; N is the number of point pairs. In this study, RSEI was first resampled to 1 km; then, the spatial weights matrix was generated. Specifically, the inverse distance was used to conceptualize the spatial relationship, the euclidean distance was adopted as the distance method, the exponent was set to 1, and the number of permutations was set to 999 [46]. Finally, I g and I l were acquired based on the spatial weights matrix.

F. GEOGRAPHICAL DETECTOR METHOD
The geographical detector method can measure the spatially stratified heterogeneity of the geographic variable Y and explore how factor X explains the spatial pattern of Y [38], [47]. Different from traditional methods, the geographical detector method can analyze both categorical and discretized continuous variables; besides, it can investigate factors interaction result. The correlation degree between geographic variable Y and factor X is measured by the q-statistic value, which is expressed as Eq. (9) [39].
where h indicates the number of geographical variable Y stratified by factor X (a categorical variable or discretized continuous variable); N and N h represent the number of total grids and grids in stratum h; σ 2 and σ 2 h are the variances of geographical variable Y at the whole region and stratum h. The value of q ranges from 0 to 1. A higher q value indicates that Y has a stronger spatially stratified heterogeneity, and factor X can explain 100 × q% of the spatial pattern of Y (http://www.geodetector.org). In this study, Y represents the RSEI value in 2001 and 2015; X represents each driving factor. Each factor X needs to be categorized using a user-defined method or a categorization algorithm. Here, we choose four discretized methods, which are equal interval (EI), natural breaks (NB), quantile (QU), and geometrical interval (GI). Besides, each method at five, six, and seven groups was performed to acquire the optimal q statistic value [48].
The interaction detector is one promising tool to identify the combined effect of any two different driving factors, which has five categories (Eq. (10)- (14)).
Single factor nonlinear weaken : Nonlinear − weaken : where q is the factor detector value, X 1 and X 2 were the any two driving factors.    The local Moran's I result displayed the spatial distribution maps in 2001 and 2015 ( Figure 5). Typically, it includes five spatial clustering models, which are high-high aggregation (HH), high-low aggregation (HL), low-high aggregation (LH), low-low aggregation (LL) and not significant (NS). Based on Figure 5, the EEQ of JJJ was mainly composed of HH and LL. In both years, the LL aggregation areas were mostly found in the northwestern, southern-central, and eastern parts. The HH aggregation areas were mainly located in the northern regions and partly in the western areas. During 2001 ∼ 2015, the area of the LL aggregation in the northwestern part decreased and was mainly turned into HH aggregation, indicating the improvement of the ecological environment in these regions. In the central and eastern parts of JJJ, however, many isolated LL regions were gradually aggregated. Figure 6(a) and Figure 6(b) display the NTEDI value and spatial distribution of different grades. Among them, natural breaks (Jenks) method was applied to divide NTEDI into five grades. We found that the region with higher NTEDI value, including high and second high grades, mainly distributed in the northwestern part, except for one located in the eastern coastal area. As for those regions with a lower NTEDI value, including second low and low grades, they were mostly distributed in the central Hebei plain. To the medium grade, it mainly distributed in the Yanshan and Taihang mountains area.

D. EEQ AND NTEDI ANALYSIS IN DIFFERENT ESR
We analyzed the RSEI and NTEDI values of each ESR (Table 1). Specifically, compared with 2001, in 2015, the average RSEI value in BWHE, YTME, and BCRE all increased; however, the change rates were different; the order of three ESRs' change rate was as follows: BWHE (26.19%) > YTME (5.88%) > BCRE (2.70%). CHPE was the only region with a negative change rate (-5.41%). The same order was also acquired in the average NTEDI value. Combined with Figure 6(b), we could further find that the high and second high NTEDI grades were mostly distributed in BWHE, while second low and low NTEDI grades were mainly located in CHPE. Besides, YTME was mostly covered with the medium NTEDI grade. As for BCRE, all five NTEDI grades distributed in this region. Table 2 shows the optimal q statistic value. From the perspective of the whole region, the second industry proportion had explained the highest variation of JJJ's eco-environment (average q value was 0.4255), followed by the slope (average q value was 0.4149), the industrial Sulphur dioxide emission (average q value was 0.3954), and annual average temperature (average q value was 0.3809). Besides, except for the annual precipitation's q value decreased, all left factors' q value increased. Among it, the q value of land use had increased with the highest degree (174.63%), followed by population density (142.41%), elevation (141.59%), and gross regional product (103.04%).

E. FACTOR DETECTION ANALYSIS FROM THE WHOLE AND DIFFERENT ESR VIEW
To each ESR, the explanatory power showed a difference. Specifically, to BWHE, soil type, annual precipitation, slope, and land use were the top four factors that occupied the high variations of this region's eco-environment, with the average q value of 0.4999, 0.3888, 0.2999, and 0.2941, respectively. From the view of the change rate of q value in a different year, the annual average temperature had the highest decreased degree (−59.19%), other factors with a decreased degree had land use, second industry proportion, industrial Sulphur dioxide emission, and soil type. To those left factors with an increased degree, elevation had the highest increased degree (75.23%), followed by gross regional product (63.89%) and slope (39.48%). To YEME, elevation, annual average temperature, and population density were the top three factors with the average q value of 0.4271, 0.3805, 0.3677, respectively. The explanatory power of land use displayed the highest decreased degree (−15.89%). Besides, the gross regional product and soil type had a low decreased degree (−6.98% and −0.85%). Moreover, the q value of elevation, annual precipitation, and the annual average temperature had a higher increased degree, which was 35.25%, 32.15%, 24.32%, respectively. As for CHPE, the annual average temperature had the highest average q value (0.3239), followed by the second industry proportion and industrial Sulphur dioxide emission (0.3062 and 0.3005). Based on the change rate of q value, except for the annual precipitation and gross regional product had a negative change rate (−49.71% and −3.67%). In those factors with a positive change rate, elevation had the highest increased degree (581.25%), followed by the slope (295.00%), land use (185.76%), and population density (150.97%). To BCRE, the annual average temperature, second industry proportion, industrial Sulphur dioxide emission, and gross regional product were the top four factors with the high average q values, which were 0.4839, 0.4776, 0.4521, and 0.4330, respectively. From the change rate point of view, five factors had a negative change rate, which was the elevation (−66.42%), industrial Sulphur dioxide emission (−44.68%), second industry proportion (−43.99%), annual average temperature (−34.29%), and gross regional product (−8.53%). In the rest factors with a positive change rate, the slope had the highest increased rate (186.36%), followed by population density (74.71%), the annual precipitation (34.66%), soil type (31.05%), and land use (28.88%).

F. INTERACTION DETECTION ANALYSIS FROM THE WHOLE AND DIFFERENT ESR VIEW
The interaction of second industry proportion and industrial Sulphur dioxide emission of CHPE in 2001 displayed an irrelevant relationship. All the left interaction detection results in both years revealed that any two driving factors' interaction could enhance the explanatory power, whether in the whole region or the ESR. Table 3 listed the top three interaction detection results. We found that there existed a difference in a different year. Besides, all results represented a bi-or nonlinear enhancement.
Specifically, from the perspective of the whole JJJ, in both years, soil type, elevation, annual average temperature, second industry proportion, and industrial Sulphur dioxide emission together explained the spatial variations of JJJ's eco-environment. In 2001, soil type was an important factor in interacting with the other three factors, while in 2015, elevation was displayed as one important factor. From the regional scale, different ESR showed a difference. To BWHE, soil type, elevation, annual average temperature, annual precipitation, and population density had high explanatory power. Among it, in 2001, soil type acted as an important factor; in 2015, elevation represented as a significant factor. As for YEME, elevation, second industry proportion, annual precipitation, industrial Sulphur dioxide emission, and annual average temperature were important factors in explaining the eco-environmental variations, among them, elevation displayed as one important factor in both years. To CHPE, the second industry proportion, elevation, annual precipitation, annual temperature, industrial Sulphur dioxide emission, and gross regional product were key factors. In 2001, the second industry proportion was an important factor, while in 2015, elevation represented as a significant factor.  As for BCRE, annual average temperature, population density, gross regional product, soil type, second industry proportion, annual precipitation together explained the variations of regional EEQ. Among them, in 2001, the annual average temperature was found as a key factor. In 2015, the population density was displayed as an important factor. Besides, in 2001, the order of the average interaction value of top three interaction results was as follows: BCRE (0.7176) > YTME (0.6318) > BWHE (0.6139) > Whole (0.5264) > CHPE (0.3595). To 2015, the sequence was YEME (0.7031) > Whole (0.6884) > BWHE (0.5980) > BCRE (0.5566) > CHPE (0.4979). Moreover, in some regions, one factor had a low ability to explain the eco-environmental variations while interacted with other factors; it represented a remarkable increase. For instance, in BWHE, in 2015, the q value of elevation was 0.0985. However, the explanatory power increased sharply to 0.6003 when interacted with annual precipitation. Furthermore, in 2001, to CHPE, the q values of elevation and second industry proportion were 0.0256 and 0.2557, respectively. However, the interaction of two factors displayed a nonlinear enhancement. These indicated that driving factors from one aspect might act as a bridge enhancing driving factors from other aspects [29].

A. COMPARISON OF RSEI AND EI
In this study, we adopted the ecological index (EI), which was promoted by the Ministry of Ecology and Environment of China in 2015, to compare its results with RSEI [49]. Table 4 shows the average values of the two indexes at the whole and city level. Among it, RSEI value has rescaled to 0 ∼ 100 to compare EI at the dimension. From the perspective of the whole JJJ, two indexes all displayed an upward trend, further validating that the JJJ's EEQ had improved. From the view of the city level, eight cities showed the same change trend in both EI and RSEI, while five cities showed a different one. The reason for this might be the difference in the data acquisition and processing. The RSEI was derived from remote sensing datasets; therefore, it could reflect the JJJ's EEQ anywhere. However, to EI, it required both remote sensing and statistical data. Previous studies had concluded that the RSEI could reflect regional eco-environment changes more effectively than EI; even there existed a small difference [8], [19]. Hence, it was feasible to evaluate regional EEQ using RSEI.
Overall, based on both indexes, the eco-environment of JJJ had improved during 2001 ∼ 2015. Existing researches also reported similar conclusions. For instance, Wang et al. analyzed the land-use dynamics of JJJ. They found that deforestation, grass destruction, and grain cultivation gradually declined [50]. Indeed, the implementation of numerous projects by the government has contributed a lot [51]. For example, ''Returning Farmland to Forest (grass) Project,'' ''Three-North Shelter Forest Program,'' and ''Beijing-Hebei Ecological Water Resources Protection Forest Project,'' which were begun in 2000, 1979, and 2009, respectively [5]. Taking Zhangjiakou city as an example, the high NTEDI grade was mainly distributed in this city, combined with statistical data, at the end of 2015, farmland with almost 0.2257 million hectares had converted into the forest [50], [52]. To sum up, existing projects had made progress; however, more stringent plans should be formulated to acquire a better eco-environment in the future [50].

B. IMPLICATIONS OF DRIVING FACTORS ANALYSIS IN WHOLE AND DIFFERENT ESR OF JJJ
Exploring the driving factors at the whole and ESR scale will provide more detailed information in regional ecoenvironment management. In this study, the explanatory power of all driving factors in the different year have a significant difference.
From the whole perspective, the second industry proportion had a high q value, which represented its strong ability to explain the spatial variations of EEQ. Besides, its interaction with soil type and elevation enhanced explanatory power. The second industry proportion was one component of the industrial structure, while the latter had an important effect on the eco-environment [53]- [55]. Previous studies found that the increase of the second industry proportion exerted the most significant influence on the discharge of industrial solid wastes [55], it had been validated that industrial development in some developing countries could cause the deterioration of environment [56], [57]. Except for the second industry proportion, the high q value of industrial Sulphur dioxide emission also indicated that the air quality could also influence the regional EEQ [57]. The driving factors of land use and population density displayed a sharp increase during 2001 ∼ 2015, showing that anthropogenic activities exerted an increasing effect on the regional EEQ. Similar studies had proven that land-use changes and population expansion could have an impact on the variations of eco-environment [58], [59]. According to Table 3, we also found that the interaction of soil type with meteorological or topographical factors could increase the explanatory power. Namely, the variations of EEQ were not driven by one or two stable indicators, but a combined effect of numerous factors.
Compared with the detection results from the whole view, to each ESR, the q value showed a difference. To BWHE, soil type, annual precipitation, and slope had a high q value; namely, BWHE was mainly dominated by natural factors. Interaction detection revealed that two natural factors could enhance the explanatory power, for instance, the combined influence of soil type and elevation could explain 63.12% and 61.91% variations of BWHE's eco-environment in 2001 and 2015. This region belonged to the semi-arid grassland zone, and the elevation was in the range of 1300 m∼1600 m. Besides, this region was the primary area of soil formation and protection [41]. In the semi-arid climate zone, precipitation could influence the patterns of vegetation [60], [61], while vegetation was a crucial element in improving the regional eco-environment [60]. To YTME, elevation, annual average temperature, and population density had a high q value. The interaction detector revealed that the combined effect of elevation and second industry proportion could explain 63.29% and 71.75% variations of regional eco-environment. The primary land-use type in YTME was forest and grassland, with high vegetation coverage and high biodiversity. It had been reported that temperature exerted an effect on the vegetation [62], [63]. Besides, YTME located in the mountainous area; that is why the elevation was shown as an essential factor. CHPE, mainly covered by cultivated land, was not only an important food supply area but also a densely populated area in JJJ. The above analysis revealed that the EEQ of this region displayed a negative trend (−5.41%). Different from BWHE, this region was deeply driven by natural and human factors, including annual average temperature, the second industry proportion, and industrial Sulphur dioxide emission. Besides, the q value of elevation, slope, land use, and population density had increased sharply. Recently, the land-use intensity of this region was increasing, especially the continuous increase of construction land and the decrease of farmland. As for BCRE, even the annual average temperature had the highest q value; however, socio-economic and environmental factors all had a higher explanatory power. Generally, this region was close to the Bohai Sea and an area with rapid economic development [41]. Set the Binhai new district as an example, Hu et al. found that imperviousness and socio-economic factors had a high explanatory power of this region's thermal environment [40].
Overall, based on the factor and interaction detection results, it could provide more detailed information for regional eco-environment management. For BWHE, more measures should be taken to protect and improve soil formation. For example, change the mode of production of animal husbandry, ban grazing, and rest grazing, carry out house feeding and captive breeding, determine livestock by grass, and strictly control livestock carrying capacity. Besides, several existing projects should be strengthened to increase the regional vegetation coverage. What is more, the industrial structure should adjust to realize the virtuous cycle of its eco-environment. As for YTME, natural vegetation should be strictly protected; besides, overgrazing, disorderly mining, deforestation, and grassland reclamation shall be prohibited. To CHPE, we should develop ecological agriculture, gradually reduce the number of chemicals used, protect soil quality, and increase construction land moderately. To BCRE, it was suggested that the proportion of construction land should be adjusted when developing the economy. Besides, soil protection measures also need to be formulated [41], [64]- [67].

C. LIMITATIONS AND FURTHER STUDY
The spatial resolution of JJJ's eco-environmental quality was essential for the final analysis. In this study, the spatial resolution of RSEI at 1 km was acquired to perform the latter analysis. However, there still existed defects, as the higher resolution could provide more detailed information in conveying the eco-environmental variations. Nowadays, with the fast development of image pansharpening technologies, it can be expected that continuous remote sensing datasets at a high temporal and spatial resolution will be of great help to assess regional eco-environment at large scale. Besides, several conditions should be taken into consideration when utilizing RSEI. Firstly, considering the wetness information was derived from tasseled cap transformation formula, regions with large proportions of water bodies would increase the contribution of water bodies and could not reflect the true wetness information of earth surface; hence, it was not recommended to evaluate ocean's EEQ [42]; Secondly, as mentioned in Section II, crop phenology would exert an influence on the variations of regional EEQ; therefore, it was important to keep the time period at the same dimension and summer was suggested [19], [42]; Thirdly, RSEI was applied to evaluate the regional general EEQ, which was not suitable to some special cases, such as habitat analysis [19].
In addition, even the geographical detector method provided a way to detect the interaction results of two factors; it could not acquire the combined explanatory power of three or more factors. It would be further work to develop a novel method to analyze combined explanations of multi factors. Besides, ten driving factors of five aspects were detected; however, the regional eco-environment was a complex compound system; more factors need to be taken into consideration.
Generally, even there exist several limitations, we still believe that our study is meaningful. Although the spatial resolution of EEQ was 1 km, it still could reflect the spatiotemporal variations of the eco-environment. The detection of two factors interaction results has achieved a significant leap compared with traditional statistical methods. Overall, our proposed framework can provide more detailed suggestions for relevant policy-makers to formulate more specific measures to improve the JJJ's eco-environment.

V. CONCLUSION
In this paper, we evaluated the regional EEQ, explored the spatial distribution characteristics, and investigated the driving factors at the JJJ urban agglomeration and ESR scale from 2001 to 2015. Our research extended the existing knowledge of eco-environmental assessment studies from the following aspects. Firstly, we integrated RSEI and MODIS datasets to evaluate regional EEQ at the urban agglomeration scale successfully. Secondly, we analyzed the spatiotemporal distribution and spatial distribution characteristics and proposed a new index (NTEDI) to investigate the variations of JJJ's ecoenvironment at the whole and ecosystem service perspective. Thirdly, we explored the driving factors' ability to explain regional eco-environmental variations and detected a single factor's effect and combined effect of two factors in a different year by applying the geographical detector method.
Generally, the main findings of this study were as follows. Firstly, the eco-environmental quality in JJJ had improved during 2001 ∼ 2015, with the RSEI value increased from 0.43 to 0.46. Secondly, there existed strong spatial agglomeration effects of JJJ's eco-environmental quality, with the global Moran's I increased from 0.82 to 0.88. Local Moran's I revealed that, during 2001 ∼ 2015, in the northwestern regions, the LL regions showed a decreased trend and mainly turned into the HH regions, while some isolated LL regions in the central and eastern parts represented an aggregated trend. Thirdly, during 2001 ∼ 2015, the EEQ of BWHE had the highest RSEI change rate (+26.19%) and the highest NTEDI value (0.13), while the CHPE had the lowest RSEI change rate (−5.41%) and the lowest NTEDI value (-0.02). Finally, we found that driving factors displayed spatial heterogeneity. Specifically, soil type and elevation were primary factors in explaining the variations of eco-environmental quality. The interaction of natural factors (soil type, elevation, annual average temperature, and annual precipitation) explained the highest variations in BWHE. To YTME and CHPE, the interaction of topographical and socio-economic factors had high explanatory power. To BCRE, the interaction of meteorological and socio-economic factors (e.g., annual average temperature, second industry proportion, population density) accounted for the high eco-environmental quality variations. Overall, the assessment of regional EEQ and the investigation of driving factors in different ecosystem service regions could provide useful suggestions for relevant policy-makers.