Estimation of Greenhouse Gas Emissions of Taxis and the Nonlinear Influence of Built Environment Considering Spatiotemporal Heterogeneity

Abstract

The fuel consumption and greenhouse gas (GHG) emission patterns of taxis are in accordance with the urban structure and daily travel footprints of residents. With taxi trajectory data from the intelligent transportation system in Xi’an, China, this study excludes trajectories from electric taxis to accurately estimate GHG emissions of taxis. A gradient boosting decision tree (GBDT) model is employed to examine the nonlinear influence of the built environment (BE) on the GHG emissions of taxis on weekdays and weekends in various urban areas. The research findings indicate that the GHG emissions of taxis within the research area exhibit peak levels during the time intervals of 7:00–9:00, 12:00–14:00, and 23:00–0:00, with notably higher emission factors on weekends than on weekdays. Moreover, a clear nonlinear association exists between BE elements and GHG emissions, with a distinct impact threshold. In the different urban areas, the factors that influence emissions exhibit spatial and temporal heterogeneity. Metro/bus/taxi stops density, residential density, and road network density are the most influential BE elements impacting GHG emissions. Road network density has both positive and negative influences on the GHG emissions in various urban areas. Increasing the road network density in subcentral urban areas and increasing the mixed degree of urban functions in newly developed urban centers to 1.85 or higher can help reduce GHG emissions. These findings provide valuable insights for reducing emissions in urban transportation and promoting sustainable urban development by adjusting urban functional areas.

1. Introduction

The transportation industry represents a major source of worldwide greenhouse gas (GHG) emissions and is the second-largest emission sector behind only the construction sector, constituting roughly 23% of GHG emissions [1]. Over the past two decades, the travel demand of urban residents and the proportion of motorized travel have continuously surged in China. Thus, vehicle tailpipe emissions have become a major source of GHGs in urban areas. For example, in Hong Kong, China, GHG emissions from the urban transportation system are the largest source of atmospheric pollution [2]. Due to the outbreak of COVID-19, global pandemic prevention and control policies led to a reduction in GHG emissions from urban transportation [3,4]. However, by January 2023, following the complete reopening of Chinese cities, urban activities had gradually returned to previous levels [5]. This revival in urban activity led to a sustained rebound in resident travel, resulting in further increases in traffic congestion and GHG emissions. In response to the issue of GHGs from urban vehicle tailpipes, the Chinese government has continually updated vehicle emission standards, as well as electrification and fuel efficiency standards. From 2000 to 2023, China successively introduced six national motor vehicle pollutant emission standards [6]. These standards aim to tighten emission restrictions and reduce individual vehicle emissions. Simultaneously, the Chinese government is actively promoting the development of eco-friendly cities, creating more compact urban communities conducive to nonmotorized transportation, and focusing on the rational distribution of urban neighborhoods [7]. These efforts aim to further reduce urban transportation GHG emissions in the future. In recent years, there has been a gradual increase in research focusing on strategies to reduce urban GHG emissions through the encouragement of green travel modes and the optimization of urban land use structures. Several researches have shown that continuous improvement of urban land use structures and the built environment (BE) in the process of urbanization plays a crucial role in reducing urban GHG emissions [8,9,10]. The BE becomes an important factor influencing the reduction of GHG emissions from travel by urban residents.

However, when estimating GHG emissions, some studies have not made detailed calculations based on the National Motor Vehicle Pollutant Emission Standards and available vehicle trajectory data and have provided simple estimates using Euclidean distances multiplied by emission factors [7]. Moreover, many studies have investigated the nonlinear influence of different BE elements on GHG emissions, but the spatiotemporal heterogeneity of the nonlinear influence on GHG emissions requires further exploration. The aim of this study is to provide more accurate calculations of taxi GHG emissions and to offer new research perspectives on the spatial heterogeneity of the nonlinear relationship between BE and GHG emissions. Our research fills the aforementioned research gaps by addressing two key questions: (1) Which BE elements are most crucial for reducing the GHG emissions in different urban functional areas? and (2) Is there heterogeneity in the nonlinear influence of the BE on the GHG emissions between weekdays and weekends and across different urban areas? This study distinguishes between weekdays and weekends to explore the spatiotemporal heterogeneity of the nonlinear impacts of BE variables on the GHG emissions of taxis in various urban areas.

Using the China National Motor Vehicle Pollutant Emission Standard V and taxi trajectory dataset, the GHG emissions of urban taxis are accurately estimated in this study. Furthermore, based on a clustering method of traffic analysis zones (TAZs) with similar spatiotemporal emission characteristics and with reference to the study of Ding et al. [11], who distinguished between weekdays and weekends, the impact factors of GHG emissions are spatiotemporally explored. Figure 1 illustrates the framework of this study. Our study reveals the spatiotemporal differences in GHG emissions from urban taxis and the important impact factors of the BE that contribute to these differences; specifically, an in-depth study of the nonlinear associations and spatiotemporal heterogeneity between important impact factors and GHG emissions from taxis is performed. Our research findings are significant for planners, as it is essential for them to identify which BE variables have the greatest impact on GHG emissions from urban transport in each functional area of the city as they establish sustainable urban plans. Additionally, optimizing the urban BE and adjusting land use to formulate taxi GHG emission-reduction policies can reduce not only the GHG emissions generated by taxis but also, to a certain extent, can change the travel patterns and population mobility patterns of urban transportation and thus reduce the CO₂ emissions of urban residents’ trips [12].

The remainder of the paper is organized as follows: In the literature review, research on GHG emissions in the transportation sector is discussed. Section 3 describes the extraction of the research area, data preprocessing, and the selection of variables for analysis. Section 4 provides a detailed explanation of the research methodology. The results of the GHG emission measurements and the spatiotemporal heterogeneity assessment of nonlinear impact factors are presented in Section 5. Finally, we present conclusions and recommendations in Section 6.

2. Literature Review

Taxis are primarily categorized into two main types within urban intelligent transportation systems: online-hailing taxis and cruising taxis. In urban transport systems, taxis are the vehicles with the longest operating period, the largest operating range, and the most miles travelled [13]. Taxi trajectory data provide comprehensive and timely feedback on the operational characteristics and traffic status of urban transportation, comprising massive spatiotemporal information on urban human activity and mobility [14]. Taxi drivers typically possess extensive driving experience and strive for prompt order completion during their service, emphasizing the “fast-as-possible” approach. Consequently, the spatiotemporal characteristics of taxi operations can offer a relatively true reflection of urban traffic conditions and residents’ travel patterns [15]. Furthermore, taxi drivers draw upon their subjective experiences during passenger searches, rendering taxi trajectories susceptible to factors such as time constraints, congestion, road network density, and other limitations of the BE. GHG emissions attributed to taxis are generally divided into direct emissions and indirect emissions. Indirect emissions result from the increased traffic congestion caused by taxis during their operation and idling periods [16]. In this study, our primary focus is on the direct emissions produced during taxi operations. Direct GHG emissions (such as NO_x) produced by taxis account for approximately 10% of GHGs from urban transportation [17]. Wang et al. [18] argue taxis are the least emission-efficient mode of transportation. In terms of the CO₂ emissions per passenger-kilometer travelled within the residential areas of Beijing, the CO₂ produced by taxis are five times as much as those produced by bus travel. Therefore, accurately estimating GHG emissions from urban taxis and promoting emission reductions are vital tasks for improving the overall GHG emission efficiency and advancing the development of green cities.

