Spatial and Temporal Variation of GPP and Its Response to Urban Environmental Changes in Beijing

Abstract

The carbon sequestration capacity of vegetation is the key to the carbon cycle in terrestrial ecosystems. It is significant to analyze the spatiotemporal variation and influencing factors of vegetation carbon sequestration ability to improve territorial carbon sink and optimize its spatial pattern. However, there is a lack of understanding of the impact of environmental conditions and human activity on the vegetation’s carbon sequestration ability, especially in highly urbanized areas. For example, effective vegetation management methods can enhance vegetation Gross Primary Productivity, while emissions of air pollutants like O₃, CO, NO₂, and PM2.5 can suppress it. This paper mainly explores the factors influencing vegetation carbon sequestration capacity across different regions of Beijing. Based on remote sensing data and site observation data, this paper analyzed the spatiotemporal variation trend of Annual Gross Primary Production (AGPP) and the influence of environmental factors and human activity factors on GPP in Beijing from 2000 to 2020 by using the Theil−Sen’s slope estimator, Mann−Kendall trend test, and comparing Geographically Weighted Regression method (GWR) and Geographically and Temporally Weighted Regression method (GTWR). GWR is a localized multiple regression technique used to estimate variable relationships that vary spatially. GTWR extends GWR by adding temporal analysis, enabling a comprehensive examination of spatiotemporal data variations. Besides, we used land use cover data to discuss the influence of land use cover change on AGPP. The results showed that the spatial distribution pattern of GPP in Beijing was higher in the northwest and lower in the southeast, and it showed an overall upward trend from 2000 to 2020, with an average annual growth rate of 14.39 g C·m⁻²·a⁻¹. From 2000 to 2020, excluding the core urban areas, the GPP of 95.8% of Beijing increased, and 10.6% of Beijing showed a trend of significant increase, concentrated in Mentougou, Changping, and Miyun. GPP decreased in 4.1% of the regions in Beijing and decreased significantly in 1.4% of the areas within the sixth ring. The areas where AGPP significantly decreased were concentrated in those where land use types were converted to Residential land (impervious land), while AGPP showed an upward trend in other areas. CO and NO₂ are the main driving forces of GPP change in Beijing. O₃ and land surface temperature (LST) also exert certain influences, while the impact of precipitation (PRE) is relatively minor. O₃ and CO have a positive impact on AGPP as a whole, while LST and NO₂ generally exhibit negative impacts. PRE has a positive impact in the central area of Beijing, while it has a negative impact in the peripheral areas. This study further discusses opinions on future urbanization and environmental management policies in Beijing, which will promote the carbon peak and carbon neutrality process of ecological space management in Beijing. Besides, this study was conducted at the urban scale rather than at ecological sites, encompassing a variety of factors that influence vegetation AGPP. Consequently, the results also offer fresh insights into the intricate nexus between human activities, pollutants, and the GPP of vegetation.

1. Introduction

The carbon sequestration capacity of vegetation is a crucial component of the carbon cycle in terrestrial ecosystems. Investigating the drivers of GPP spatiotemporal variations and their underlying mechanisms benefits revealing the main forces of GPP variations, thus promoting ecosystem management aiming at increasing carbon sequestration of terrestrial ecosystems [1,2,3]. Gross primary production (GPP) is the total amount of organic carbon fixed by plants per unit of time through photosynthesis, which reflects the vegetation productivity status of terrestrial ecosystems. It is one of the key indicators for assessing the carbon sequestration capacity of vegetation and understanding the carbon cycle of the terrestrial ecosystem [4]. Therefore, quantifying GPP and its spatiotemporal variation is of priority in carbon cycle studies [5]. Research indicates that the spatiotemporal variation in vegetation’s ability to capture carbon is influenced by multiple factors such as climatic conditions [6], vegetation cover [7], land use patterns [8], and soil properties [9].

The current studies on factors affecting GPP changes mainly focus on meteorological conditions and land use changes, etc. [10,11]. However, in urban areas, the impact of urbanization on GPP is complex and multifaceted. On the one hand, urban expansion leads to changes in regional land cover and a decrease in vegetation coverage, which may directly reduce the productivity of vegetation [12]. On the other hand, human activities brought about by urbanization also have indirect positive impacts [13], such as the urban heat island effect, increased air pollution in the atmosphere, and more effective vegetation management and planning [14], which collectively promote the photosynthesis and carbon fixation capacity of vegetation. These indirect effects can offset the direct negative effects in some cases and even promote the growth of vegetation productivity in certain areas [15]. It has also been found that the effect of air pollution on GPP cannot be ignored. For example, O₃, which has strong oxidizing properties, is absorbed by plants and damages their cell membrane tissues, thus reducing vegetation GPP [16,17], and aerosol particles can increase GPP by increasing scattered light [18], too much particulate matter can cause plant stomatal blockage instead of causing GPP reduction. For example, Yue & Unger [19], in their analysis of the effects of O₃ and aerosols on GPP during 2002–2011, found that O₃ led to a global GPP decrease of 0.91 ± 0.44 PgCyr-1 and aerosols led to a global GPP increase of 0.05 ± 0.30 PgCyr-1. Therefore, in this study, we used four types of air pollution, O₃, NO₂, CO, and PM2.5, as the influencing factors of human activities caused by urbanization. We used Land Surface Temperature (LST), Air Temperature (T), and Preparation (PRE) as environmental factors. With these factors considered, we can further our understanding of the responses of the carbon cycle to environmental change and human activities brought about by urbanization, especially for those fast-growing cities.

