Eco-Environmental Risk Assessment and Its Precaution Partitions Based on a Knowledge Graph: A Case Study of Shenzhen City, China

Abstract

The eco-environment is under constant pressure caused by the rapid pace of urbanization and changes in land use. Shenzhen is a typical “small-land-area, high-density” megalopolis facing various dilemmas and challenges; we must understand the eco-environmental risk (ER) of rapidly urbanizing regions and promote high-quality regional development. Therefore, with the help of the Python and Neo4j platforms, this study applies the theoretical foundation of knowledge graphs (KGs) and deep learning to form the KG of an ER; with this, we sort and establish an evaluation system in two dimensions, namely social and ecological, and introduce the Monte Carlo simulation to quantify the ER in Shenzhen City and its uncertainty from 2000 to 2020 to propose sub-regional programs and targeted measures for the prevention and control of the ER. The results are as follows: The eco-environmental risk index (ERI) of the study area as a whole showed a slight increase from 2000 to 2020; at the same time, the low-risk regions were mainly located in the east and southeast, while the high-risk regions were mainly located in the west–central and northwestern parts. In addition, three sample points (points A, B, and C) were selected using the Monte Carlo method to simulate the transfer of uncertainty from the indicator weights to the assessment results. Finally, based on the quantitative results, an accurate zoning scheme for ER prevention and control was provided to the decision makers, and appropriate countermeasures were proposed.

1. Introduction

Due to the rapid advancement of land use and urbanization on a global scale over the last hundred years, frequent and complex anthropogenic disturbances have extensively and profoundly influenced regional ecosystems [1] exposing the weaknesses in “zero-risk” environmental management and resulting in serious ecological risks in many regions [2]. Furthermore, the encroachment of construction and arable land on wetlands and waters exacerbates the reduction in ecological land. Irrational land change is the main driver of ecosystem service loss [3], resulting in a dramatic increase in global and local eco-environmental risks (ERs), which seriously threaten regional ecological security [4] and human well-being [5]. Therefore, how does one determine the ecological status of a region and identify its ER? What temporal and spatial ER patterns exist? What kind of insights can such patterns provide for the use and protection of regional land resources? The answers to the above questions are not only the key to promoting studies on the interaction mechanism between regional, high-quality development and the ecological environment, but also represent important research content for the pursuit of regional sustainable development in the future.

ER assessments originated in the 1990s, and since the 21st century, they have gradually attracted great attention from scholars at home and abroad; now, they have become a research hotspot in ecology, geography, and other disciplines [6]. The early ER studies mainly focused on the evaluation of single ecological factors, such as the soil [7], pesticides [8], the atmosphere [9], etc., and the results of these studies did not fully reflect the risk status of regional ecosystems. To solve this limitation, ER assessments based on land-use change came into being. ER assessment research focuses on exploring the occurrence [3], influencing factors [10], driving mechanisms [11], and ecological effects [12] of risks; identifying and analyzing ER issues requires the effective integration of academic research, policy formulation, and eco-environmental management, which can provide a scientific basis for socio-economic and environmental decision making [13]. At present, there are two main research models: One is the evaluation of ERs based on traditional risk sources, and the basic idea is to evaluate regional ERs by applying a relative risk model (RRM) using a source receptor, exposure, and hazard analyses for risk characterization [14]; The other one is based on the theory of landscape ecology, calculating the landscape pattern index and constructing an evaluation model to evaluate the ecological risk of landscapes; this enables a more unified research paradigm to form [15]. Based on the existing research results, the construction of an assessment framework and the determination of indicator weights are the keys to improving and enhancing the ER assessment results.

The construction of an assessment framework comes first. Among the many ER assessment frameworks, the risk source–risk receptor [16] and exposure–sensitivity–adaptive capacity [17] models are common, but the boundaries between the tiers are easily confused. Therefore, to improve the scientific validity of the composite index, we must urgently clarify the logical framework and organize and classify the assessment indicators, which will help to extract and summarize valuable information. With the booming development of knowledge graphs (KGs), their powerful information integration and knowledge expression abilities undoubtedly provide a new means to solve this difficult problem [18]. A KG describes concepts, entities, and their relationships in the objective world in a structured way to organize, manage, comprehend, and express information from different sources, resembling human cognition [19]. Currently, some well-known general KGs include DBpedia [20], Wikidata [21], Neo4j [22], etc., and the KGs in finance, healthcare, business, and other areas have corresponding industry applications [23,24]. Therefore, by combining the KG of ER with the help of KG theory, we can help to establish a unified framework for the assessment of ERs. Secondly, another one of the key components is the determination of the weights of indicators. Only by assigning scientific, objective, and reasonable weights to the assessment indicators can the research object be correctly assessed and real assessment results be obtained [25]. Different methods of determining the weights lead to large differences in the assessment results, and the uncertainty of each indicator or variable can be spread and accumulated [26], which may ultimately distort the ER assessment results. At present, the main methods of uncertainty analysis are the Monte Carlo simulation, Taylor simplification, sensitivity analysis, confidence interval methods, and others [27,28]. Among them, the Monte Carlo simulation method has been applied to risk assessment, but more researchers have simulated response coefficient values in the evaluation of relative risk models, verifying that the simulated risks of each risk unit coincide with the evaluated values obtained in this study, thereby quantifying their uncertainty [29]. Therefore, this study introduces a Monte Carlo simulation model to optimize the existing research methodology based on a KG-based assessment framework.