Research on reducing GHG emissions in the transport sector, aiming to achieve carbon neutrality, has been thoroughly explored by numerous scholars in recent years. These researches primarily focus on measuring GHG emissions over the whole life cycle of road construction [19,20] and estimating GHG emissions from various modes of travel in urban transport [21,22]. Taxi GPS data possess characteristics such as wide coverage, comprehensive information, and precise spatial positioning, providing a robust foundation for research. Existing studies that utilized taxi trajectories for GHG emissions research primarily focused on the estimation and spatiotemporal analysis of urban transportation CO₂ emissions [23,24,25], measurement of emission factors for cars [26,27], investigation of CO₂ reduction through the ridesharing behaviors of ride-hailing passengers [28], and investigation of the elements that effect carbon emissions from urban residents’ travel [29]. These researches have effectively harnessed the benefits of taxi trajectory data to study micro-urban/road perspective reduction strategies and the factors influencing GHG emissions in urban transportation; however, these studies did not consider elements of the BE and nonlinear influences comprehensively.

The built environment (BE) encompasses urban land, transportation infrastructure, public service facilities, and other human-made spaces intertwined with people’s daily lives, work, and leisure [30]. After extensive research by numerous scholars since the 1990s, the BE has been conclusively conceptualized as a “5Ds” framework on the basis of density, diversity, design, distance to transit, and destination accessibility [31]. Urban planning typically precedes transportation planning, and the BE can profoundly impact transportation planning and residents’ travel behaviors. Therefore, the effects of the BE and land use on urban transport development have emerged as popular research topics. Exploring the patterns of taxi GHG emissions and their associations with the BE is crucial for planners to determine the trade-off between the urban function building and urban traffic emissions in the process of urbanization. Chen et al. [32] employed a semiparametric GWPR to explore the linear influence of the BE on urban taxi passenger volume at a 0.5 km² grid scale. The findings indicated that the influence of the BE on taxi order volume exhibits spatial heterogeneity in different urban areas, implying that the BE in distinct regions could have either positive or negative effects on taxi passenger volume. Based on data from two categories of passenger vehicles, Tan et al. [33] calculated CO₂ emissions from urban passenger transport and found spatial differences in these emissions. The main BE factors contributing to this spatial distribution difference are mixed land-use layouts and destination accessibility. Kwak et al. [34] concluded that the BE has a major impact on urban transport GHG emissions. Ashik et al. [35] investigated the relationship between the BE and commuting CO₂ emissions using the structural equation model (SEM) with vehicle ownership as an intermediate factor. In terms of the association between the BE and urban transportation GHG emissions, many studies have shown that the BE affects GHG emissions by influencing motor vehicle ownership and altering the travel patterns of residents [10]. Therefore, the studies mentioned above encompass the mechanism of the influence of the BE on GHG emissions from urban transportation. In the literature, aggregate, disaggregate, and multiscale studies have investigated the relationship between the BE and GHG emissions. Aggregate studies have assessed GHG emissions and the BE at an aggregate scale, such as a metropolis, district, or block, to investigate their associations [36]. Disaggregate studies have employed individual- or household-scale resident travel surveys and questionnaire data [37]. Furthermore, Yang [38] combined the neighborhood scale, subdistrict scale, and 1 km² grid scale to comprehensively investigate the influence of the BE on commuting-based GHG emissions at multiscale. Although these studies encompassed various perspectives and scales, they did not explore the spatiotemporal heterogeneity of the relationship between the BE and GHG emissions in different urban areas.

In regard to the study approach for assessing the influence of the BE on transportation GHG emissions, studies have employed both linear and nonlinear models to analyze influencing factors. Linear models employ linear regression analysis to determine the magnitude of impact of various factors by examining the corresponding variable coefficients. These models use spatial econometric models and their enhanced models. Nonlinear models mainly use ensemble learning algorithms including GBDT and XGBoost to analyze nonlinear associations and thresholds between BE elements and urban transportation. Both linear and nonlinear models have respective advantages and disadvantages. Linear models can be used to identify spatial heterogeneity in impacts, and the coefficients of independent variables often display statistical significance, with interpretable p values. However, linear models overlook the nonlinear relationships between dependent and independent variables. The advantage of nonlinear models is that they can be used to more comprehensively explore the associations between dependent and independent variables. Ding et al. [39] argued that ignoring the effect of nonlinear associations may lead to erroneous conclusions and fail to uncover the impacts of some key features on predictor variables. Therefore, we use a nonlinear model to investigate the associations between the BE and the GHG emissions to comprehensively explore the influence mechanism of the BE on the GHG emissions from taxis.

From the above literature review, it is evident that current research on the association between the BE and GHG emissions from urban transportation lacks analyses based on nonlinear models, which are needed to consider the spatiotemporal heterogeneity of the nonlinear influence of the BE on taxi GHG emissions at the urban level on weekdays and weekends. Many studies have focused primarily on various impact factors at the urban spatial scale, failing to consider the differences in GHG emissions at different spatiotemporal scales at the TAZ level and to investigate the underlying reasons for these differences. Therefore, we find that there are certain inadequacies concerning knowledge of the nonlinear impact mechanisms of GHG emissions from urban functions areas at various temporal and spatial scales. Due to the complex characteristics of GHG emissions from taxis, which involve various travel modes, such as commuting, leisure, and social activities, studying taxi GHG emissions provides a basis for comprehensively assessing urban GHG emission characteristics and the corresponding impact factors.

3. Data and Variables

3.1. Data Processing

3.1.1. Research Area Definition and Traffic Analysis Zone (TAZ) Subdivision

The dataset in this study is the taxi operation GPS trajectory data in Xi’an, China. Six districts of the main city of Xi’an are chosen as the research area. The longitude range of this area is 108.787–109.274°, and the latitude range is 34.164–34.458°; this area covers approximately 826 km².

For simplicity, the research area is divided into TAZs based on the theory of TAZ subdivision. Taxi trajectories are absent in smaller-scale urban functional areas such as parks, school zones, and neighborhoods. Therefore, this study posits that when investigating GHG emissions from taxis, the subdivision scale of TAZs should be moderately expanded beyond that of conventional TAZs. This expansion should account for factors such as urban road networks, infrastructure layouts, and the distribution of the BE, integrating the principles of traditional TAZ division theory. Figure 2 displays the research area and the ultimate segmentation of the 268 TAZ subdivisions, and

Table 1. Statistical description of the TAZ areas.

Statistical Indicators	Mean	Min	25% Quantile	Median	75% Quantile	Max
TAZs area /km2	4.58	0.25	2.67	3.83	5.22	27.86

presents the statistical indicators of the area of each TAZ.

3.1.2. Data Preprocessing Process

The taxi trajectory dataset employed in this research originates from the Xi’an Taxi Management Bureau. Our research employs data from 2019 to estimate GHG emissions to minimize the effects of COVID-19 waves during the post-2020 lockdown. GPS trajectory data for taxis operating in Xi’an, China, were collected between 28 February 2019 and 6 March 2019. The structure and field descriptions of the dataset are given in

Table 2. The data formats.

Field	Format	Field Description
ID	A75589 or UD33606	Vehicle registration number
IN_DATE	2019/3/2 00:32:03	Data entry time
GPS_TIME	2019/3/2 00:00:30	Trajectory point recording time
LNG	108.888170	Longitude
LAT	34.268991	Latitude
HEIGHT	436	Altitude/m
SPEED	50	Speed/km·h−1
EFF	0 or 1	0 record invalid; 1 record valid
CAR_STATE	4, 5, 7, 9	4 Empty vehicle status; 5 Passenger loaded status; 7 Engine off; 9 Abnormal

.