Beijing, as a highly urbanized city in China and a pilot area for sustainable and low-carbon development, has seen the majority of GPP research focus on ecologically sensitive areas such as the Mentougou district. This focus has led to a neglect of the spatial and temporal dynamics of GPP and the factors influencing it. When assessing the impact of environmental factors on GPP, it is imperative to consider the spatial dependence characteristic to prevent biased evaluation outcomes. The urban environment and its suburban counterpart interact through wind patterns, suggesting that GPP in suburban areas may be influenced by urban various factors [10,20,21]. The Geographically Weighted Regression (GWR) model proposed by Brunsdon et al. [22] is a local multiple regression technique used to estimate variable relationships that vary spatially and can effectively explore the impact of multiple factors on GPP. Furthermore, studies have indicated varying impacts of air pollutants and meteorological factors on vegetation growth at different stages of urbanization. The Geographically and Temporally Weighted Regression (GTWR) model introduced by Huang et al. [23] takes into account the temporal dimension to address spatiotemporal heterogeneity, which has considerable applicability in this problem. Although GWR and GTWR have been widely used in various studies, they are less frequently utilized in examining the determinants of GPP.

In this study, we used Theil Sen’s slope estimation, Mann Kendall (M-K) trend test, GWR, and GTWR methods to investigate the spatiotemporal changes of GPP in Beijing and explore the impact of environmental and human activity factors on GPP using AGPP data from 2000 to 2020 and meteorological and air pollutant data from 2014 to 2020. Additionally, we also employed land use/land cover (LULC) data to analyze the impact of LULC transitions on Annual Gross Primary Productivity (AGPP). Through this comprehensive approach, we aim to optimize the spatial structure and facilitate the city’s goal of achieving carbon neutrality.

2. Data and Methods

2.1. Study Area

Beijing, located at 39°56′ N latitude and 116°20′ E longitude, extends geographically from 115°25′ E to 117°30′ E and from 39°26′ N to 41°03′ N. The city experiences a temperate monsoon climate, with hot, wet summers and cold, dry winters. Surrounded by the Taihang and Yanshan mountain ranges, Beijing is positioned near the Bohai Sea to the east and the Yellow River to the south. This climate and location support a rich plant life and contribute to the city’s agricultural potential. Beijing, extending over 16,410.54 square kilometers and comprising 16 districts, boasts a varied landscape that includes both mountainous regions and plains. The mountainous areas, which constitute around 10,200 square kilometers or 62% of the city’s total area, and the plains, spanning 6200 square kilometers, together create Beijing’s rich geographical diversity. The city’s average elevation is 43.5 m; the plains typically range from 20 to 60 m in elevation, while the mountains peak over 1000 m. In the context of this study, areas above 500 m are classified as mountains, and those below are classified as plains. The plains of Beijing are mainly used for farming and urban development, whereas the mountainous areas are largely covered by forests and grasslands. Over the past twenty years, Beijing has seen a significant increase in urbanization, with the rate increasing from 77.52% in 2000 to 86.6% in 2020, and it is expected to reach 99.46% by 2030 [24].

2.2. Data Sets

In this paper, we choose the dataset of annual gross primary productivity in China’s terrestrial ecosystems during 2000–2020, published by Fan et al. [25]. This dataset was constructed based on the annual GPP observation dataset obtained from ChinaFLUX [26] long-term observation data and publicly available data, and an annual GPP(AGPP) assessment model constructed by simulating annual GPP per unit leaf area using a random forest regression tree model, thereby accurately depicting the geographical and temporal fluctuations in GPP across China. We utilized the CNLUCC dataset, a high-precision national-scale multi-temporal land use/land cover database derived from manual interpretation of Landsat imagery, as our source for surface coverage data. Additionally, we incorporated Land Surface Temperature (LST) data, monthly precipitation (PRE) and temperature records (T), and Chinese air pollutant data.

Table 1. Research data in this study.

Data	Spatial Resolution	Temporal Resolution	Period	Data Source
AGPP	500 m	Year	2000–2020	https://www.scidb.cn/en/detail?dataSetId=b496b208f51e44fcaf326e8b0f792c34, accessed on 9 February 2024
LST	1 km	8-day	2014–2020	https://ladsweb.modaps.eosdis.nasa.gov/, accessed on 9 February 2024
LUCC	30 m	5-year	2000–2020	https://www.resdc.cn/DOI/doi.aspx?DOIid=54, accessed on 7 March 2024
Air Temperature	1 km	Month	2014–2020	http://loess.geodata.cn, accessed on 11 February 2024
Precipitation	1 km	Month	2014–2020	http://loess.geodata.cn, accessed on 11 February 2024
PM2.5	1 km	Year	2014–2020	https://weijing-rs.github.io/product.html, accessed on 11 February 2024
CO	10 km	Year	2014–2020	https://weijing-rs.github.io/product.html, accessed on 11 February 2024
NO2	10 km	Year	2014–2020	https://weijing-rs.github.io/product.html, accessed on 11 February 2024
O3	10 km	Year	2014–2020	https://weijing-rs.github.io/product.html, accessed on 11 February 2024