Overall, the existing studies have made great progress in theory and methodology, but there are still some shortcomings: (1) ER evaluation has non-uniform assessment frameworks; and (2) the different weight allocation methods lead to different evaluation results, which represent uncertainty in the final results. Therefore, this study develops a KG of the ER, which helps to establish a unified ER evaluation framework; at the same time, a mature Monte Carlo simulation model is introduced to facilitate indicator weighting in the ER evaluation, and the weights of each indicator are analyzed with the corresponding uncertainty, thus increasing the accuracy of the evaluation results. Moreover, as a key region in the Guangdong–Hong Kong–Macao Greater Bay Area, the ecology and economy of Shenzhen can be developed in a coordinated manner, so that it becomes a pioneering and exemplary location [30]. Therefore, this study overcomes the shortcomings of the existing studies and also strengthens the spatial quantitative study of the ER in Shenzhen, an economically developed region, which is important for solving regional development problems and even supporting national development. The main objectives are: (1) to grasp the temporal and spatial dynamics of the ER in Shenzhen, an economically developed region, from 2000 to 2020; and (2) to implement zoning management based on the results of an ER assessment to improve the ER issues and to further promote high-quality sustainable development.

2. Study Area and Data Source

2.1. Overview of the Study Area

Shenzhen City is located in the southern part of Guangdong Province, China, with a geographical location of 113°46′~114°37′ E, 22°27′~22°52′ N, and a topography of high in the southeast and low in the northwest, with an average altitude of 70–120 m (Figure 1). The total land area of the study area is 1953 km², and the total length of the coastline is 260.5 km, with natural ecological reserves such as the Dapeng Peninsula National Geopark and the mangrove forest in Shenzhen Bay.

2.2. Data Source and Processing

The data used in this study mainly include the administrative boundaries, the land-use data, the socio-economic statistics, and the composition of the geospatial data (

Table 1. Data categories and sources.

Category	Data	Data Sources	Website Source
Spatial data	China Basin and River Network Dataset	Resource and Environmental Science Data Center of the Chinese Academy of Sciences	http://www.resdc.cn/ (accessed on 2 June 2023).
	Vector boundary data
	GDP Spatial Distribution Kilometer Grid Dataset (Resolution: 1 km)
	Population Spatial Distribution Kilometer Grid Dataset (Resolution: 1 km)
	Digital elevation data (Resolution: 30 m)	Geospatial data cloud	http://www.gscloud.cn/ (accessed on 3 May 2023).
	Land-use data (Resolution: 30 m)
	Normalized Difference Vegetation Index (NDVI) (Resolution: 30 m)
	Landsat remote sensing images	USGS website	http://glovis.usgs.gov (accessed on 14 March 2023).
	Chinese soil data set in the Harmonized World Soil Database (HWSD)	National Tibetan Plateau Data Center	http://data.tpdc.ac.cn (accessed on 1 February 2023).
	Point of interest data	Baidu Map Company, Beijing, China	http://map.baidu.com (accessed on 15 February 2023).
Statistical data	Dataset of daily values of surface climate data in China (V3.0)	China Meteorological Data Network	http://data.cma.cn/ (accessed on 9 July 2023).
	Food production	China Economic and Social Big Data Research Platform	http://www.data.cnki.net (accessed on 14 April 2023).
	Shenzhen Water Resources Bulletin	Shenzhen Municipal Government Portal	http://www.sz.gov.cn/ (accessed on 14 April 2023).
	Shenzhen Statistical Yearbook	China Yearbook Network	http://www.yearbook.cn/ (accessed on 15 April 2023).

). Meanwhile, the data were processed using the cropping and resampling tools of ArcGIS 10.2 and uniformly converted into raster data with an accuracy of 1 km × 1 km to ensure the consistency of the spatial accuracy of the data. The sources of the spatial and statistical data used in this study are listed in Table 1.

3. Methods

This study is divided into three parts: firstly, the KG of ER was constructed, the indicators were extracted under the “socio-ecological system” framework, and the ER assessment index system was established; secondly, the Monte Carlo simulation model was introduced to simulate the weights and carry out the evaluation of the ER in the study area; and finally, the temporal and spatial distributions of the results of the ER assessment were characterized, and precaution partitions were carried out based on the overlay analysis results (Figure 2).

3.1. Building an Evaluation Framework Based on KGs

Based on the literature (CNKI and Web of Science databases) and web resources, with the help of the Python and Neo4j platforms, we apply the theoretical foundation of KGs and deep learning to represent the entities in the objective world and their relationships with each other through a network-like structure and graphs, which ultimately contributes to the formation of a KG used to avoid ERs. Secondly, constructing the KG of the ER involves schema layer construction, corpus establishment, knowledge extraction, and storage, etc. The specific operations are as follows:

Section 1: schema layer construction for the KG. The ultimate purpose of the KG of the ER constructed in this study is to serve as a risk assessment indicator system. Alongside this, from the perspective of the hierarchical framework, the KG schema layer of the ER is divided into six levels from top to bottom according to the framework of assessment theory: the object, property, methodological modeling (risk assessment process), status (risk source, impact factors, and risk issues), activity (risk governance activities), and relationship layers. By categorizing and sorting out the keywords and relationships in the ER domain, the KG schema layer of the ER is summarized (Figure 3).

Section 2: the semantic annotation for the KG. The data selected for this part of the study were taken from Acta Ecologica Sinica, Environmental Science, the Journal of Natural Resources, Transactions of the Chinese Society of Agricultural Engineering, Science of The Total Environment, Ecological Indicators, the Journal of Environmental Management, the Journal of Cleaner Production, and other specialized publications (all obtained from the journal databases or the respective official websites) and the Baidu Encyclopedia website. Currently, there is no publicly accessible semantic annotation set related to ERs, and many types of semantic annotation software do meet the users’ expectations [31]. The semantic annotation should start at the schema layer of the KG, which contains the entity and relation annotations. Relationships are expressed as triples (E1, R, and E2), and R characterizes the relationship between entities E1 and E2 [32].