Based on data from the Shaanxi Road Transportation Service Development Center (https://jtyst.shaanxi.gov.cn/dlys, accessed on 11 July 2024) from 2019, the number of registered taxis in Xi’an is 16,526, of which the number of electric vehicles is 7406. Since electric vehicles do not produce GHG emissions during operation, we exclude electric vehicles (vehicles containing string “D” in the second position of the ID, and ID length is seven characters) from the study. To ensure the integrity and accuracy of this study, it is crucial to preprocess the data collected by vehicle-mounted devices. This is necessary due to the likelihood of encountering outliers, redundant information, and noise in the dataset. The preprocessing steps, as depicted in Figure 3, are designed to remove these irregularities and maintain the quality of the data.

After Figure 3 processing, the number of fuel-consuming taxi vehicles with reliable data on a daily basis during the study period is 6808, 6807, 6780, 6783, 6878, 6670, and 6552, respectively.

3.2. Variable Selection

We aggregate the calculated GHG emissions data using Arcgis 10.8 and GeoPandas, following the method outlined in Section 4.1.2, into the predivided TAZs. Considering the different areas of TAZs, we use the mean value of GHG emissions from taxis on weekdays and weekends per km² within each TAZ as the dependent variable. The specific calculation steps and results are shown in Section 4.1.2 and Section 5.1.

We categorize and integrate the “5Ds” elements of the BE as independent variables. Existing research on the delineation boundaries of “5Ds” elements in BE analysis shows ambiguity and uncertainty [40]. Some dimensions overlap; for instance, Sun and Quan [41] consider bus stop density as one of the indicators measuring destination accessibility, whereas Weng et al. [42] use both bus stop density and metro station density as elements of distance to transit. Therefore, in the selection of the variables “distance to transit” and “destination accessibility”, we apply the approach of Ji et al. [43] and select “metro/bus stop density” as an independent variable. Additionally, considering the presence of taxi stops within the research area, the final variable is determined as “metro/bus/taxi stops density”. The variables we select encompass a range of factors within the 268 TAZs in the research area, such as the density of various types of functional places, the road system, and the density of public transportation infrastructure sites. In previous research, the influence of road network density on CO₂ emissions from urban transportation has shown both positive and negative associations. Wang et al. [44] argued that a higher road network density improves accessibility within a region, which decreases CO₂ emissions in the region, while Shu and Lam [45] argued that a high road network density is more likely to cause vehicle congestion and thus enhance CO₂ emissions in the region. Our research investigates and validates these phenomena from the perspective of the TAZ scale. The road network dataset employed for the study are from OpenStreetMap and are in the WGS-84 coordinate system. The data for independent variables are sourced from the Gaode Map Open Platform, which consists of 203,277 records of point of interest (POI) data. Andrade et al. [46] indicated that using POI data to characterize a city’s functional zones and land use is a reasonable approach. In comparison to land use maps, POI data can provide a more comprehensive representation of various types of land use for different functional areas.

The formula for Road network density in TAZ numbered i is:

$${\rho }_i={\displaystyle \raisebox{1ex}{$L_i$}\!\left/ \!\raisebox{-1ex}{$S_i$}\right.}$$

(1)

where ${\rho }_i$ is the density of the road network of the TAZs numbered i, $L_i$ is the length of all roads in the TAZs numbered i, and $S_i$ is the area of the TAZs numbered i.

Mixed degree of urban function is the entropy of all types of POIs; its computational formula as follows:

$$H_i=-{\sum }_j^k{p}_{\mathit{ij}}\mathit{log}{p}_{\mathit{ij}}$$

(2)

where $H_i$ is mixed degree of urban function of TAZs numbered I, $H_i\ge 0$; $p_{\mathit{ij}}$ represents the proportion of category j POIs in TAZs i relative to the total number of POIs, 0 ≤ $p_{\mathit{ij}}$ ≤ 1. A larger $H_i$ indicates a greater abundance of functional types within TAZs and a smaller difference in the number of each functional type.

The raw data obtained in Arcgis10.8 integrate in TAZs to calculate the values of the independent variables.

Table 3. Statistical analyses for BE variables.

Category	Variable	Description	Mean	Min	Max	Std.dev
Design	Road network density	——	9.489	0	33.331	4.960
	Intersection density	Intersections count per square kilometer within TAZ	15.090	0	80.391	11.826
Diversity	Mixed degree of urban function	——	1.722	0	2.259	0.438
Density	Restaurant density	Restaurants count per square kilometer within TAZ	68.0128	0	498.241	100.569
	Scenic area density	Scenic areas, park, temples, etc. count per square kilometer within TAZ	1.272	0	9.234	4.511
	Communal facility density	Communal facilities count per square kilometer within TAZ	2.737	0	28.059	4.025
	Shopping use density	Malls, supermarkets count per square kilometer within TAZ	110.822	0	816.667	155.529
	Education facility density	Schools, universities, and other tutoring centers count per square kilometer within TAZ	13.989	0	229.735	26.820
	Finance use density	Insurance company, bank, etc. count per square kilometer within TAZ	5.797	0	89.024	10.314
	Sport use density	Stadium, sports field count per square kilometer within TAZ	5.983	0	57.339	9.276
	Service of life density	Moving companies, courier collection points, telecom business halls, etc. count per square kilometer within TAZ	56.331	0	355.759	79.359
	Residential density	Residential buildings count per square kilometer within TAZ	9.600	0	58.283	13.530
	Accommodation service density	Hotels count per square kilometer within TAZ	11.317	0	201.889	21.983
	Government use density	Government agencies count per square kilometer within TAZ	6.798	0	90.353	11.731
	Medical Institutions density	Hospitals, pharmacies count per square kilometer within TAZ	9.065	0	70.748	13.633
Distance to transit/ Destination accessibility	Metro/bus/taxi stops density	Bus, metro, and taxi stations count per square kilometer within TAZ	17.517	0	108.807	22.976

displays the detailed descriptions and statistical analyses of each variable.

4. Methodology

4.1. Vehicle GHG Emissions Estimation Model

The GHGs emitted from the transportation sector, which are measured internationally, include not only CO₂ but also CO, CH₄, and NO_x [47]. The “bottom-up” methods provided by the IPCC start from the micro perspective, taking various types of transport as the main body for energy-consumption statistics, i.e., employed the mileage of various transport modes and the energy consumption per kilometer measure the total energy consumption [48].

According to the test from bench and vehicle emissions, the results of the COPERT model based on the “bottom-up” method are closer to the light vehicle emission situation in the urban area of China [49]. The COPERT model was developed by the European Environment Agency and is applicable to vehicles meeting European emission standards. Historically, China’s motor vehicle pollutant emission standards have generally followed European emission standards. The emission standards I–V correspond to the European motor vehicle emission standards I–V [50]. Therefore, the COPERT model has been used by many scholars for the measurement of transportation-based GHG emissions in China [51]. The COPERT model is chosen to estimate the GHG emissions of operating taxis in the research area.

4.1.1. COPERT Model Assumptions

The COPERT model is calibrated with extensive experimental dataset for the parameters related to various types of pollutant emissions by combining numerous indicators, such as motor vehicle GHG emission standards, vehicle category, and fuel category. The parameters calibrated by the COPERT model are selected to measure the GHG emissions from urban taxi operation in this study at the microscale, and the basic assumptions in the measurement process are as follows: (1) Environmental climate change and weather changes are assumed to have negligible effects on the calibration of the parameter. (2) Prolonged vehicle use and mechanical damage may affect the calibration of parameters, so the duration of use and damage to the vehicle are not considered in this study. (3) The impacts of individual driving behaviors and driving habits on the model parameters are not considered. (4) The COPERT model classifies the emission pattern of engines into cold-start emissions or hot emissions. Because records with an invalid trajectory state at start-up are removed during the data analysis process in this study, only hot emissions from taxis are considered. (5) Vehicle emission standards are environmental regulations aimed at reducing tailpipe emissions from motor vehicles. They entail stricter emission regulation through measures such as regulating pollution control devices, enhancing emission-reduction technologies in vehicles, and setting limits on pollutant emissions. These standards play a crucial role in reducing vehicular emissions, addressing air pollution, and combating haze. The higher the emission standard level is, the less GHGs are emitted from the vehicle model based on that emission standard compared to the GHGs emitted if the previous standard is used.

The research area of Xi’an, China, has been implementing China’s motor vehicle pollutant emission standard V since 1 January 2017, and the implementation of China’s motor vehicle emission standard VI began on 1 July 2019. The study period of this paper is from 28 February to 6 March 2019, so we assumed that standard V is followed for estimating taxi energy emissions during the study period.

4.1.2. Model Estimation Process

In this paper, the GHG hot emissions from taxi engines are estimated as follows:

$$E_{q,r}=F_{q,r}{d}_r$$

(3)

where E_q,r is the emission of GHG q in taxi trajectory segment r, d_r is the mileage travelled by the taxi in segment r, and F_q,r is the vehicle emission factor for GHG q in taxi trajectory segment r. In this model, the vehicle GHG emission factors are related to vehicle category, emission standard, fuel category, engine technology standard, type of roadway on which the vehicle travels, and vehicle operating speed, of which the first five categories can be determined by the calibration parameters of the model, and the vehicle operating speed has to be obtained through the calculation of trajectory data. We use the speed of the segment to calculate the emission factor, given by:

$$F_{q,r}={\displaystyle \frac{{\alpha }_q{v_r}^2+{\beta }_q{v}_r+{\gamma }_q+\frac{{\delta }_q}{v_r}}{{\epsilon }_q{v_r}^2+{\zeta }_q{v}_r+{\eta }_q}}$$

(4)

where v_r (km/h) is the mean speed of the taxi on segment r, α_q, β_q, γ_q, δ_q, ε_q, ζ_q, η_q are the parameters of the model, and q is the four types of GHGs. We select the GHG emission calculation parameters for the Euro 5 motor vehicle emission standards that correspond to the taxi types in Xi’an, as outlined in

Table 4. COPERT model parameters.

	αq	βq	γq	δq	εq	ζq	ηq
CO	0.00045	−0.10208	6.87693	10.38385	0.00162	−0.43756	30.33733
NOx	−0.00031	0.10306	0.23906	−0.33928	0.03454	1.98601	1.26376
CH4	0.00000	0.00000	2.87000	0.00000	0.00000	0.00000	1000.00000
CO2	0.34656	−18.30053	1513.54426	0.00000	0.00080	0.09133	3.51264

The parameters of the CO₂ are derived by converting the national energy consumption factor (MJ/km) of the People’s Republic of China from the Emissions Report Overview 1.4 on the COPERT website, with a conversion factor of 72.24 g/MJ.

.

In the GHG emission calculations for TAZs, continuous trajectories of taxis are assigned to all TAZ interiors using the spatial intersection approach in Arcgis 10.8. The A-B trajectories and C-F trajectories (shown in Figure 4) in three of the TAZs numbered 8, 222, and 228 are used as examples for illustration. Trajectory segment AB is located entirely within TAZ number 222, so emissions from trajectory segment AB are assigned to TAZ 222. The C-F trajectory segment spans three TAZs, namely, 222, 8, and 228, so the trajectory segment CF is split into three trajectory segments, namely, CD, DE, and EF, which are assigned to TAZs 222, 8, and 228, respectively. The real travelled mileage in the assigned trajectory segments and the operating time are calculated from the dataset, and running speeds are computed for CD, DE, and EF. With CO₂ emissions as an example, the CO₂ emissions for the AB trajectory are calculated as follows:

$$E_{C{O}_2,\mathit{AB}}=F_{C{O}_2,\mathit{AB}}{d}_{\mathit{AB}}$$

(5)

The CO₂ emissions for the CF trajectory are as follows:

$$E_{C{O}_2,\mathit{CF}}=E_{C{O}_2,\mathit{CD}}+E_{C{O}_2,\mathit{DE}}+E_{C{O}_2,\mathit{EF}}=F_{C{O}_2,\mathit{CD}}{d}_{\mathit{CD}}+F_{C{O}_2,\mathit{DE}}{d}_{\mathit{DE}}+F_{C{O}_2,\mathit{EF}}{d}_{\mathit{EF}}$$

(6)

The CO₂ emissions in TAZs numbered 8, 222, 228, respectively:

$${\mathit{TAZ}}_{{\mathit{CO}}_2,8}=E_{C{O}_2,\mathit{DE}}$$

(7)

$${\mathit{TAZ}}_{{\mathit{CO}}_2,222}=E_{C{O}_2,\mathit{AB}}+E_{C{O}_2,\mathit{CD}}$$

(8)

$${\mathit{TAZ}}_{{\mathit{CO}}_2,228}=E_{C{O}_2,\mathit{EF}}$$

(9)

4.2. GBDT Model

An ensemble GBDT model is employed to identify the important impact factors in this research, including nonlinear factors, that influence GHG emissions from taxis. An ensemble algorithm merges multiple machine learning models to build a more powerful model, and it has been shown that ensemble models are effective for processing large amounts of regression data. The GBDT is an ensemble model based on decision trees that has been increasingly applied to research related to urban transportation [52]. Wu et al. [53] encourage planning scholars to apply GBDT to identify the nonlinear relationship and threshold between BE and traffic CO₂ emissions. In contrast to other machine learning algorithms, GBDT can identify and rank the influences of features on predictor variables. GBDT merge many weak learner trees with shallow depths and improve modeling performance by increasing the number of trees in a continuous iteration. This approach is suitable for datasets in which binary features coexist with continuous features [54].

GBDT can analyze the nonlinear association between a single feature and the predictor variable and overcome the problem of multicollinearity in traditional linear models [55]. Furthermore, modern deep learning algorithms are similar to “black boxes”, and many cannot be used to assess the importance and trend of each input feature. Compared to deep learning algorithms, the advantage of a GBDT is that it can deal with high-dimensional multi-feature data and does not have to perform feature selection; after model training is completed, the level of feature importance is determined. In this study, because we first conduct spatiotemporal clustering in 268 traffic zones, the use of the GBDT model allows us to avoid standardizing the data for each cluster individually.

GBDT uses the negative gradient of the loss function as an approximation of the residuals from the previous round of learners, constructs a new decision tree in the direction of the gradient where the residuals decrease, and optimizes the general loss function by approximating the residuals through the negative gradient. Mathematically, let X = {x₁,x₂,…,x_n} represent the independent variables in this study, primarily consisting of 16 indicators related to road systems and land use, while y serves as the dependent variable, i.e., taxi GHG emissions per unit km² within each TAZ. $\widehat{y}\left(x\right)$ is an approximate function of y expressed in terms of X. The training set is set to D = {(x₁,y₁), (x₂,y₂), …, (x_N,y_N)}. The aim of GBDT is to set a function $\widehat{y}\left(x\right)$ with M decision trees such that the loss function $L\left[y,\widehat{y}\left(x\right)\right]={[y-\widehat{y}(x\left)\right]}^2$ is minimized. $\widehat{y}\left(x\right)$ can be formulate as follows:

$$\widehat{y}(x)={\sum }_{m=1}^M{\widehat{y}}_m(X)={\sum }_{m=1}^M{\theta }_{m}T(X,A_m)$$

(10)

where A_m is the mean of split positions and terminal nodes in an individual tree $T\left(X,A_m\right)$, ${\theta }_m$ is the optimal solution parameter that minimizes the loss function locally or globally by model convergence, and the estimation of the parameter is carried out by gradient boosting, as shown in the following procedure.

First, initialize ${\widehat{y}}_0\left(x\right)$ as follows:

$${\widehat{y}}_0\left(X\right)={\mathit{arg}\mathit{min}}_{\theta }{\sum }_{i=1}^{N}L\left(y_i,\theta \right)$$

(11)

Second, iterate over each independent variable in the dataset and compute its negative gradient u_im as the residual.

$$u_{\mathit{im}}=-{\left[{\displaystyle \frac{\partial L(y_i,\widehat{y}(x_i))}{\partial \widehat{y}\left(x_i\right)}}\right]}_{\widehat{y}\left(x_i\right)={\widehat{y}}_{m-1}\left(x_i\right)}$$

(12)

Based on u_im the optimal $\theta$ can be obtained by fitting a decision tree $T\left(X,A_m\right)$:

$${\theta }_m={\mathit{arg}\mathit{min}}_{\theta }{\sum }_{i=1}^{N}L(y_i,{\widehat{y}}_{m-1}(x_i\left)+\theta T\right(X,A_m\left)\right)$$

(13)

By updating the equations with the above Equations (12) and (13) and adding the learning rate a (0 < a < 1) to prevent overfitting in the approximation function to prevent overfitting, the final model is obtained as:

$$\widehat{y}\left(X\right)={\widehat{y}}_{m-1}\left(X\right)+a{\theta }_{m}T\left(X,A_m\right)$$

(14)

GBDT is sensitive to the setting of parameters, so it is necessary to carry out appropriate hyperparameter tuning operations during use. The main parameters for hyperparameter tuning in this study are the tree numbers (n), tree depth (max_depth), and learning rate. The weekday model and the weekend model are tuned separately to select the appropriate GBDT parameters. In the hyperparameter tuning step, the appropriate range of n is first selected and then traversed for different learning rates. Then, max_depth is determined to reduce the complexity of each tree. Finally, the parameter values that yield the smallest output of the loss function L are selected. The final hyperparameter tuning results are optimal for the weekday dataset with the following parameters: learning rate = 0.002, max_depth = 2, and n = 5000. For the weekend dataset, the optimal hyperparameters are as follows: learning rate = 0.001, max_depth = 3, and n = 20,000. Therefore, in our weekday model, we set a maximum of 5000 trees with a depth of 2 and a learning rate of 0.002; on weekend model, we set a maximum of 20,000 trees with a depth of 5 and a learning rate of 0.001. The test-R² values of the weekday model and weekend model are approximately 0.82 and 0.79, respectively.

After training the GBDT model and performing hyperparameter tuning, we estimate the feature importance of each feature (i.e., an independent variable) in relation to the response variable (i.e., the dependent variable) in the prediction and apply the feature importance score (FIS) to rank the importance of the factors influencing the GHG emissions of taxis, reflecting the nonlinear effect of each type of impact factor on the GHG emissions of taxis. The squared importance of each feature $I_{x_i}^2$ can be calculated by averaging the squared importance of all additive trees as follows:

$$I_{x_i}^2={\displaystyle \frac{1}{M}}{\sum }_{m=1}^M{I}_{x_i}^2\left(T_m\right)$$

(15)

$I_{x_i}^2$ i.e., FIS. After controlling for other features variables, the GBDT can plot a partial dependence plot (PDP) to visualize the nonlinear association between individual features variables and the GHG emissions of taxis [56].

5. Results and Discussion

5.1. Estimation of GHG Emissions

Using the model, methodology, and data mentioned above, the emissions of CO, NO_x, CH₄, and CO₂ from taxis in the research area from 28 February 2019 to 6 March 2019, were estimated, as shown in Figure 5.

The day with the highest emissions was Friday, 1 March 2019; the total amount of GHG emissions from the four types of GHGs during taxi operation was 418,325.40 kg, and the emissions of CO, NO_x, CH₄, and CO₂ were 561.13 kg, 69.14 kg, 6.09 kg, and 415,380.01 kg, respectively. The day with the lowest emissions was Sunday, 3 March 2019; the total amount of GHG emissions during taxi operations was 395,596.56 kg, and the CO, NO_x, CH₄, and CO₂ emissions were 531.94 kg, 65.69 kg, 5.76 kg, and 394,993.17 kg, respectively. The main reason for the variation in GHG emissions throughout the week is the differing total travel kilometers by taxis. For instance, on the day with the highest total GHG emissions, Friday, March 1, taxis in the research area covered a total distance of 2,128,714 km. In contrast, on the day with the lowest total emissions, Sunday, 3 March, taxis travelled 2,006,552 km, significantly less than on Friday. Furthermore, we illustrated the spatial distribution differences by plotting heat maps for the highest and lowest emission days using CO₂ as an example in Figure 6. It can be observed that there are significant differences in taxi CO₂ emissions in the northern and southwestern parts of the research area. To study the temporal distribution patterns of GHG emissions from taxis, we divided each day of the study period into 24 periods on an hourly basis, and the emission patterns of the four types of GHGs in each period are shown in Figure 7.

The hourly emission patterns of the GHGs display similar trends throughout the duration of the study in Figure 7. Overall GHG emissions are divided into three phases: midnight emission peak (23:00–0:00), A.M. emission peak (7:00–9:00), and midday emission peak (12:00–14:00). The midday emission peak occurs primarily on weekdays due to the behavior of short commute trips that may occur during weekday lunch breaks. There is a midnight peak in taxi emissions in Xi’an, attributed to the significant nighttime recreational activities. The reason for the low taxi GHG emissions in Xi’an during the P.M. peak period is that from 16:00 to 18:00, when taxi drivers change shifts, the number of vehicles operating decreases. The midnight emissions on weekdays are less than those on weekends, and the A.M. and midday peak emissions on weekdays are significantly greater than those on weekends.

In analyses of the total amount of GHG emitted from vehicles, the unit emission factor of each GHG in a trajectory segment, i.e., the amount of GHG emitted per kilometer travelled by a vehicle, is often considered. Since the emission trends of the four types of GHGs are generally consistent and the international focus on GHG emissions in the transportation sector is primarily related to CO₂, researchers both domestically and abroad have conducted extensive studies on the unit emission factor of CO₂. Additionally, as demonstrated in the previous calculations, CO₂ has the highest emission volume among the four GHGs. Therefore, in the subsequent sections, we focus on CO₂ as an illustrative example. We use CO₂ unit emission factors as an example of unit emission factors for 24 periods in the research area. The CO₂ unit emission factor correlates with the vehicle travel speed, the lower or higher the speed is, the greater the emission factor. When the taxi operating speed is low in areas of urban traffic congestion, frequent vehicle acceleration, deceleration, idling, and mode switching cause a sharp increase in the unit emission factor of vehicles.

The CO₂ unit emission factor of taxis in the research area was calculated as 0.1964 kg/km, among which the emission factor for weekdays is 0.1957 kg/km and the emission factor for weekends is 0.1979 kg/km, suggesting that the unit emission factor for weekends is greater than that for weekdays. Our results are slightly lower than the average value of 0.215 kg/km for gasoline-powered passenger cars in Shanghai, which aligns with China’s motor vehicle pollutant emission standard IV.

$$f_{{\mathit{CO}}_2,t}={\displaystyle \raisebox{1ex}{${\sum }_{r=1}^n{e}_{{\mathit{CO}}_2,r,t}$}\!\left/ \!\raisebox{-1ex}{${\sum }_{r=1}^n{d}_{r,t}$}\right.}(t=0,1,\dots ,23)$$

(16)

Equation (16) is used to calculate the hourly taxi CO₂ unit emission factor during the study period 28 February 2019 to 6 March 2019, where $f_{{\mathit{CO}}_2,t}$ is the CO₂ emission factor at hour t, $e_{{\mathit{CO}}_2,r,t}$ is the CO₂ emission of trajectory segment r at hour t, and $d_{r,t}$ is the distance travelled by the taxi in trajectory segment r at hour t. The variation of the calculated CO₂ emission factor per hour for taxis is shown in Figure 8.

Notably, the peak CO₂ unit emission factor for each day occurs in the P.M. peak period, and the lowest value occurs at 00:00 midnight, but due to the high demand and trajectories of taxi travel in the research area at night, the CO₂ emissions from taxis at midnight are characterized by a small emission factor and high emission levels. On weekdays, the CO₂ unit emission factor shows clear peaks from 8:00–9:00 A.M. and 18:00–19:00 P.M.; on weekends, the CO₂ unit emission factor shows a rising trend before the P.M. peak period, with a clear peak at midday, and the emission factor on Saturdays is significantly greater than that on Sundays. Among the P.M. peak periods, the largest CO₂ unit emission factors are found during the Friday P.M. peak period and Saturday P.M. peak period, and the CO₂ unit emission factor is significantly greater on Friday evening from 20:00–23:00 than at other times. The maximum CO₂ unit emission factor in the A.M. peak period occurs on Monday. Figure 8 shows obvious heterogeneity in the temporal characteristics of CO₂ unit emission factors for taxis in the research area on weekdays and weekends.

5.2. Spatiotemporal Characteristics of GHG Emissions in TAZs

Based on the TAZs established in Section 3.1.1, the CO₂ emission data calculated for all taxi trajectories during the 7 days of the study period are aggregated. Figure 9 shows a characteristic heat map of the distribution of CO₂ emissions from taxis for each TAZ for A.M. peak, midday peak, P.M. peak, and midnight peak periods on weekdays and weekends. The notable phenomenon is that the white TAZ numbered 11 is surrounded by orange or red TAZs. This TAZ refers to the Daming Palace National Heritage Park located in the central urban area, where taxis are prohibited from entering at all times; they can only operate on the roads surrounding the park. On both weekdays and weekends, TAZs in the northern part of the research area, as well as TAZs in the southern region, exhibit high levels of taxi CO₂ emissions during the A.M. peak, midday peak, and P.M. peak periods.

During the A.M. peak period, the areas where the difference in CO₂ emissions from taxis between weekends and weekdays is large are located in the TAZs near South Second Ring Road in the city center, where large commercial centers and tourist attractions are located; therefore, the weekend CO₂ emissions of taxis are significantly greater than those on weekdays in these TAZs. The distributions of CO₂ emissions from taxis on weekdays and weekends during the midday peak period are basically the same, with CO₂ emissions from taxis on weekdays being slightly greater than those on weekends; the corresponding high-emission areas are mainly distributed along Metro Line 2. During the P.M. peak period, CO₂ emissions from taxis are slightly lower on weekdays than on weekends. Due to the many types of travel options available during the P.M. peak off-work commuting travel periods and the high level of traffic congestion, people in areas with high population densities and better public transport infrastructure try to avoid using taxis to travel, and the majority of the residents choose to travel by public transport, resulting in lower emissions from taxis. The trips during the P.M. peak period on weekends are mostly for family leisure and recreation, with more residents choosing taxi travel at this time than on weekdays. During the midnight peak period, the high CO₂ emission TAZs are located in the south-central part of the research area, and the distribution characteristics on weekdays and weekends are relatively similar. These areas include scenic areas with famous night-time tourist attractions, such as Yong Ning Gate, Tang Dynasty Never Nights City, and Big Wild Goose Pagoda, as well as night-time recreational facilities, such as bars. When public transport ceases operations at midnight, taxi travel becomes the first choice for tourists travelling in these areas.

From the above analyses, we can clearly see significant spatiotemporal distinctions in GHG emissions within each TAZ during the peak hours. In the next part of this study, the reasons for such differences are further analyzed. Since there are differences in land use, BE, and infrastructure in the TAZs, all these factors may have an impact on the spatiotemporal distribution of GHG emissions from taxis. Next, we use a clustering algorithm to cluster TAZs with relatively similar spatiotemporal GHG emission characteristics and investigate the factors that cause spatiotemporal heterogeneity in GHG emissions from taxis in TAZs based on the clustering results.

5.3. Spatiotemporal Clustering Analysis of TAZs Based on GHG Emissions

In this study, an unsupervised machine learning algorithm, the k-means algorithm, is used for cluster analysis to classify TAZs into different clusters based on their emission patterns. The k-means algorithm is widely used in clustering studies related to TAZs [57,58]. We divided the study period into daily, 24-h periods and continued our investigation using CO₂ emissions as an example. Due to variations in the area of each TAZ, we calculated the unit area CO₂ emissions (kg/km²) for each TAZ over 24 h, defined as the CO₂ emissions (CE) per TAZ. We then used the CE results to construct feature vectors for clustering. First, the CE for each TAZ in a 24-h period is calculated, and a 24-dimensional vector V_i = {CE_1i,CE_2i,…,CE_24i} is constructed to denote the emission characteristics for each TAZ across a 24-h period. Second, the number of clusters is determined, and the cluster count is verified with validity metrics. In this study, the elbow method is employed, with the sum of squared error (SSE) as the criterion for evaluating the clustering results. For number of clusters ranges from 1 to 10, as the number of clusters increases to 4, the SSE gradually stabilizes and decreases, as shown in Figure 10a. Therefore, it is established that the number of clusters is 4. Finally, the k-means clustering is applied to partition each 24-dimensional vector into different clusters based on the CO₂ emission patterns of the TAZs. This approach is used to reveal the distribution of TAZs exhibiting similar temporal and spatial GHG emission characteristics.

Figure 10b shows the clustering results. Among these results, there are 54 TAZs in cluster 1 and 141 TAZs in cluster 2, which represent low-emission areas. Cluster 1 is distributed in subcentral urban areas, while cluster 2 is distributed on the outskirts of the city. On the other hand, cluster 3 and cluster 4 are high-emission areas, with 40 TAZs in cluster 3 and 33 TAZs in cluster 4. Cluster 3 is primarily located in the city center, forming a strip-like distribution from south to north, while cluster 4 is mainly found in newly developed urban center areas. Figure 10c shows the distribution characteristics of CO₂ emissions for all TAZs in clusters 1–4. It can be observed that TAZs with similar emission characteristics are grouped into the same cluster. For TAZs with similar spatiotemporal GHG emission characteristics after clustering, we apply the method in Section 4.2 to analyze the impact factors on weekdays versus weekends for different clusters, to identify the important factors that result in differences, and to analyze the influencing mechanisms.

5.4. Nonlinear Impact Analysis of BE on GHG Emissions from Taxis

Based on the clustering results from Section 5.3, utilizing the independent variables described in Section 3.2 and following the principles and parameters of the GBDT model outlined in Section 4.2, we distinguish between weekdays and weekends to construct the GBDT model and analyze the nonlinear effect of BE on the CEs of the TAZs in clusters 1–4.

Table 5. FIS ranks of independent variables in the weekday and weekend models.

Feature	Weekday								Weekend
	Cluster 1	Ranking	Cluster 2	Ranking	Cluster 3	Ranking	Cluster 4	Ranking	Cluster 1	Ranking	Cluster 2	Ranking	Cluster 3	Ranking	Cluster 4	Ranking
Intersection density	0.027	13	0.079	4	0.065	5	0.064	6	0.028	8	0.003	8	0.041	9	0.085	6
Road network density	0.044	6	0.130	3	0.087	4	0.230	1	0.091	4	0.078	4	0.087	4	0.139	1
Mixed degree of urban functions	0.067	5	0.029	8	0.034	12	0.072	4	0.042	6	0	9	0.018	16	0.078	7
Restaurant density	0.040	10	0.076	5	0.043	11	0.028	12	0.006	16	0.102	3	0.069	5	0.049	9
Scenic area density	0.027	12	0.007	15	0.024	15	0.064	5	0.010	14	0	9	0.028	14	0.098	3
Communal facility density	0.041	9	0.013	12	0.048	9	0.041	9	0.021	10	0	9	0.042	8	0.047	10
Shopping use density	0.031	11	0.026	9	0.056	7	0.017	16	0.020	11	0	9	0.130	3	0.030	12
Metro/bus/taxi stops density	0.244	1	0.252	1	0.137	2	0.027	13	0.379	1	0.350	2	0.192	1	0.095	4
Education facility density	0.042	8	0.011	13	0.169	1	0.025	14	0.025	9	0	9	0.042	7	0.010	16
Finance use density	0.102	3	0.005	16	0.062	6	0.024	15	0.146	2	0	9	0.039	11	0.010	15
Sport use density	0.042	7	0.022	10	0.051	8	0.033	11	0.028	7	0	9	0.047	6	0.027	13
Service of life density	0.020	15	0.018	11	0.028	14	0.039	10	0.013	13	0.004	7	0.036	12	0.124	2
Residential density	0.159	2	0.231	2	0.103	3	0.126	2	0.122	3	0.041	1	0.140	2	0.019	14
Accommodation service density	0.072	4	0.062	6	0.031	13	0.063	7	0.046	5	0.048	5	0.030	13	0.090	5
Government use density	0.018	16	0.032	7	0.015	16	0.094	3	0.008	15	0.005	6	0.020	15	0.055	8
Medical Institutions density	0.024	14	0.007	14	0.044	10	0.055	8	0.016	12	0	9	0.040	10	0.043	11

presents the feature-importance rankings of the independent variables for CEs in the four clusters on weekdays and weekends. A higher numerical FIS of the independent variable indicates that this variable is more important in the model and has a greater impact on taxi CEs. Although we obtained rankings and FISs for all 16 variables, due to the large number of independent variables and the small impact of some of them on the dependent variables, we refer to the study of Yang et al. [59] and select the four independent variables with the highest FISs for each ranking as the important impact factors in this study.

Overall, the important impact factors for CEs on weekdays and weekends in cluster 1–3 are generally the same but exhibit slight differences. Cluster 4 has different important impact factors on weekdays and weekends, except for road network density, which ranks first on both weekdays and weekends in terms of FIS. Metro/bus/taxi stops density, residential density, and road network density appear most frequently as important impact factors, consistent with the conclusions of Yang [38]. An et al. and Guo et al. [60,61] also reported that the locations of public transport stops in a transport network have an important influence on the travel of residents. Thus, establishing rational distributions of the residential population and public transport stops is important for reducing urban GHG emissions.

In cluster 1, for both the weekday and weekend models, the metro/bus/taxi stops density is an important impact factor affecting the CEs of taxis. The reason for this is that travel demand and attraction are generally greater in areas with higher metro/bus/taxi stops densities and residential densities, and cruising taxi drivers actively travel to these areas to find passengers when their taxis are unoccupied, which results in greater CEs in these areas. The finance use density has a slightly greater influence on weekends than on weekdays. The only different important impact factor in two models on cluster 1 is the accommodation service density (rank 4, 0.072) of the weekday model versus the road network density (rank 4, 0.091) of the weekend model. Notably, due to the heavier traffic congestion on weekends than on weekdays, the road network density has a greater effect on the CEs of taxis. Cluster 2 is located on the outskirts of the city. In both the weekday and weekend models, metro/bus/taxi stops density and residential density consistently rank among the top two factors in terms of FIS values. Especially on weekends, the impact of residential density reaches 0.410. Road network density in cluster 2 is an important impact factor on both weekdays and weekends, and the impact of this factor is greater on weekdays (rank 3, 0.130) than on weekends (rank 4, 0.078). This finding is slightly different from the findings for cluster 1. In addition, on weekdays, intersection density (rank 4, 0.079) is one of the important impact factors impacting the CEs in cluster 2, and the other important impact factor on weekends is restaurant density (rank 3, 0.102), suggesting that the restaurant sector in these areas has a greater attraction to taxi travels on weekends. Notably, the FIS of the independent variables such as shopping use density and education facility density in cluster 2 on weekends is 0, indicating that on weekends, these independent variables do not have a direct impact on the CEs of taxis in cluster 2. This is also related to the fact that the areas in cluster 2 are far from the city center and not attractive for taxi trips on weekends. Education facility density (rank 1, 0.169) is the most important impact factor in the weekday model of cluster 3. Due to the high pedestrian and vehicular traffic near education facilities on weekdays, there is a high demand for transportation, and this scenario is prone to causing traffic congestion, resulting in slow-moving taxi operations. Hence, it becomes a critical factor impacting the CE in these areas. Road network density has an FIS of 0.087 for both weekdays and weekends, indicating that it has the same important influence on the dependent variable in cluster 3 in both the weekday and weekend models. On weekends, shopping use density (rank 3, 0.130) emerges as an important impact factor. Due to increased leisure and recreational activities among residents during weekends, areas with high shopping use density experience elevated demand and attraction for taxi travel, leading to an increase in CEs. The top four independent variables in cluster 4 for weekday impact factors are road network density (rank 1, 0.230), residential density (rank 2, 0.126), government use density (rank 3, 0.094), and mixed degree of urban functions (0.072). The top four independent variables for weekend impact factors are road network density (rank 1, 0.139), service of life density (rank 2, 0.124), scenic spot density (rank 3, 0.098), and metro/bus/taxi stops density (rank 4, 0.095). The influence of road network density is greater on weekdays than on weekends, and except for road network density, the remaining three important impact factors are different between weekdays and weekends in cluster 4.

Next, we visualize the nonlinear impacts of the important impact factors on taxi CEs in TAZs on both weekdays and weekends using partial dependence plots (PDPs). Assuming that the other variables remain constant, these plots illustrate the impact mechanisms of individually important impact factors in relation to the CEs of taxis, as depicted in Figure 11, Figure 12, Figure 13 and Figure 14. It is evident from the plots that variables with high FIS exhibit pronounced nonlinear relationships with taxi CEs in TAZs.

Figure 11, Figure 12, Figure 13 and Figure 14 show the PDPs for clusters 1–4 for weekdays and weekends considering the nonlinear associations between the important impact factors and taxi CEs. A significant positive or negative nonlinear relationship between important impact factors and the dependent variable can be found. Despite variations in FISs, the impact mechanisms of the same important impact factors on the dependent variable are largely consistent between the weekday and weekend models within the same cluster; the trends of the PDPs remain largely similar. Examples include road network density in cluster 1, finance use density in cluster 1, and residential density in cluster 2. However, it should be noted that in some clusters, there are disparities in the PDP trends for the same important impact factors between weekdays and weekends, such as for road network density in cluster 2. Additionally, there is a threshold for the impact of the independent variable on the dependent variable, i.e., when a certain important impact factor reaches the threshold, the change in the value of the dependent variable grows rapidly or stagnates. In the next part of this study, we analyze in detail the impact threshold of the independent variable on the dependent variable.

As shown in Figure 11, in cluster 1, the trends of the effects of metro/bus/taxi stops density, finance use density, and residential density on CEs are generally consistent between the weekday and weekend models. When the metro/bus/taxi stops density is less than 8 (Figure 11a,e), there is no significant impact on CEs. However, as the metro/bus/taxi stops density increases from 8 to 12, the intensity of taxi CEs rapidly increases, stabilizing at approximately 420 between 12 and 30, and it then rapidly increases and stabilizes for metro/bus/taxi stops densities greater than 30. The same conclusion is reached for two important impact factors, finance use density and residential density, which rapidly increase from 2–7 and 3–5, respectively, after which a threshold is reached and no further growth is observed. According to the weekend model of cluster 1, road network density has a negative impact on the CEs of taxis (Figure 11h), similar to the findings of Wang et al. [44]. They argued that an increase in road network density enhances the accessibility and connectivity of an area, leading residents to choose nonmotorized transportation, thereby reducing CO₂ emissions. Figure 11h shows that when the road network density increases to 12.5–15, taxi CEs decrease significantly; when the road network density is higher than 15, taxi CEs no longer decrease. The accommodation service density in the weekday model also has a nonlinear negative effect on the taxi CEs in cluster 1 (Figure 11d).

In the weekday and weekend models of cluster 2 (Figure 12), an increase in residential density yields a rapid increase in taxi CEs within the 0–2.5 range, after which emissions stabilize in the 2.6–6 range (Figure 12b,e). This phenomenon is attributed to certain residential areas prohibiting taxi entry. When the residential density within a TAZ exceeds 2.5, taxis are constrained to picking up, dropping off, and searching for passengers on the outer roads of the residential area, preventing further increases in CEs. Similarly, in the weekend model of cluster 2, when the restaurant density and road network density are within the ranges of 0–10 and 4.5–7, respectively, rapid increases in the CEs of taxis occur. Subsequently, further increases in these variables no longer impact CEs. The influence of metro/bus/taxi stops density on taxi CEs exhibits roughly the same trend in the weekday and weekend models. Notably, for cluster 2, unlike the negative effect for cluster 1, road network density has a positive influence on the dependent variables (Figure 12c,h), which is in accordance with the results of Shu and Lam [45], who found that an increase in road network density leads to an increase in CO₂ emissions by increasing the number of vehicles entering the area, which increases the congestion of the road network. The famous “Braess Paradox” [62] also supports this conclusion. In the weekday model, the effect of another important impact factor, intersection density, on the dependent variable shows an inverted U-shaped trend (Figure 12d), as dense intersections increase road connectivity. When the density of intersections is less than 15, an increase in intersection density increases traffic flow, and CE increase. When the intersection density is higher than 15, the connectivity of the road network increases, and more people choose nonmotorized travel, which leads to a gradual decrease in CO₂ emissions. The calculated threshold value of the effect of the intersection density on CEs is equal to that obtained in the study of Guo et al. [63].

In cluster 3, as the metro/bus/taxi stops density increases from 0 to 35, the dependent variable rapidly increases from 200 to 700. Further increases beyond this range no longer affect the variable (Figure 13b,e). A similar situation is observed for the impact of residential density on CEs (Figure 13c,f). According to the weekday model of cluster 3, as the education facility density increases from 0 to 22, taxi CEs decrease from 410 to 340 (Figure 13a). This is because cluster 3 is located in the city center and includes numerous educational facilities. In these areas, residents may choose to walk or use nonmotorized transportation for short commutes on weekdays. Additionally, parents often use private cars for student drop-off and pick-up, resulting in a lower demand for taxi services in these areas on weekdays. In addition, there is a threshold for the effect of education facility density on CEs, and the dependent variable no longer decreases when it is greater than 21. Road network density has a different effect on CEs in the weekday model than in the weekend model, and in the weekend model, when the road network density is greater than 10, taxi CEs tend to stabilize and no longer increase (Figure 13f). In contrast, in the weekday model, emissions continue to fluctuate and grow after the road network density surpasses 10 (Figure 13d). These variation trends are similar to the nonlinear effects of changes in the road network density on CO₂ emissions from online-hailing taxis, as observed in the study by Gao et al. [7]. Moreover, in the weekend model, shopping use density exhibits a stepwise increase, rapidly rising within the intervals of 60–90 and 200–230.

In cluster 4, the changes in the influence of the road network density on the values of dependent variable in the weekday model and the weekend model are different, but the trends are the same, displaying growth followed by a decrease and eventual plateau (Figure 14a,e). When the road network density exceeds 15, the road capacity reaches saturation, and the traffic volume in the area no longer increases. Therefore, the CEs do not increase further. Cluster 4 is located in newly developed urban center areas, and in the weekday model, the increase of government use density and residential density makes the taxi CEs display stepwise growth (Figure 14b,c). When the mixed degree of urban functions increases from 0.7 to 1.85, the taxi CEs increase, especially with rapid growth in the interval of 1.75–1.85 (Figure 14d). When the mixed degree of urban functions increases further, the taxi CEs start to decrease, indicating that taxi CEs can be reduced when the mixed degree of urban functions reaches a certain threshold. This proves that in newly developed urban center areas, when the mixed degree of urban functions reaches a certain threshold, it can reduce the travelling demand of taxis, thus reducing the emissions. This conclusion is similar to the findings of Wu et al. [53], who concluded that there are upper and lower thresholds for land use mixing, suggesting that some of the newly developed urban center areas in a city need to integrate urban functions to control GHG emissions from urban transportation. Regarding the weekend model, the CEs experience rapid growth when the service of life density reaches the range of 80–120 (Figure 14f). Additionally, because Xi’an is a renowned tourist city, the scenic area density serves as an important influencing factor in the weekend model. CEs rapidly increase when the scenic area density is within the range of 0–2, but beyond this threshold, the CEs remain unchanged (Figure 14g). Notably, for cluster 4 in the weekend model, there is a complementary association between public transportation and taxi demand. As the metro/bus/taxi stops density increases, residents in this cluster tend to choose regular bus and subway services more frequently on weekends, leading to a decrease in taxi CEs (Figure 14h). This observation suggests that in newly developed urban center areas, appropriately increasing the density of public transportation stations, such as bus and subway stops, can contribute to reducing GHG emissions from taxis.

6. Conclusions and Recommendations

In this study, we utilize taxi trajectory data and the COPERT model to estimate the emissions of four GHGs in the research area. Subsequently, TAZs with similar CO₂ emission patterns are clustered using the k-means algorithm, followed by an in-depth analysis from the perspective of the BE to explore the causes of spatiotemporal variations in taxi GHG emissions across different areas. Finally, we select 16 variables of the BE based on “5Ds” elements and employ the GBDT model to construct weekday and weekend models to analyze the factors that influence taxi CEs for different emission pattern clusters. We identify the top four most impactful factors with the highest FISs for each cluster in both weekday and weekend models and analyze the nonlinear relationships between these important impact factors and taxi CEs.

The main conclusions and policy recommendations are elaborated as follows: (1) The total GHG emissions from taxis in Xi’an exhibit higher levels on weekdays compared to weekends, with the highest CO₂ unit emission factors occurring during P.M. peak hours and showing significant differences between weekdays and weekends. (2) Metro/bus/taxi stops density, residential density, and road network density are the most common important impact factors in the TAZs. (3) The independent variables with high FISs have a significant nonlinear relationship with taxi CEs, and there is a threshold for the influence of each important impact factor on the dependent variable. For cluster 2 (outskirts of the city) in the weekday model, the impact of intersection density on taxi CEs exhibits an inverted U-shaped trend. The CEs of taxis reach a maximum threshold when the intersection density exceeds 15. (4) Despite different FIS values for the same important impact factors across weekday and weekend models within clusters, the PDP trends remain consistent, indicating that the same independent variables influence taxi CEs similarly on both weekdays and weekends, such as residential density in cluster 1 and road network density in cluster 4. (5) The influence of the independent variables on taxi CEs varies across clusters. For cluster 1 (subcentral urban area) in the weekend model, an increase in road network density has a negative effect on taxi CEs, whereas in city center areas (cluster 3) in the weekend model, this effect is positive. In newly developed urban center areas (cluster 4), as the mixed degree of urban functions increases, the taxi CEs show a trend of initially increasing and then declining.

The above research findings indicate that the association between the BE and taxi GHG emissions is complex. The impact of BE elements on taxi GHG emissions varies across different functional zones within a city. Increases in the road network density in subcentral urban areas, the metro/bus/taxi stops density and education facility density in city centers, and the mixed degree of urban functions in newly developed urban center areas to 1.85 or above contribute to reductions in GHG emissions from taxis to a certain extent. Since taxi operations reflect the characteristics of urban transport and residents’ travel, the results indicate that the reduction of urban transport GHG emissions should begin with optimizing the urban structure, enhancing the mixed degree of urban functions, and adjusting the functionality of urban areas. These findings are crucial for developing comprehensive emission-reduction strategies within intelligent transportation systems.

This study has certain limitations that require further exploration. First, our study does not obtain data information on hybrid electric vehicles (HEVs) and does not account for scenarios where taxi trajectories pass through TAZ boundaries, which may lead to discrepancies between the calculated results and real-time data from environmental monitoring stations. Second, we use the GBDT model to investigate the nonlinear characterization of a single independent variable in relation to the dependent variable, ignoring possible correlations among the independent variables, which could be further investigated in the future by examining the joint impacts of two or three independent variables on the dependent variable. Finally, future studies should aim to conduct a finer-scale study of GHG emissions, considering HEVs from traffic flows at the road network level to provide a more accurate basis for analyses and identify strategies for urban GHG emission reduction.