delineates the specifics of these data sources, including their spatial and temporal resolutions. Due to the restricted period for which some pollutant data are available, this paper only examines the effects of seven environmental factors on AGPP from 2014 to 2020.

The AGPP data and LST data underwent preliminary processing in the ArcGIS 10.2 software environment, through the process including band extraction, reprojection, splicing, cropping, invalid value removal, unit conversion, etc., to finally obtain the annual average data of GPP and surface temperature in Beijing. At the same time, annual data were obtained for temperature and precipitation data using the image element statistics tool, and annual air pollutant data were obtained for air pollutant data in Beijing using the mask extraction and raster calculator tools.

To resolve the discrepancies in data resolution within our study, we implemented a systematic approach to standardize the data. Initially, we utilized ArcGIS to resample the data for the seven influential factors to a 500-m resolution, aligning it with the resolution of the AGPP data. Next, we determined the dataset’s minimum spatial resolution, finding it to be 10 km, and adopted this as our sampling unit. For each unit, we identified the pixel location with the median AGPP value and extracted the corresponding data for the seven influencing factors from that pixel. These values were then used as the representative dataset for each unit, enabling subsequent analysis and processing to be conducted at the unit level. This approach effectively addressed the issue of inconsistent data resolutions and minimized the potential distortions in GWR and GTWR analyses that can arise from resampling.

2.3. Methods

2.3.1. Theil−Sen’s Slope Estimator

When a time series exhibits a linear trend, the actual slope can be determined through a non-parametric method [27]. The slope estimation of the data is calculated using the following formula:

$$\beta =\mathrm{mean}\left({\displaystyle \frac{{\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_{\mathrm{i}}}{\mathrm{j}-\mathrm{i}}}\right),\forall \mathrm{j}>\mathrm{i}$$

(1)

where $x_j$ and $x_i$ represent data points at times  $j$ and $i$ respectively, with ($j$ > $i$). The median of the $N$ values of $Q_i$ constitutes a Sen’s slope estimator. A positive β signifies an ascending trend in the time series data, while a negative β points to a descending trend. This method, grounded in non-parametric statistics, is a robust approach for trend calculation. It is known for its computational efficiency, resilience to measurement errors, and resistance to outliers, making it a preferred choice for analyzing trends in extensive time series datasets.

2.3.2. Mann−Kendall Trend Test

The Mann-Kendall trend test evaluates time-series data for the presence of monotonic ascends or descends in the Y-axis values [28]. In this research, we employ the M-K test to identify such trends within the AGPP data series. The Mann–Kendall test for a time series is supposed to compute as follows:

$$\begin{array}{c}\mathrm{sgn}\left({\mathrm{x}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{j}}\right)=\left\{\begin{array}{cc}1,& {\mathrm{x}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{j}}>0\\ 0,& {\mathrm{x}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{j}}>0 \\ -1,& {\mathrm{x}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{j}}>0\end{array}\right.\end{array}$$

(2)

$$\mathrm{S}=\sum _{i=1}^{N-1} \sum _{j=i+1}^N \mathrm{sgn}\left(x_i-x_j\right)$$

(3)

where $N$ represents the dataset’s size, and $x_i$ and $x_j$ are the data values at times $i$ and $j$. A negative $S$ suggests a downward trend, while a positive $S$ suggests an upward trend. The trend’s statistical significance is evaluated using the variance of the MK statistic, applicable when $N>10$, and is computed as follows:

$$\mathrm{Var}\left(\mathrm{s}\right)={\displaystyle \frac{n\left(n-1\right)\left(2n+5\right)-{\sum }_{i=1}^p q_i\left(q_i-1\right)\left(2q_i+5\right)}{18}}$$

(4)

where $p$ indicates the number of concurrent groups within the dataset, and qi refers to the count of data points in the kth group. If the grouping variable is missing, the summation is excluded from the equation. The test statistic’s standard Z-value of the test statistic can be calculated using the subsequent formula:

$$Z=\left\{\begin{array}{l}{\displaystyle \frac{\mathrm{S}-1}{\sqrt{\mathrm{var}\left(\mathrm{S}\right)}}} S>0\\ 0;S=0\\ {\displaystyle \frac{\mathrm{S}-1}{\sqrt{\mathrm{var}\left(\mathrm{S}\right)}}}S<0\end{array}\right.$$

(5)

If the Z-value surpasses the critical value of significance level (α), the trend is identified as statistically significant.

This approach is advantageous as it is not contingent on data adhering to a normal distribution, nor does it necessitate linearity in trends. It remains unaffected by missing values and outliers, which has led to its extensive application in assessing trend significance within long-term time series. For this study, we utilize the “pyMann-Kendall” Python package developed by Hussain et al. [29] to conduct the M-K trend test.

2.3.3. Mann−Kendall (M-K) Mutation Test

The advantage of this method is that it is not only simple to calculate but also can determine the time when the mutation starts and point out the mutation area [30].

Assuming random independence in the time series, the statistics are calculated as below:

$$U_{Fk}={\displaystyle \frac{S_k-E\left(S_k\right)}{\sqrt{\mathrm{var}\left(S_k\right)}}}\left(k=\mathrm{1,2},\cdots ,n\right)$$

(6)

where $E\left(S_k\right)$ is the mean of $S_k$ and $\mathrm{var}\left(S_k\right)$ is the variance of $S_k$.

The above process is repeated according to the reverse order $x_n$, $x_{n-1}$,…, $x_1$ of the time series $x$, and the statistical variable ${UB}_k$ ($k$ $=n, n-1,..., 1$) is obtained. At the same time: ${UB}_k$ = −${UF}_k$. The ${UB}_k$ and ${UF}_k$ are plotted as curves. If the two curves intersect and the U value at the intersection is less than 1.96, that point is identified as the sequence’s mutation point, and the test confidence level is α = 0.05.

2.3.4. Geographically Weighted Regression Method

Geographically Weighted Regression enhances the traditional regression framework [22] by allowing the model’s parameters to vary spatially rather than remaining constant across the entire dataset. So that the model can be expressed as

$$Y_i={\beta }_0\left(u_i,v_i\right)+\sum _k {\beta }_k\left(u_i,v_i\right)X_{ik}+{\epsilon }_{i } i=1,\dots ,n$$

(7)

where $\left(u_i,v_i\right)$ signifies the spatial coordinates of point $i$, ${\beta }_0\left(u_i,v_i\right)$ is the intercept value, and ${\beta }_k\left(u_i,v_i\right)$ refers to the set of parameter values at point $i$. In contrast to the global model, which employs constant coefficients across different regions, this model permits spatial variation in parameter estimates, thus potentially enabling a more accurate representation of regional-specific impacts.

In order to calibrate the model, data near point $i$ are assumed to more significantly influence the estimation of parameters ${\beta }_k\left(u_i,v_i\right)$ than data that are more distant from the observation point $i$:

$$\widehat{\beta }\left(u_i,v_i\right)={\left[X^{T}W\left(u_i,v_i\right)X\right]}^{-1}{X}^{T}W\left(u_i,v_i\right)Y$$

(8)

where $W\left(u_i,v_i\right)$ is an $n\times n$ matrix with diagonal elements representing the geographic weights of the observation data for point $i$, while the off-diagonal elements are zero. This weight matrix is determined for every point $i$ where parameters are calculated.

The weight matrix in Geographically Weighted Regression signifies the varying significance of each data observation for parameter estimation at location i, with closer observations typically carrying more weight. Consequently, each point i’s estimate is associated with a distinct weight matrix.

$$W_{ij}=\exp \left(-{\displaystyle \frac{d_{ij}^2}{h^2}}\right)$$

(9)

$$d_{ij}=\sqrt{{\left(x_i-x_j\right)}^2+{\left(y_i-y_j\right)}^2}$$

(10)

where $h$ is a non-negative parameter called the bandwidth, which causes the influence to decay with distance. We utilized the Corrected Akaike information criterion [31] to identify its optimal value. $d_{ij}$ is the distance between locations $i$ and $j$, and ($x_i$, $y_i$) and ($x_j$, $y_j$) denote the coordinates of the respective points.

2.3.5. Geographically and Temporally Weighted Regression

The GTWR model offers a robust estimation approach by incorporating time and space data into its weight matrices, effectively capturing non-stationary spatial and temporal dynamics. This feature addresses the GWR model’s limitation of handling only cross-sectional data [23]. Thus, the model can be formulated as:

$$Y_{\mathrm{i}}={\beta }_0\left(u_{\mathrm{i}},v_{\mathrm{i}},t_{\mathrm{i}}\right)+\sum _{k=1}^d {\beta }_{\mathrm{k}}\left(u_{\mathrm{i}},v_{\mathrm{i}},t_{\mathrm{i}}\right)x_{\mathrm{i}\mathrm{k}}+{\epsilon }_{\mathrm{i}}$$

(11)

where $i$ ranging from 1 to n. $Y_{\mathrm{i}}$ denotes an n × 1 vector of explanatory variables, ${\beta }_0$ is the constant coefficient, and ($u_{\mathrm{i}},v_{\mathrm{i}},t_{\mathrm{i}}$) specifies the geographic and temporal coordinates of the observation. ${\beta }_{\mathrm{k}}\left(u_{\mathrm{i}},v_{\mathrm{i}},t_{\mathrm{i}}\right)$ is the unknown parameter for the kth factor at these coordinates, with x_ik representing the n × k matrix of explanatory variables. Parameters are estimated using a locally weighted least squares approach, where nearby observations for a given point are assigned higher weights, and the difference between observed and fitted values is used to assign lower weights. The method minimizes the sum of squared differences to derive the parameter estimates.

Constructing the GTWR model hinges on defining the spatial weight matrix. Typically, the space-time weight matrix is formulated as $W(u_{\mathrm{i}},v_{\mathrm{i}},t_{\mathrm{i}})\hspace{0.17em}$ = $diag$($w_{\mathrm{i}1}$,  $w_{\mathrm{i}2}$, …,  $w_{\mathrm{i}\mathrm{n}}$), with diagonal element $W_{\mathrm{i}\mathrm{j}}$ representing the spatiotemporal decay function. In this paper, according to the research of Huang et al. [23], a Gaussian function is utilized as the weighting function.

$$W_{\mathrm{i}\mathrm{j}}=\exp \left(-{\left({\displaystyle \frac{d_{ij}}{h}}\right)}^2\right)$$

(12)

where $h$ signifies the bandwidth that modulates the influence decay with distance, with its optimal value identified using the Corrected Akaike information criterion [31]. $d_{\mathrm{i}\mathrm{j}}$ denotes the space-time gap between points i and j. To normalize disparities across dimensions, the GTWR model automatically adjusts the spatial scale parameter $\lambda$ and the temporal scale parameter $\mu$. The spatial distance $d^S$ and the temporal distance $d^T$ are integrated to form a composite spatiotemporal distance $d^{ST}$, thereby enabling the formulation of a spatiotemporal distance function.

$$d_{\mathrm{i}\mathrm{j}}^{ST}=\sqrt{\lambda \left[{\left(u_{\mathrm{i}}-u_{\mathrm{j}}\right)}^2+{\left(v_{\mathrm{i}}-v_{\mathrm{j}}\right)}^2\right]+\mu {\left(t_{\mathrm{i}}-t_{\mathrm{j}}\right)}^2}$$

(13)

where $\lambda$ = 0 implies the absence of spatial effects, reducing the space-time distance to a function proportional to the time distance, thus converting the GTWR model into a basic time-weighted regression. Similarly, if $\mu$ is zero, the model disregards temporal effects, effectively making it a GWR model. However, when both λ and μ are non-zero, the GTWR model constructs a space-time weight matrix as ($u_{\mathrm{i}},v_{\mathrm{i}},t_{\mathrm{i}}$) = diag ($w_{\mathrm{i}1}$,  $w_{\mathrm{i}2}$, …,  $w_{\mathrm{i}\mathrm{n}}$), and $W_{\mathrm{i}\mathrm{j}}$ can be calculated as below:

$$W_{\mathrm{i}\mathrm{j}}=\exp \left\{-\left({\displaystyle \frac{\lambda \left[{\left(u_i-u_j\right)}^2+{\left(v_i-v_j\right)}^2\right]+\mu {\left(t_i-t_j\right)}^2}{h^2}}\right)\right\}$$

(14)

The GTWR model applies local regression to each observed unit, estimating parameters at different time points as geographical positions change. This approach better captures the spatial dependence and spatiotemporal differences of each driving factor. Compared with GWR, GTWR accounts for temporal non-stationarity and the integrated effects of space and time, rendering it apt for analyzing complex spatiotemporal data that requires consideration of both spatial and temporal factors. Due to the short time frame of our data, ranging from 2014 to 2020, and constraints in data availability, we will conduct GWR and GTWR analyses. This will allow us to compare the outcomes and choose the most appropriate method for our study.

In this study, we function the GTWR ADDIN by Huang et al. [23] in ArcGIS 10.8 to get and compare the results of GWR and GTWR models.

3. Results

3.1. Spatiotemporal Distribution and Variation of AGPP in Beijing

This study used AGPP data from SCIDB to analyze the spatiotemporal distribution and variation of AGPP in Beijing. This AGPP data contains some no-data areas and some invalid data areas that are less than zero. The result is shown in Figure 1. The geographic distribution of AGPP across Beijing reveals a discernible trend, marked by higher values in the northwestern areas and lower values in the southeastern areas. From 2000 to 2020, AGPP in Beijing ranges from below 0 to 3000 g C·m⁻²·a⁻¹, and the multi-year annual mean GPP varies from −116.18 to 2527.71 g C·m⁻²·a⁻¹ (Figure 1). The mean AGPP in mountainous areas was significantly larger than the mean AGPP in the plains. In the northwestern mountainous areas, such as the Jundu Mountains, AGPP generally exceeded 1000 g C·m⁻²·a⁻¹, while in the plain areas in the southeast, AGPP was generally lower, mostly below 1000 g C·m⁻²·a⁻¹.

Figure 2 shows spatial and temporal variation trends of AGPP, AGPP increased in most regions in Beijing, with a growth rate ranging from −55.35 and 71.98 g C·m⁻²·a⁻¹ (Figure 2a,b). Significantly decreasing areas are mainly concentrated within Beijing’s sixth ring road (Figure 2b). Figure 2c shows the 21-year time series of annual mean GPP. Figure 2c presents the 21-year time series of annual mean GPP, showing an overall increase from 2000 to 2020 at a rate of 17.41 g C·m⁻²·a⁻¹. The apex of AGPP values was registered in the year 2019. Conversely, the nadir was noted in 2003.

In the M-K curve of the AGPP mutation test in Figure 2d, the UF curve is greater than 0 except for the case of less than 0 in 2000–2004, indicating that the average annual Gross Primary Productivity has shown a steady upward trend since 2004, and the upward trend has been significant since 2011. The UF and UB curves had only one crossing point in 2013, and this crossing point was outside the critical line with a critical value of 1.96, so it was hard to judge that 2013 was the time when the mutation began.

3.2. The Influence of Urban Environmental Factors on AGPP

To analyze the influence of urban environmental factors on AGPP, GWR and GTWR models were tested with the original data sets. Before fitting GWR and GTWR, a multicollinearity test is necessary to eliminate redundant indicators. This study used the Variance Inflation Factor (VIF) as an indicator and screened redundant factors with a threshold of 5 [32]. The results of the multicollinearity test are shown in

Table 2. Variable Inflation Factor (VIF) Analysis of Remaining Factors.

Factors	VIF
CO	1.416183
LST	2.419883
NO2	3.076362
O3	1.932922
PRE	1.646616

. We removed the redundant factor with the highest VIF value step by step until the VIF value of all factors is not greater than 5. The results of the multicollinearity test for the remaining factors are shown in Table 2.

According to the results in Table 2, we removed PM2.5 and T. Then we fitted GWR and GTWR, and the results are reported in

Table 3. Comparison results of GTWR and GWR models.

	GTWR	GWR
Year	2014–2020	2014	2015	2016	2017	2018	2019	2020
Bandwidth (°)	0.114996	0.121176	0.124996	0.16544	0.237059	0.124996	0.124996	0.314818
ResidualSquares	84,427,900	7,762,450	8,123,490	9,794,280	15,555,900	8,693,000	9,987,170	14,828,000
AICc	13,191.3	1926.78	1930.96	1904.79	1935.11	1938.39	1953.52	1914.67
R2	0.602178	0.727502	0.731239	0.655626	0.526772	0.674349	0.724807	0.397604
R2 Adjusted	0.599978	0.716514	0.720401	0.64174	0.507691	0.661218	0.713711	0.373314

. The detailed results of GWR model can be seen in the Supplementary Materials.

Compared to the GTWR fitting results, the GWR fitting results exhibit notable year-to-year variations, which significantly diminishes their representativeness and reliability for in-depth analysis. Although the AICc value associated with GWR is considerably smaller than that of GTWR, this discrepancy is primarily due to the distinct data complexities inherent to each model. GTWR, which incorporates higher data complexity and greater model sophistication, naturally results in a larger AICc value. Therefore, the AICc values produced by GWR and GTWR should not be directly compared as they reflect different levels of model intricacy. A smaller AICc value for GWR does not necessarily imply that its fitting results are superior to those of GTWR. Consequently, the GTWR fitting results are deemed more appropriate and reliable for conducting a thorough analysis. Besides, we evaluated the significance of the GTWR fit using a p-value threshold of 0.05. The fitting results, along with their significance, are presented in Figure 3.

Figure 3 presents the GTWR fitting results derived from the original dataset, offering a comprehensive overview of how urban environmental factors influence AGPP in the Beijing area. The analysis reveals that CO and NO₂ are the most critical factors determining AGPP, with O₃ and LST also exerting certain influences, while the impact of precipitation (PRE) is relatively minor. Specifically, CO and O₃ generally have a positive effect on AGPP across most areas of Beijing, with their negative impacts being limited to a few regions in the northeastern and southwestern areas. In the northwest of the city, the positive effects of CO and O₃ are more pronounced. NO₂ and LST generally exhibit negative impacts on AGPP. LST only shows a positive effect on AGPP in a few areas of the southwest and east, with negative effects in the remaining regions and a more obvious negative impact in the southeast. NO₂ shows a positive effect on AGPP in a few areas of the northeast, southwest, and northwest, with negative effects in other regions and a more pronounced negative impact in the northwest. PRE has the least impact on AGPP, and its impact varies in different regions of Beijing. It has a positive effect on AGPP in most central areas of Beijing while having a negative effect in some northern and southwestern areas. Overall, the study results emphasize that CO and NO₂ play the most significant role in shaping the spatial pattern of AGPP in Beijing, with O₃ and LST also having certain influences. The impact of PRE is minimal and generally negligible. In terms of significance, the results indicate that the Geographically and Temporally Weighted Regression (GTWR) fits all factors with statistical significance across all domains, providing strong credibility for the model’s predictive power.

4. Discussion

4.1. Spatiotemporal Distribution of AGPP in Beijing

In this paper, we observed that the spatial distribution pattern of AGPP in Beijing is generally higher in the northwest and lower in the southeast, which is consistent with the findings of Gao et al. [33]. There may be two factors contributing to this phenomenon: elevation and human activities. On the one hand, the elevation may cause differences in vegetation types, material metabolism, functional structure, and other vegetation itself, and on the other hand, human activities in the mountainous areas in northwestern Beijing and the plain areas in southeastern Beijing are not equally active, and human activities and urban development in the plain areas are relatively more, which may lead to the shrinkage of vegetation area and the degradation of vegetation growth environment. For example, a study by Lü et al. [34] demonstrated that from 2000 to 2010, 979.74 km² of vegetation in Beijing experienced severe degradation, 725.90 km² experienced moderate degradation, and 301.65 km² was mildly degraded. These degradation phenomena are directly related to the increase of urban development land. The phenomenon of lower AGPP in the area around the urban center is also consistent with the findings of Li [35] and Liu [36], and this result is mostly related to a large number of human construction activities in urban development. In addition, this paper also found that most of the AGPP in key ecological protection areas in Beijing showed a significant upward trend, which suggests that the increase of AGPP in Beijing may be related to the implementation of national and urban policies of returning farmland to forest and closing mountains for forestry.

In terms of temporal variation, this paper found that the AGPP in Beijing exhibited a general upward trend from 2000 to 2020, which is consistent with the findings of Zhang Xinzhu [37] and Du Wenli et al. [38] Previous studies have shown that the combined effects of climate change, CO₂ fertilization effects, nitrogen deposition, and ecological restoration efforts promote an increase in AGPP, while factors such as drought events and urban development land expansion may cause a decrease in AGPP. In the case of Beijing, major events and ecological restoration work may be the main reasons for the fluctuating rise in AGPP from 2000 to 2020. The years 2004 and 2008 had the greatest ups and downs in AGPP in Beijing and were the two years when ecological work was vigorously promoted to host the Beijing Olympics under the Beijing Urban Master Plan (2004–2020). Under the guidance of the Beijing Urban Master Plan (2004–2020) and the government’s attention, the conditions for vegetation growth in Beijing were greatly improved, and AGPP surged.

4.2. Influence of LUCC on AGPP in Beijing

Land Use and Land Cover Change (LUCC) can exert an obvious influence on AGPP. For example, the conversion of vegetated areas into impervious surfaces can lead to a reduction in AGPP. Consequently, we also examined the effects of LUCC on AGPP using the CNLUCC dataset. To facilitate analysis, we have consolidated the original CNLUCC classification system into five primary categories: Cropland, Greenland, Water area, Residential land (impervious land), and Unused land. We have compiled land use/land cover data for Beijing spanning the years 2000 to 2020. The result is shown in Figure 4.

As depicted in Figure 4, it is evident that there has been a significant conversion of various kinds of land into impervious land within the central region of Beijing, whereas transformations among other kinds of land cover are relatively minor. Building on this observation, we conducted an analysis of land use alterations in Beijing over the period from 2000 to 2020 and compared the result with the trend of AGPP during the same period. The comparative outcomes are illustrated in Figure 5.

The result shows that areas with a negative trend in AGPP are predominantly found where the surface cover types have been converted to impervious surfaces, while AGPP remains stable or increases in other regions. The result indicates that urban sprawl exerts an obvious negative influence on AGPP. Therefore, for cities undergoing rapid expansion, the proper planning of urban growth areas and urban growth strategies is vital for the city’s green transformation and sustainable development.

4.3. Influence of Environmental Factors on AGPP in Beijing

In terms of environmental factors, LST has significant and varying impacts on AGPP in Beijing. LST exerts a negative influence in the northern, central, and southeastern regions, where high surface temperatures likely limit plant productivity. However, in certain areas, such as those extending eastward and westward from the central parts, LST’s impact shifts to a more positive influence, suggesting that moderate temperatures may benefit plant growth in these regions. This result highlights temperature as a key driving factor for vegetation productivity in Beijing. Within an optimal range, higher temperatures promote photosynthesis and boost plant productivity. However, excessive heat increases respiration and transpiration, causing dehydration and reducing productivity. These findings align with existing research [39,40,41].

The impact of precipitation (PRE) on AGPP in Beijing is relatively less significant, and its impact varies in different regions of Beijing. The negative impact of precipitation on AGPP in Beijing is attributed to the region’s temperate semi-humid continental monsoon climate, which is marked by summer precipitation and a dearth of precipitation in spring and winter. This seasonal precipitation pattern induces water stress in vegetation during early growth stages, consequently affecting AGPP [42,43]. Additionally, extreme weather events, such as rainstorms and droughts, can exacerbate the decline in vegetation productivity. Rainstorms can lead to soil erosion and nutrient depletion, while droughts impede plant water uptake and photosynthesis [43]. As shown in Figure 6, the PRE in the Beijing area fluctuates greatly, and the intense fluctuations in PRE have caused obvious disturbance to vegetation growth. In the central regions of Beijing, human activities are capable of regulating water resources. Therefore, PRE exerts a positive influence on the AGPP. On the periphery, however, human activities exhibit a limited ability to manage water resources effectively. Consequently, erratic precipitation patterns and the occurrence of extreme weather phenomena result in a negative impact of PRE on AGPP. In essence, the suppressive effect of precipitation on Beijing’s AGPP is likely due to the seasonal distribution of precipitation and the frequency of extreme weather events.

Therefore, in the context of environmental considerations, the utilization of appropriate water resource management strategies to align the water-heat relationship, enhance urban water resource regulation capabilities, and mitigate extreme events such as floods and droughts can significantly enhance vegetation productivity in Beijing and contribute to the advancement of sustainable urban development.

4.4. Influence of Pollutant Factors on AGPP in Beijing

In terms of pollutant factors, CO and O₃ generally have a positive effect on AGPP throughout most of Beijing, with O₃ having a less obvious effect. However, in small areas in the northern, southwestern, and southeastern regions, CO and O₃ exhibit negative effects. NO₂ generally exhibits negative impacts on AGPP while showing a positive effect on AGPP in a few areas of the northeast, southwest, and southeast.

CO can affect vegetation AGPP by affecting the antioxidant and drought resistance of plants [44]. The strong oxidizing properties of O₃ and NO₂ also damage plant structure and have an inhibitory effect on AGPP [16,17]. At certain concentrations, NO₂ can stimulate plant growth, but at high concentrations, it becomes toxic, adversely affecting plant development [45]. There is a certain gap between the model fitting results and existing research results; the reason could be as follows: O₃ and CO levels can also indicate the intensity of human activities in this study, and vegetation may benefit from positive human interventions, such as more effective vegetation management strategies, which may offset some of the negative impacts of pollutants. This could result in a positive influence of CO and O₃ on AGPP across most areas. Conversely, areas exhibiting negative impacts are likely due to severe local pollution of O₃ and CO, where the detrimental effects of pollution on AGPP surpass the positive effects of human activities. In terms of NO₂, unlike O₃ and CO, the main source of NO₂ in China is industrial emissions [46]. Consequently, the correlation between NO₂ and human activities is relatively low, and the results fitted by GTWR are closer to the natural effects of NO₂. As a result, NO₂ generally has a negative impact on AGPP, with only a small fraction of areas exhibiting a promotional effect on AGPP.

Therefore, employing more effective vegetation management practices can significantly mitigate the impact of pollutants on AGPP. Focusing on key polluted areas and implementing strategies to reduce pollutant emissions, coupled with appropriate vegetation management, can contribute to the enhancement of AGPP in the Beijing region and facilitate the achievement of sustainable urban development. At present, Beijing has implemented comprehensive initiatives in pollution control, yielding substantial outcomes. As depicted in Figure 7, while O₃ concentrations are anticipated to persist in their upward trajectory due to the rising number of motor vehicles, both CO and NO₂ exhibited a pronounced decreasing trend over the period from 2014 to 2020. This trend underscores the substantial success of Beijing’s environmental management strategies.

Moreover, Beijing has executed a multitude of strategies focused on the realms of urban greening and forest conservation in pursuit of sustainable urban development. For example, by implementing projects to enhance the ecological functions of urban green spaces, improve the quality and efficiency of plain forests, and accurately improve the quality of mountainous forests, the resilience of urban ecosystems has been strengthened. Besides, the Saihanba Mechanical Forest Farm’s desertification control has significantly raised forest coverage from 11.4% to 82%, offering critical ecological services to the Beijing Tianjin Hebei region. These measures and achievements reflect Beijing’s emphasis and commitment to the construction of ecological civilization.

5. Conclusions

In conclusion, this study highlights the significant spatial and temporal patterns of AGPP in Beijing, with higher values in the northwest and lower values in the southeast. Elevation and human activities, particularly urban development, play crucial roles in shaping these patterns. In terms of LUCC, the expansion of impervious land will have a significant negative impact on AGPP. Therefore, for cities experiencing rapid growth, the strategic planning of urban expansion is essential for facilitating a green transition and ensuring sustainable urban development.

Among the factors analyzed, CO and NO₂ are the most critical factors determining AGPP, with O₃ and LST also exerting certain influences, while the impact of precipitation (PRE) is relatively minor. Land Surface Temperature (LST) generally impacts AGPP negatively, whereas, in certain regions characterized by favorable temperatures, it enhances plant growth. The effect of PRE varies in different areas of Beijing, but it is not obvious. Irregular precipitation and extreme situations can both lead to negative impacts of precipitation on AGPP. Pollutants like CO and O₃ have a positive effect on AGPP in most areas, suggesting that effective vegetation management can offset their negative impacts. However, in severely polluted areas, the detrimental effects of these pollutants surpass the benefits of human activities. Unlike O₃ and CO, NO₂ exerts a detrimental effect on the AGPP due to its origins in industrial processes, which exhibit a low correlation with human activities. In conclusion, integrated water resource management, pollution control, and ecological restoration are crucial for enhancing AGPP and sustaining urban development in Beijing. These strategies are essential for fostering vegetation growth, achieving carbon neutrality, and ensuring environmental sustainability in the city. However, the impact of O₃ on AGPP in this study shows some differences from existing research findings. This discrepancy has not been fully examined and may be influenced by other environmental factors in Beijing. Further research and exploration are needed to understand this issue better.