In this study, both the literature (unstructured data) and other network data (semi-structured data) are converted into a text format, and the entities and relationships are annotated using a semantic annotation software tool developed by Baidu Brain (https://console.bce.baidu.com/easydata/app/, accessed on 7 July 2023). Among them, the semantic annotation process is as follows: (1) Annotate the entities. The annotators select the entity type after selecting the text, and then click the annotation button to complete the entity annotation; (2) Annotate the relationships. First, set the head and tail entities on the entity list, and then set the relationship type and click the annotate relationship button to annotate a relationship.

Section 3: knowledge extraction for the KG. To be able to clearly and intuitively see the characteristics of the dataset, the semantically tagged file data are further transformed into JSON format data. Furthermore, this study proposes a combination of BERT with a graph model structure and Global Pointer module. The model structure consists of three main parts: the BERT model, the graph module, and the Global Pointer joint-decoding module [33]. A detailed model is shown in Figure 4. Firstly, BERT is used to obtain dynamic word vectors rich in contextual information; then, the graph model is used to obtain local and non-local information about the word vectors; and finally, the Global Pointer decoder is used to convert the ternary extraction into a quintuple extraction, which solves the problem of ternary overlapping in relational extraction to a certain extent [34].

The BERT model employs a bi-directional transformer-encoding technique to capture the contextual features of the sequence, which enables Auto Encode, i.e., extracting valid data from the pre- and post-information of a sentence and pre-processing them efficiently [35]. In Figure 4, E₁, E₂, …, E_N are the embeddings of the bi-directional encoder representations from the transformers; each word in the sequence is obtained by token, segment, and position embedding. T₁, T₂, …, T_N are the targets of the bi-directional encoder representations from the transformers; they are the sequence vectors with rich semantic features obtained after feature extraction by the bi-directional transformers.

The Global Pointer model [36] is used for relational extraction, which is, in principle, a three-element set of SPO extraction, but in an experiment, it is essentially a “five-element group (sh, st, p, oh, and ot)” extraction, where sh and st are the head and tail positions of the subject entity, and oh and ot are those of the object entity, respectively. The model can be constructed using a probability map by simply designing a quintuple scoring function, and then enumerating and calculating all the quintuples, where the number of quintuples is calculated as follows:

$$m\times \frac{n(n+1)}{2}\times \frac{n(n+1)}{2}=\frac{1}{4}m{n}^2{(n+1)}^2$$

(1)

where the total number of p is m and the sentence length is n.

To accurately evaluate the performance strengths and weaknesses of the model [36], we adopt three basic evaluating criteria (consisting of precision (%), recall (%), and the F1 score (%)) for entity relationship extraction to evaluate the performance of the model. The formulas for each evaluation parameter are represented as Equations (2)–(4). In them, TP represents a correctly predicted positive sample, FP represents an incorrectly predicted positive sample, and FN represents an incorrectly predicted negative sample. The effectiveness of the BERT-based Graph Global Pointer model improves as the number of iterations increases (

Table 2. P, R, and F1 values of the model (Unit: %).

Times	$\mathit{P}\mathit{r}\mathit{e}\mathit{c}\mathit{i}\mathit{s}\mathit{i}\mathit{o}\mathit{n}$	$\mathit{R}\mathit{e}\mathit{c}\mathit{a}\mathit{l}\mathit{l}$	$\mathit{F}1(\mathit{F}1-\mathit{s}\mathit{c}\mathit{o}\mathit{r}\mathit{e}$)
Epoch = 50	0.82	0.81	0.81
Epoch = 100	0.85	0.86	0.85
Epoch = 200	0.86	0.86	0.87

).

$$Precision=\frac{TP}{tp+fp}\times 100\%$$

(2)

$$Recall=\frac{TP}{TP+FN}\times 100\%$$

(3)

$$F1-score=\frac{2Precision\times Recall}{Precision+Recall}\times 100\%$$

(4)

Section 4: knowledge storage of the KG. Neo4j, a graph database system that uses graph data structures for storage, significantly improves the data search performance and has become the primary tool for knowledge storage [37]. Therefore, in this study, the KG of the ER is stored in the graph database Neo4j. The main purpose of this is to make use of the LOAD CSV method in the Cypher language that is provided by the Neo4j database. Specifically, first, the parsed entity node and relationship data are saved as .csv files and placed in the import folder of Neo4j; then, the nodes and relationships are imported using the LOAD CSV statement of the Cypher language; finally, the relationships between the entities are stored in the Neo4j graph database with Cypher statements to form the KG of the ER.

3.2. Methodology for ER Assessment

This study constructs an ER assessment framework based on the KG of an ER and applies a Monte Carlo model to quantify the uncertainty that the weight determination process introduces to the assessment results by constructing a sample of weights. The specific process is as follows: ① The weights of the indicators were determined using the analytical hierarchy process (AHP) and the deviation maximization, entropy value, mean square deviation, and coefficient of variation methods to construct the weight samples; ② The triangular probability distribution of the weights was formed by the minimum, maximum, and mean values of the weight indicators; ③ Based on the weight samples and indicator normalization results, a Monte Carlo simulation was run using Python, and its mean value after 500 iterations was the ERI. It should be added that when five weight values are obtained and the amount of data is small and only the upper and lower bounds and the most probable values of the parameters are available, the probability density function of the triangular distribution [38] can be built from the data, and the Monte Carlo simulation can be performed. However, the formulas involved in the above process are as follows:

Standardization of indicators: The units for assessing the indicators of ER are diverse and have a certain impact on the results. Therefore, this study combines the attributes of the indicators with the dimensionless processing of the indicators [39]. The specific formula is as follows:

$$\mathrm{Positive}\mathrm{indicator}:Y_{ij}=\frac{X_{ij}-\mathit{min}(Xi)}{\mathit{max}(X_i)-\mathit{min}(X_i)}$$

(5)

$$\mathrm{Negative}\mathrm{indicator}:Y_{\mathrm{ij}}=\frac{\mathit{max}(X_i)-X_{ij}}{\mathit{max}(X_i)-\mathit{min}(X_i)}$$

(6)

where $X_{ij}$ stands for the actual value of the indicator in the i-th year to evaluate ER, $\mathit{min}(Xi)$ stands for the minimum value of each index to evaluate ER, $\mathit{max}(X_i)$ stands for the maximum value of each index to evaluate ER, and $Y_{ij}$ is the normalized value. After standardization, the values were between 0 and 1.

Quantifying the ER results: the indicators are normalized and their weight values are determined, and the subsystem risk index in the ER system of the study area and the composite index of the system are obtained by normalizing the values of the indicators and weighting them with the corresponding weights, using the following formula:

$${ERI}_i={\sum }_{j=1}^n{W}_j\times y_j$$

(7)

where ${ERI}_i$ stands for the composite quantitative index for year i; $W_j$ stands for the weight value of the j-th evaluation indicator; and $y_j$ stands for the standardized value of the j-th evaluation indicator.

3.3. Spatial Characterization Methods for ER

The Moran’s I index (I) and LISA index ($I_i$) were used to analyze the spatial correlation of ER in the study area. The Moran’s I can reflect the degree of similarity of the values of attributes of spatial proximity or similar units, with the formula as follows:

$$I=\frac{n{\sum }_i{\sum }_j{w}_{ij}\left(Y_i-\overline{Y}\right)\left(Y_j-\overline{Y}\right)}{\left({\sum }_{i\ne j}{w}_{ij}\right){\sum }_i{\left(Y_i-\overline{Y}\right)}^2}$$

(8)

The value of I is taken between [−1, 1], where $Y_i$, $Y_j$ stands for the value of the variable in the closest matching spatial units; $w_{ij}$ stands for the spatial weight matrix; $\overline{Y}$ stands for the mean value of the attribute values. If $I>0$, the observations of the study unit tend to be spatially clustered and have a positive spatial correlation; if $I<0$, the observations are spatially discrete and have a negative spatial correlation; and if $I=0$, the observations are spatially uncorrelated.

The LISA index reflects the degree of difference and significance between a region and its proximity zones, with the following formula:

$$I_i=\frac{Y_i-\overline{Y}}{S^2}{\sum }_{j\ne i}^{n^{\prime }}\left(Y_j-\overline{Y}\right)$$

(9)

$$S^2=\frac{1}{n^{\prime }}\sum {\left(Y_i-\overline{Y}\right)}^2$$

(10)

where $n^{\prime }$ stands for the sample size or number of study units and $S^2$ stands for the variance of the statistic. When $I_i>0$, it implies that a high value is surrounded by a high value or a low value is surrounded by a low value, and can be expressed as a high–high (H–H) agglomeration or a low–low (L–L) agglomeration. When $I_i<0$, it implies that high (low) regions are surrounded by low (high) regions, which can be expressed as a high–low (H–L) agglomeration or a low–high (L–H) agglomeration. It implies that the observed area is not associated with the adjacent regions and is indicated as “not significant”.

4. Results

4.1. ER-Based KG Used to Build an Assessment Framework

Through combing the KG of the ER, we found that there are big differences in the indicator system frameworks of different evaluation dimensions, and the theoretical approach to ER evaluation has gradually improved from determining the evaluation procedure framework to constructing various quantitative evaluation models. The current integrated evaluation takes into account the probability, intensity, and distributional heterogeneity of the spatial and temporal dimensions of the risk and source analyses. It has evolved from a small-scale evaluation of a single risk source and receptor to a large-scale comprehensive evaluation of multifarious sources of risk and receptors of risk. In this evaluation, we not only paid attention to external coercion, but we also considered the spatial structure, vulnerability, and other attributes of the risk receptor itself. Currently, other academics pay more attention to the anthropogenic sources of the risk and risk receptors that are heavily influenced by human activities, such as urban and rural ecosystems contaminated by heavy metals and persistent organic pollutants.

Furthermore, the main assessment dimensions that appeared more frequently in the KG of the ER were “potential–connectivity–resilience”, “risk sources–risk receptors”, “production–living–ecology”, “landscape pattern–ecological process”, “exposure–sensitivity–adaptive capacity”, and “pressure–state–response”. However, the above assessment dimensions tend to neglect the influence of social systems on ecological and environmental risks when sorting out the indicators. Furthermore, urban ecosystems, typically human–land systems, are the products of mutual trade-offs among multiple elements of socio-ecological systems; complex features, such as the diversification of risk sources and receptors, the determination of administrative boundaries, and the spatial heterogeneity of internal functional partitioning, make it impossible to follow the evaluation dimensions of general regional ecosystems in urban ecosystems (Figure 5). For the study area, which belongs to the Guangdong–Hong Kong–Macao Greater Bay Area, long-term resource exploitation and rapid urban development have changed the land-use composition, structure, and function, which have caused certain disturbances or responses in the region [40]. In response to these issues and in conjunction with the KG of the ER, two characteristic attributes of socio-ecological systems were used as the criterion layer for ER assessment and indicator screening (

Table 3. ER Assessment Indicator System.

Object	Subsystems	Indicators	No.	Indicator Character
ER	Ecological subsystem	Soil erodibility factors (K)	E1	+
		Surface temperature (LST)	E2	+
		Average annual rainfall	E3	+
		DEM	E4	+
		Carbon storage (Calculated by the InVSET model)	E5	−
		Soil conservation (Calculated by the InVSET model)	E6	−
		Water conservation (Calculated by the InVSET model)	E7	−
		Habitat quality (Calculated by the InVSET model)	E8	−
	Social subsystems	Population density	S1	+
		Disturbance density of transport facilities	S2	+
		Industrial density	S3	+
		Food production per unit area	S4	−
		Economic density	S5	−
		Density of parkland	S6	−

). Among them, the social subsystem mainly contains human activities that affect the ER (S1, S2, and S3) and responses to the ER during urban development (S4, S5, and S6), while the ecological subsystem contains natural background quality-threatening ERs (E1, E2, E3, and E4) and background quality responses to the ER (E5, E6, E7, and E8). Therefore, according to the KG of the ER, this study combed the indicators (Figures S1 and S2 in Section S1 of the Supplementary Materials) from the socio-ecological subsystem dimension and constructed an ER assessment framework (Table 3).

4.2. Spatial and Temporal Heterogeneity of ERs

4.2.1. Overall Changes

According to the triangular probability distribution of the weights (

Table 4. The triangular probability distribution of each indicator weight.

Indicator	Max	Min	Mean	Indicator	Max	Min	Mean
E1	0.064	0.015	0.037	S1	0.101	0.009	0.052
E2	0.102	0.023	0.062	S2	0.089	0.049	0.068
E3	0.095	0.014	0.063	S3	0.088	0.012	0.054
E4	0.084	0.015	0.055	S4	0.115	0.018	0.071
E5	0.106	0.019	0.064	S5	0.533	0.035	0.168
E6	0.131	0.017	0.072	S6	0.135	0.021	0.078
E7	0.156	0.014	0.074
E8	0.186	0.018	0.081

) and the ER evaluation model, the average ERI of each evaluation cell was measured, and the spatial coverage of the ER at the grid scale from 2000 to 2020 was obtained (Figure 6), which showed that the spatial and temporal differentiations were very obvious. Meanwhile, using the natural breakpoint method, the ER was classified into five levels, corresponding to these risk zones: the lowest-risk (ERI ≤ 0.35), lower-risk (0.35 < ERI ≤ 0.45), medium-risk (0.45 < ERI ≤ 0.55), higher-risk (0.55 < ERI ≤ 0.60), and highest-risk areas (ERI > 0.60).

On the time scale, the mean ERI in the study area from 2000 to 2020 was 0.495 (2000), 0.501 (2010), and 0.538 (2020). Moreover, the highest ERI showed an increasing trend (from 0.705 in 2000 to 0.751 in 2020), and the proportion of highest-risk areas gradually increased (7.79% in 2000; 13.99% in 2010; and 20.47% in 2020) (

Table 5. Area proportion of each level of ER.

Year	Lowest-Risk Area	Lower-Risk Area	Medium-Risk Area	Higher-Risk Area	Highest-Risk Area
2000	18.72%	17.08%	28.73%	27.68%	7.79%
2010	20.27%	17.88%	21.07%	26.79%	13.99%
2020	16.11%	20.46%	17.57%	25.39%	20.47%

). Furthermore, the 2000 and 2010 ERIs were dominated by medium- and higher-risk areas, while the 2020 ERI was dominated by higher- and highest-risk areas. Among them, the proportion of high-risk areas (higher- and highest-risk areas) increased by 10.39%, the proportion of low-risk areas (lowest- and lower-risk areas) changed less and remained around 35% overall, while the area proportion of medium-risk areas decreased by 11.16%. In addition, in terms of spatial distribution, the high-risk areas in 2000 were mainly concentrated in the coastal zone of the study area, covering several districts (the Bao’an, Nanshan, Futian, and Luohu Districts); in 2010, they were concentrated in the western part of the study area (the Bao’an, Nanshan, Guangming, and Longhua Districts); and in 2020, the distribution of areas in 2010 was similar, but the final distribution was much more concentrated. The continuous expansion of construction land area in these areas has generated an ecological pressure, leading to a continuous increase in the ER (Figure 6).

Throughout the study period, Moran’s I of the ERI was greater than 0, showing a positive spatial correlation, i.e., the ERs interacted with each other and were spatially similar (

Table 6. Moran’s I of the ERIs.

Value	2000	2010	2020
Moran’ s I	0.756	0.785	0.791

). Overall, Moran’s I of the ERI exhibited an increasing trend over the last 20 years, which indicated that the degrees of spatial aggregation and heterogeneity of ER increased due to socio-economic development and land-use change. The “H-H” agglomeration units (“hotspots”) account for 30.29% (2000), 34.27% (2010), and 32.19% (2020) of the total study area from 2000 to 2020, respectively. In 2000, the “L-L” aggregation unit was mainly concentrated in the southeastern part of the study area, while in 2010 and 2020, the “L-L” aggregation unit showed an expanding trend and was mainly concentrated in the southeastern and northwestern parts of the study area (Figure 5). Overall, the number of “H-H” increased, and then decreased, with an overall decreasing trend, while the number of “L-L” increased. The decrease in the number of non-significant spatial aggregation units indicates more local spatial aggregation of the ERs.

4.2.2. Analysis of Subsystem Changes

First, we analyzed the social subsystem. The distribution of the risk in the social subsystem for the period 2000–2020 gradually decreased, and the maximum, mean, and minimum values increased by varying degrees. The increase in the maximum and minimum values was smaller, with the maximum value increasing from 0.413 (2000) to 0.445 (2020) and the minimum value increasing from 0.101 (2000) to 0.129 (2020). Between 2000 and 2010, the spatial distribution of the risk in the social subsystem was slightly limited, with high values mainly found in the southeastern part of the study area and low values mainly found in the core development areas of the districts and counties; the spatial distribution of the risk in the social subsystem in 2020 was significantly different from that in the previous two periods, with gradually fewer low-value areas (more concentrated in the southwest of the study area) and gradually more high-value areas (Figure 7).

The overall mean indicator values in the social subsystem showed varying degrees of change over the period 2000–2020. For the positive indicators (S1, S2, and S3), the mean valued showed fluctuating changes. The mean values of S1 were 0.194 (2000), 0.196 (2010), and 0.192 (2020), while the mean values of S3 were 0.161 (2000), 0.218 (2010), and 0.148 (2020). On the other hand, the mean values of S2 increased, with mean values of 0.132 (2000), 0.138 (2010), and 0.165 (2020). For the negative indicators (S4, S5, and S6), the changes in the mean values of S5 and S6 were larger. The mean values for S5 were 0.466 (2000), 0.487 (2010), and 0.856 (2020); for S6, they were 0.884 (2000), 0.829 (2010), and 0.941 (2020). It should be added that with socio-economic and regional development, the GDP (S5) and the density of parklands (S6) have increased substantially, especially in the southwestern part of the study area, but the disparity between regional development also increased, and there are many low-value zones, which make the overall value after standardization high, and this is what leads to an increased risk in the social subsystems (Figure 7).

Second, we analyzed the ecological subsystem. The distribution of the risk in the ecological subsystem from 2000 to 2020 gradually widened, with the maximum, mean, and minimum values showing different degrees of change. Among them, the mean of the risk value of the ecological subsystem decreased a little, with the mean value decreasing from 0.280 (2000) to 0.278 (2020), whereas the maximum and minimum values showed some changes, with the maximum value increasing from 0.410 (2000) to 0.428 (2020) and the minimum value decreasing from 0.144 (2000) to 0.119 (2020). The spatial distribution of the risk in the ecological subsystem from 2000 to 2020 had a high degree of similarity, with the high-value areas mainly concentrated in the western, central, and northeastern parts of the study area, while the low-value areas were mainly concentrated in the southeastern part of the study area (Figure 7).

The overall mean values of the indicators in the ecological subsystem also showed different degrees of change from 2000 to 2020. For the positive indicators (E1, E2, E3, and E4), the changes in their mean values were small. The mean values of E2 decreased: 0.593 (2000), 0.527 (2010), and 0.504 (2020); the mean values of E3 decreased in a fluctuating manner: 0.472 (2000), 0.409 (2010), and 0.452 (2020); and the mean values of the other two indicators (E1 and E4) were stable. Therefore, to a certain extent, E2 and E3 contributed more to the assessment of the results. For the negative indicators (E5, E6, E7, and E8), the magnitude of change in their means was large. The change in the mean value of E6 was not significant, while the means of E5, E7, and E8 increased. The mean values of E5 were 0.345 (2000), 0.367 (2010), and 0.375 (2020); the mean values of E8 were 0.502 (2000), 0.545 (2010), and 0.565 (2020); and the mean values of E7 were 0.432 (2000), 0.452 (2010), and 0.481 (2020). Therefore, to some extent, E5, E7, and E8 contributed more to the assessed outcomes (Figure 7).

4.2.3. Uncertainty Analysis of ER Assessment Results in the Study Area

In the Monte Carlo simulation, the uncertainty of the weights was gradually transferred to the assessment results during calculation. As a result, there was uncertainty in this evaluation result. Moreover, the ERI was calculated using the ecological and social subsystem risk indices, and the uncertainty was directly transferred through them.

The areas with the highest ERI, as shown in the assessment results, were mainly located in the southwest and northwest of the study area (point A). This is because the area rapidly developed and has weak ecosystem services, which is one of the reasons for the high-level ERI. The ERI indices for point A ranged from 0.58~0.64 (2000) to 0.47~0.53 (2010) and 0.65~0.69 (2020), and their mean values for 500 iterations were 0.614 (2000), 0.519 (2010), and 0.675 (2020), respectively. The adj. R² values were 0.806 (2000), 0.805 (2010), and 0.907 (2020), respectively. In addition, the risk indices for the ecological and social subsystems of point A also showed uncertainty. The mean values of the 500 iterations of the risk index for the ecological subsystem of point A were 0.309 (2000) (range: 0.27–0.33), 0.312 (2010) (range: 0.27–0.33), and 0.359 (2020) (range: 0.32–0.38). Their adj. R² values were 0.816 (2000), 0.831 (2010), and 0.842 (2020), respectively. In contrast, the averages of 500 iterations of the risk index for the social subsystem at point A were 0.298 (2000) (range: 0.28–0.32), 0.225 (2010) (range: 0.19–0.24), and 0.312 (2020) (range: 0.31–0.32). Their adj. R² values were 0.821 (2000), 0.866 (2010), and 0.913 (2020) (Figure 8 and Figure S3 in Section S2 of the Supplementary Materials), respectively.

The assessment results show that the lower areas of the ERI were mainly located in the southeastern part of the study area (point B); this is because this area is used for ecological purposes, the landscape type is dominated by woodland, the ecosystem service function level is high, and human activities cause minimal impacts, which was also one of the reasons for the reduction in the ERI. The ERIs for point B were 0.27~0.34 (2000), 0.32~0.40 (2010), and 0.35~0.41 (2020); their means of 500 iterations were 0.314 (2000), 0.365 (2010), and 0.387 (2020); and the adj. R² values were 0.845 (2000), 0.747 (2010), and 0.751 (2020), respectively. Meanwhile, the mean values of 500 iterations of the risk index for the ecological subsystem at point B were 0.133 (2000) (range: 0.09–0.15), 0.165 (2010) (range: 0.14–0.19), and 0.175 (2020) (range: 0.14–0.20). The adj. R² values were 0.828 (2000), 0.842 (2010), and 0.846 (2020). The averages of 500 iterations of the risk index for the ecological subsystem at point B were 0.193 (2000) (range: 0.17–0.22), 0.198 (2010) (range: 0.17–0.22), and 0.213 (2020) (range: 0.18–0.23), respectively. The adj. R² values were 0.909 (2000), 0.865 (2010), and 0.917 (2020) (Figure 8 and Figure S3 in Section S2 of the Supplementary Materials), respectively.

The results of the assessment showed that the areas with a medium ERI were mainly located in the central part of the study area (point C). This is because this area is a peripheral area for urban development, where ecosystem service functioning and human activities were some of the causes affecting the ER. The ERIs for point C were 0.37~0.41 (2000), 0.44~0.50 (2010), and 0.55~0.60 (2020); their means of 500 iterations were 0.394 (2000), 0.474 (2010), and 0.581 (2020); and their adj. R² values were 0.930 (2000), 0.814 (2010), and 0.867 (2020), respectively. Meanwhile, the averages of 500 iterations of the risk index for the ecological subsystem at point C were 0.186 (2000) (range: 0.15–0.21), 0.287 (2010) (range: 0.24–0.30), and 0.305 (2020) (range: 0.27–0.33). The adj. R² values were 0.753 (2000), 0.839 (2010), and 0.817 (2020), respectively. In contrast, the averages of 500 iterations of the risk index for the social subsystem at point C were 0.211 (2000) (range: 0.201–0.253), 0.209 (2010) (range: 0.18–0.22), and 0.264 (2020) (range: 0.24–0.28). The adj. R² values were 0.923 (2000), 0.922 (2010), and 0.894 (2020) (Figure 8 and Figure S3 in Section S2 of the Supplementary Materials), respectively.

4.3. Precaution Partitions to Circumvent ER

ER evaluation is the process of collecting, collating, and expressing scientific information, which is conducive for environmental decision making [13]. In addition, in the context of the development of information technology, the promotion of accurate ER knowledge is not only an intrinsic requirement for ecological civilization construction, but it is necessary for the promotion of high-quality sustainable development in the region [41]. In this study, the rank of each ER in the study area from 2000 to 2020 was obtained by superimposing the distribution graphs. As shown in Figure 9, the top three areas transferred were: medium risk to higher risk (9.66% of the total area), higher risk to the highest risk (9.55% of the total area), and the lowest risk to lower risk (7.06% of the total area). It can be seen that most of the ER grades in the study area shifted from low to high grades, and the ER was intensified, which, in turn, led to serious ecological problems. In addition, according to the needs of ER prevention and control, the study area was divided into key, strict, and general control areas via the process of ER level transfer, as described below.

(1)

The key control area of the ER: This area is mainly where the ER shifted from a low to a high grade, and the ER deteriorated. The key control areas are mainly located in the western part of the study area, where ecological land use is scarce and scattered, and construction land use accounts for a large proportion and covers an ever-expanding area, resulting in ecological and environmental pressures; the ER continues to rise. In response to such problems, our future work will focus on improving land use, rationalizing industrial arrangements, and transforming the mode of economic development. Furthermore, establishing a proper ecological compensation system will be essential to promote the natural environment and sustainable socio-economic development [42].

(2)

The strict control zone of the ER: This zone is mainly where the grade of ER has not shifted, and the ER has not improved significantly. It should be noted that the strict control zone is targeted at areas with high ER levels (highest- and higher-risk areas). Strict control zones are mainly located in the southwest and southeast of the study area and should play a leading role in territorial spatial planning, where the incremental development of construction land is accompanied by the renewal of stock, to minimize the negative impacts on the ecology and maintain the current level of ecological risk, while promoting economic development [43]. In addition, it is essential to improve the capacity and quality of environmental facilities to make full use of their functions and potential applications, such as the development of ecotourism, the improvement of nature reserves, etc. This management model promotes both the ecological protection of the specific environment and the growth of the regional economy, creating a virtuous circle of benefits.

(3)

The general control area of the ER: This area is mainly where the ER shifted from higher to lower grades, with some improvement in the ER. It should be added that the general control area focuses more on the areas where an improvement in rank is not significant. The general control zone is mainly located in the central part of the study area, with a small amount in the southeastern and western parts of the study area. To mitigate some of the negative impacts associated with socio-economic development, targeted ecological protection and restoration should be carried out in this region. In addition, the ecological advantages of the site need to be maintained to prevent further deterioration of the ER, while strengthening the infrastructure [44].

Figure 9. Risk level transfer of the ER and precaution partitions.

Figure 9. Risk level transfer of the ER and precaution partitions.

5. Discussion

5.1. Exploration of the ER Assessment Framework

Firstly, with the rapid development of big data and artificial intelligence technologies, more and more research has been devoted to the development of knowledge-driven ER assessments, which monitor and explain the development of ERs and produce large amounts of academic data [45]. However, the existing research results [46] have tried to show the feasibility and effectiveness of the KG in terms of its application in the ecological field [32] in multiple tasks and referenced a typical research framework from construction to application. Therefore, by constructing the KG of the ER, this study comprehensively and deeply grasps the interrelationships and evolutionary processes between various elements in the field of ER and provides effective data and decision-making support for the establishment of an ER evaluation index system.

Through the KG of the ER, it was found that when the ER is assessed for ecological categories, such as wetlands [47], rivers [13], and lakes [48], the assessment framework is mainly constructed from the risk source–risk receptor perspective, focusing on the impact of the natural ecological environment on the ER [16]; when the ER is assessed for ecological categories, such as croplands [6] and cities [49], the assessment framework is mainly constructed using the dimensions of disturbance and exposure of the ER, taking into account the impact of natural ecology and human activities [17]. Although the above assessment frameworks can effectively characterize the ER, they often confuse the classification or attribution of the indicators. To address this problem, this study considers the urban system, as a typical complex land-use system, as a product of multiple trade-offs between the social and ecological systems [50], and the urbanization process exacerbates the conflicts between the ecological environment and socio-economic development. Therefore, this study introduces the socio-ecological system framework [51], which is conducive to elucidating the interactions between the social and ecological subsystems and revealing the non-linear, multi-level constitutive characteristics of the system.

Finally, human or economic activities on urban land pose a challenge to the ER in fast-growing cities [13]. These urban systems are highly vulnerable and face multiple threats [3,52]. The concentration of people, industries, and infrastructure in urban areas increases the risk of natural hazards, which highlights regional differences in ER resilience [45]. This study area is a highly urbanized region, which focuses more on the disturbance or response to ERs generated during urban development. Therefore, this study screens the dimensional indicators of ecosystem service function (E5–E8), climate (E3), human activities (S1), urban facility development (S2, S3, S5, and S6), etc. In conclusion, the ER assessment indicators constructed through the KG in this study overlap with the findings of Na et al. [41], Guo et al. [44], Zhao et al. [2], etc.

5.2. Exploration of ER Assessment Methods and Results

Establishing an ER assessment method that is applicable to cities has been a research hotspot in this field [53]. Compared with traditional ER assessments, the method in this study can better determine the weights of indicators. The key to the quantification of ERs is to construct reasonable and appropriate weights, and the different distributions of the weights of each indicator will directly affect the evaluation result [25]. In the existing research results, more methods, such as the entropy value approach [54], analytic hierarchy process (AHP) [55], fuzzy comprehensive evaluation (FCE) [17], etc., are used to determine the weights of the indicators. However, among the many weighting methods, how does one choose the appropriate method to truly reflect the level of ER? To address this problem, this study introduces a Monte Carlo simulation model to simulate the weights for ER evaluation in the study area, while the impact of the uncertainty of the weights on the evaluation results was analyzed. Compared with subjective weighting, objective weighting, and combined subjective–objective weighting methods [56], this methodology effectively simulates the transfer of uncertainty from the indicator weights to the quantitative results, and at the same time, can integrate ER uncertainty and spatial information, making the ERI quantitative results more objective and reliable.

Secondly, based on the quantitative results of the ER from 2000 to 2020, the overall risk level increased. Among them, the average value of the social subsystem risk index is higher, and its change is more significant, which is also the key factor leading to the increase in the risk level. The areas with high risk levels are more concentrated in the areas with more human activity or dispersed ecological lands, which affects the final quantification results. Furthermore, the areas with many ecosystem services are too fragile, so the social subsystems are at a higher risk. Therefore, the rational coordination of social and ecological subsystems is still an important step to avoid these risks. In addition, based on the results of the Monte Carlo simulation operations, we selected three sample sites (A–C) to analyze the ER levels and their uncertainties (95% confidence intervals) and found that the predominantly impervious region with a high level of urban development has high ERIs; the region with a high level of vegetation cover and fewer urban development activities has low ERIs; and the region between the former and the latter levels has medium ERIs. For the urban systems, the ERI is highly correlated with the quality of the ecosystem and the level of urban development [57]. In conclusion, the methodological system constructed shows a lower overall risk level in areas with ecological land use or less development and fewer construction activities than that in the core areas of urban construction and development, which is somewhat similar to the pattern described by Lu et al. [13] and Shen et al. [58].

5.3. Shortcomings and Prospects

This study examines the level of ER in Shenzhen over the last 20 years at the grid scale and develops improvement measures based on the zoning scheme. Compared with the existing research results, there have been some improvements in the judgment perspective and research methods. However, there are some shortcomings. First, although this study combines the social and ecological dimensions to construct the evaluation index system of the ER based on a KG, the index system can only alleviate the problem of inconsistency in the existence of evaluation indexes to a certain extent, and the results of the study are still unavoidably affected by the indexes. Second, this study quantifies the uncertainty of the ER assessment process, but lacks an in-depth exploration of the division of ER levels and the driving mechanism in the risk evolution process. Finally, this study mainly focuses on the evolution process of ER ranking when determining the zoning scheme, but the management measures still need to be improved.

Given these shortcomings, this study still needs to explore many aspects: for example, the evaluation index system of the ER and its processing method still deserve in-depth analysis and selection; to determine whether there is a more scientific and effective method for ER rating; and the natural breakpoint method of GIS also deserves to be investigated. Additionally, we must determine whether the spatial resolution of the other data can be further improved, and so on. This study will further improve the evaluation framework of ER, better reflect the interaction between the ecological and socio-economic factors, and formulate more integrated and comprehensive zoning programs and management measures. It is believed that as the research progresses, more in-depth evaluation mechanisms and methods will be gradually established and improved to provide more scientific guidance for the realization of regional high-quality development.

6. Conclusions

Based on the KG of the ER, this study constructs the indicator system of ER in the study area from two dimensions: social and ecological. Furthermore, the study quantifies the uncertainty due to the determination of weights with the help of the Monte Carlo simulation model, divides different levels of risk control areas according to the quantification results, and puts forward risk avoidance measures in a targeted manner. The conclusions are as follows:

(1)

The average values of ERI from 2000 to 2020 are 0.495 (2000), 0.501 (2010), and 0.538 (2020), and its highest value shows an increasing trend. In terms of spatial distribution, the high-risk area in 2000 was mainly concentrated in the coastal zone of the study area, and the areas distributed in 2010 and 2020 had some similarity, but the distribution results in 2020 were more concentrated. In addition, the mean value of the risk index of the social subsystem was higher than that of the ecological subsystem in the ER’s assessment system.

(2)

In order to better elucidate the uncertainty in the ER assessment, three sample sites (A–C) were selected to describe the means and intervals of change of their ERI, risk index of ecological subsystem, and risk index of social subsystem from 2000–2020. To some extent, it shows the uncertainty generated by different weights of indicators and its transmission process.

(3)

Based on the quantitative results of ER and its transfer characteristics in the study area from 2000 to 2020, it was divided into key control areas, strict control areas, and general control areas of ER, and targeted improvement measures were proposed. Although the study area is a high-speed urbanization area, particular attention is paid to guaranteeing the coordination of the ecological and social subsystems.