**E**ff**ect of Landscape Elements on the Symmetry and Variance of the Spatial Distribution of Individual Birds within Foraging Flocks of Geese**

**Mads Bech-Hansen 1,\*,**†**, Rune M. Kallehauge 1,**†**, Dan Bruhn 1, Johan H. Funder Castenschiold 1, Jonas Beltoft Gehrlein 1, Bjarke Laubek 2, Lasse F. Jensen <sup>2</sup> and Cino Pertoldi 1,3**


Received: 23 July 2019; Accepted: 23 August 2019; Published: 2 September 2019

**Abstract:** Behavioural instability is a newly coined term used for measuring asymmetry of bilateral behavioural traits as indicators of genetic or environmental stress. However, this concept might also be useful for other types of data than bilateral traits. In this study, behavioural instability indices of expected behaviour were evaluated as an indicator for environmental stress through the application of aerial photos of foraging flocks of geese. It was presumed that geese would increase anti-predator behaviour through the dilution effect when foraging near the following landscape elements: wind turbines, hedgerows, and roads. On this presumption, it was hypothesized that behavioural instability of spatial distribution in flocks of geese could be used as indicators of environmental stress. Asymmetry in spatial distribution was measured for difference in flock density across various distances to disturbing landscape elements through the following indices; behavioural instability of symmetry and behavioural instability of variance. The behavioural instability indices showed clear tendencies for changes in flock density and variance of flock density for geese foraging near wind turbines, hedgerows, and roads indicating increasing environmental stress levels. Thus, behavioural instability has proven to be a useful tool for monitoring environmental stress that does not need bilateral traits to estimate instability but can be applied for indices of expected behaviour.

**Keywords:** environmental stress; behavioural instability; biomonitoring tool; unmanned aerial vehicles (UAVs); disturbance

#### **1. Introduction**

#### *1.1. Behavioural Instability*

A study by Pertoldi et al. [1] coined the term behavioural instability as an indicator of genetic or environmental stress based on asymmetry in a bilateral behavioural trait, e.g., clockwise and counter-clockwise directional movement. The authors applied the conventional indices used for the estimation of developmental instability in directional movements. However, the concept of behavioural instability can be considered in a broader sense [1], and the concept might prove to be a viable biomonitoring tool in other fields of research, e.g., studies on environmental impact on wildlife. Furthermore, behavioural instability might be useful for other types of data with an expectation of symmetry other than bilateral traits.

#### *1.2. Environmental Stressors of Geese*

An important factor influencing the habitat use of water birds is the potential disturbance by landscape elements such as wind turbines [2], hedgerows, and roads [3,4]. Disturbance might be caused by anthropogenic activity or these landscape elements offering hiding spots for predators [5–7], which essentially can result in displacement and habitat loss. Presuming geese foraging near disturbing landscape elements experience increased stress levels as a result of the potential increased predation risk, anti-predator behaviour through the dilution effect would increase [8]. The dilution effect is defined as increasing group size or density as a strategy for reducing individual risk of being targeted by a predator at the cost of reduced foraging efficiency, e.g., due to increased competition [8].

#### *1.3. Behavioural Instability of Symmetry*

Based on the presumption of geese increasingly relying on the dilution effect when foraging near disturbing landscape elements, it would be possible to analyse stress as a function of flock densities at various distances to these landscape elements. Thus, it would be possible to show an asymmetrical distribution by plotting the density of the birds versus the distance to an obstacle in a linear regression as illustrated in Figure 1b. If the obstacle is inducing anti-predator behaviour, the slope (a) of the linear regression will become negative (a < 0) indicating a decreasing flock density with distance. Deviations from a symmetrical distribution (a - 0) will be considered as an estimator of behavioural instability and will be referred to as behavioural instability of symmetry (BSYM).

#### *1.4. Behavioural Instability of Variance*

Pertoldi et al. [1] also noted that different genotypes have different perceptions of stressors, and that in a suboptimal environment the perception of stress will become more differentiated among genotypes [1]. If this finding is valid for the way in which the birds behave in the presence of stress, then heterogeneity can be expected where the average distance of the residuals from the regression line is not constant at different distances from the landscape element as illustrated in Figure 1c. This variation can be considered as another estimator of behavioural instability, with a higher variance of the residuals being interpreted as a smaller capacity to predict the behaviour of the individuals in the presence of stressors. This estimator is later referred to as deviation from a behavioural instability of variance (BVAR).

#### *1.5. Aims of the Investigation*

Based on the mentioned theories, we hypothesized that behavioural instability indices of expected behaviour can be used as a tool for measuring environmental stresses in wildlife. This will be evaluated through the analysis of the behavioural instability of the spatial distribution of flocks of geese relying on two different measures: behavioural instability of symmetry (BSYM) (Figure 1b) and behavioural instability of variance (BVAR) (Figure 1c). These indices were tested using aerial photos of flocks of geese foraging (Figure 1b) near the following landscape elements; wind turbines, hedgerows, and roads, which are all known to cause the disturbance of geese.

**Figure 1.** Overview of behavioural instability indices used to monitor one flock of pink-footed geese foraging near a road: (**a**) Aerial photos of a flock of individual geese are georeferenced in geographic information systems (GIS) and geographic positions of individual geese are geotagged (illustrated as circles). Flock density is quantified by Kernel density estimates using the geographic positions of identified birds. Kernel density estimates are illustrated as white (low density) to dark red (high density); (**b**) Density estimates of individuals as a function of distance to the nearest landscape element, in this case, roads. The mean distance is illustrated as a red vertical line. A linear regression of density (LRD): Density Estimate = a × Distance + b (illustrated as a grey line) is fitted to measure behavioural instability through asymmetry of density between distances left and right of the mean distance. Symmetrical density would result in a non-significant slope (a = 0). Contrarily, an anti-predator behaviour induced by the road would yield a negative slope (a < 0). Residual distances between observed densities and LRD is illustrated as blue lines; (**c**) Residual distances of density (Figure 1b) as a function of distance to the nearest landscape elements, in this case roads. A linear regression of residuals (LRR): Residual = a × Distance + b (illustrated as a grey line) is fitted to measure behavioural instability through asymmetry of variance in density between distances left and right of the mean distance. A constant variance in density across distance would result in a non-significant slope (a = 0). Contrarily, increased variance would yield a negative slope (a < 0).

#### **2. Materials and Methods**

Aerial photos of foraging geese were used to evaluate changes in anti-predator behaviour by examining the spatial distribution of flocks of different distances to wind turbines, hedgerows, and roads. This was done to evaluate the two indices: behavioural instability of symmetry (BSYM), and behavioural instability of variance (BVAR), as indicators of environmental stressors. The studied species were barnacle goose (*Branta leucopsis*) (19 flocks, with a total of 18,925 individuals), pink-footed goose (*Anser brachyrhynchus*) (23 flocks, with a total of 26,313 individuals), and greylag goose (*Anser anser*) (5 flocks, with a total of 353 individuals).

Spatial distribution was measured based on distances of individuals to said landscape elements, as well as from flock density estimates extracted from geographic information systems (GIS) (Figure 1a). Thus, differences in spatial distribution will be tested between flocks of geese are observed at different distances from the different landscape elements, through the analysis of behavioural instability of symmetry (BSYM) (Figure 1b) and behavioural instability of variance (BVAR) (Figure 1c).

#### *2.1. Data Collection*

In the period March 2017 to February 2018 aerial photos of foraging flocks of geese were collected in Northern Jutland, Denmark, mainly around Klim and Nørrekær Enge wind farm, using a UAV of the model DJI Phantom 4 Pro Quadcopter. Data were collected during daytime on 26 different days, chosen based on the season (winter migration of geese) and weather forecast (no rain and low wind speed). All UAV overflights were performed at an altitude of 100 m either flown manually using the DJI GO 4 Drone application (Version 4.0.6), with the video camera capturing aerial photos vertical downwards in 4K resolution or flown automatically using DroneDeploy (Version 2.69). Additionally, data collection using the UAV were conducted following the recommended flight altitude (100 m) and take-off distance from the studied flocks (~500 m) suggested by Bech-Hansen et al. [9] as a means to prevent initial disturbance of the birds.

#### *2.2. Data Extraction*

Aerial photos were imported into QGIS (Version 2.8.20) using the geo-referencing plugin GDAL (Version 3.1.9) as large ortho-mosaics created using Autostitch (Demo version) or DroneDeploy limited free services. Birds were then identified, given UTM coordinates, and, if possible, given a species-specific id. A total of 45,591 birds were identified from 47 flocks of geese. Aerial photos were not corrected for barrel distortion effects as it was considered of minor influence.

For each flock, a density heatmap was produced in QGIS from the location of identified birds using the Heatmap plugin (Version 0.2), based on the Kernel density estimation [10], with radius set to 5 m, and cell size x and y set to 1 m each, thus, resulting in comparable density estimates across all studied flocks. Density estimates for each identified bird were extracted using the QGIS plugin Point Sampling Tool (Version 0.4.1), while the distance between the identified birds and nearest wind turbine, hedgerow, and road, were measured using the NNJoin plugin (Version 1.3.1) in QGIS.

#### *2.3. Data Analysis*

All analyses were conducted based on bird distances to the three landscape elements: wind turbines, hedgerows, and roads. Based on previous studies, birds at distances of more than 600 m from the nearest associated landscape structure were assumed not to be affected by these structures [2–4,6] and were therefore excluded from the analysis.

Density estimates of all individual birds were plotted as a function of bird distances to nearest landscape element and fitted with a linear regression (Figure 1b), which is, henceforth, referred to as linear regression of density (LRD) [Equation (1)].

$$\text{Density Estimate} = \mathbf{a} \times \text{Distance} + \mathbf{b} \tag{1}$$

An LRD with negative significant slope (a < 0) indicates an increase in flock density in the direction of the studied landscape element and vice versa. A negative slope also indicates an asymmetrical distribution of the flock density and consequently an increased BSYM in the direction of the studied landscape element.

The absolute values of the residual distances of the LRD [Equation (2)] as previously mentioned, were plotted versus the distance to the obstacles in a linear regression of residuals (LRR) [Equation (3)] (Figure 1c).

$$\text{Residual} = \text{Observed density} - \text{predicted density} \,\text{(LRD)}\tag{2}$$

$$\mathbf{b} \text{ Residual} = \mathbf{a} \times \mathbf{Distance} + \mathbf{b} \tag{3}$$

A slope of LRR different from 0 (a - 0) indicates heterogeneity of the residuals. A negative slope (a < 0) indicates an increasing variance of the residuals' absolute distance from LRD in the direction of the studied landscape element and vice versa. Therefore, an increased variance can be considered as an increased BVAR in the direction of the studied landscape element and vice versa.

The determination coefficients (r2), the slope and the intercept of all LRDs and LRRs have been estimated as well as significance for each slope (later noted as asterisks; \*: *p* < 0.05, \*\*: *p* < 0.01, \*\*\*: *p* < 0.001).

#### **3. Results**

#### *3.1. Correlation between Distance to Landscape Elements and Behavioural Instability of Symmetry (BSYM)*

#### 3.1.1. Wind Turbines

The barnacle geese and pink-footed geese showed a significant increase of BSYM with decreasing distance to the wind turbines (r<sup>2</sup> = 20.53%) \*\*\* and (r<sup>2</sup> = 10.51%) \*\*\* respectively (Table 1 and Figure 2a,d). The greylag geese showed the same trend, but the regression was not significant (r2 = 0.45%, *p* > 0.05) (Table 1 and Figure 2g).

#### 3.1.2. Roads

The greylag geese showed a significant increase of BSYM with decreasing distance to the roads (r2 = 8.08%) \*\*\* (Table 1 and Figure 2h). Whereas, barnacle geese and pink-footed geese showed (with r<sup>2</sup> < 5%) the opposite trend with a significant decrease of BSYM with decreasing distance to the roads (r2 = 0.04%) \*\*\* and (r<sup>2</sup> = 1.01%) \*\*\* respectively (Table 1 and Figure 2b,e).

#### 3.1.3. Hedgerows

Barnacle geese and pink-footed geese showed a significant (with r<sup>2</sup> < 5%) increase of BSYM with decreasing distance to the hedgerows (r<sup>2</sup> = 1.46%) \*\*\* and (r<sup>2</sup> = 2.13%) \*\*\* respectively (Table 1 and Figure 2c,f).

The greylag geese showed (with a non-significant regression and r<sup>2</sup> < 5%) the opposite trend with a significant decrease of BSYM with decreasing distance to the hedgerows (r2 = 0.18%, *p* > 0.05) (Table 1 and Figure 2i).

#### *3.2. Correlation between Distance to Landscape Elements and Behavioural Instability of Variance (BVAR)*

#### 3.2.1. Wind Turbines

The barnacle geese and pink-footed geese showed a significant increase of BVAR with decreasing distance to the wind turbines (r<sup>2</sup> = 14.69%) \*\*\* and (r<sup>2</sup> = 22.18%) \*\*\* respectively (Table 1 and Figure 3a,d). Whereas, the greylag geese showed (with r2 < 5% and *p* > 0.05) the opposite trend with a decrease of BVAR with decreasing distance to the wind turbines (Table 1 and Figure 3g).

**Table 1.** Species: Barnacle goose (BG), pink-footed goose (PINK), and greylag goose (GREY); distances from the landscape elements: wind turbines (dwi), roads (dro), and hedgerows (dhe), n = number of measurements. Coe fficient of determination (r2) slope (a), intercept (b) and significance level (*p*) of both linear regression of density (LRD) (which regress the density estimates of individual birds as a function of bird distances to nearest landscape element) and of linear regression of residuals (LRR) (which regress the absolute values of residual distance from LRD as a function of bird distances to nearest landscape element). The r2 values above 5% (r2 > 5%) are in bold.

39

of bird distances to the nearest landscape element. The following species were regressed: barnacle goose (BG), pink-footed goose (PINK), and greylag goose (GREY). Distances from the landscape elements; wind turbines (dwi), roads (dro), and hedgerows (dhe). The r2 values of LRR above 5% (r2 > 5%) are in bold.

#### 3.2.2. Roads

The greylag geese showed a significant increase of BVAR with decreasing distance to the roads (r2 = 10.67%) \*\*\* (Table 1 and Figure 3h). Whereas, barnacle geese and pink-footed geese showed a significant increase (with r<sup>2</sup> < 5%) increase of BVAR with decreasing distance to the roads (r2 = 2.35%) \*\*\* and (r<sup>2</sup> = 2.64%) \*\*\* respectively (Table 1 and Figure 3b,e).

#### 3.2.3. Hedgerows

The barnacle geese showed a significant increase of BVAR with decreasing distance to the hedgerows (r<sup>2</sup> = 5.20%) \*\*\*, (Table 1 and Figure 3c). Whereas, pink-footed geese showed a significant (with r<sup>2</sup> < 5%) increase of BVAR with decreasing distance to the hedgerows (r2 = 2.74%) \*\*\*, (Table 1 and Figure 3f).

Lastly, the greylag geese showed (with a non-significant regression and r2 < 5%) an increase of BVAR with decreasing distance to the hedgerows (r2 = 0.46%, *p* > 0.05) \*\*\*, (Table 1 and Figure 3f).

#### **4. Discussion**

We have shown that the spatial distribution of foraging geese changes when they get closer to a disturbing landscape element (Table 1, Figures 2 and 3) and clear tendencies were observed indicating an increasing anti-predator cohesion of flocks of geese with decreasing distance to disturbing landscape elements. The change in spatial distribution is visible both for the behavioural instability of symmetry (BSYM) and behavioural instability of variance (BVAR), which will be discussed below.

There are clear tendencies, with few exceptions (with r2 < 5%), for negative slopes of the linear regression of density (LRD), which indicates an increasing BSYM with decreasing distance to the landscape elements (Table 1 and Figure 2). There is also a clear tendency with few exceptions (with r<sup>2</sup> < 5%) for a negative slope of the linear regression of residuals (LRRs), which means an increasing BVAR with decreasing distance to the landscape elements (Table 1 and Figure 3). Thus, both indices show asymmetry of spatial distribution, implying environmental stress in flocks of geese induced by foraging near the studied landscape elements, which indicates these measurements to be useful tools for monitoring environmental stress. We have chosen the 5% threshold arbitrarily; however, it is also notable that the same negative trends were observed for other LRD and LRR although with r2 below 5% (Figures 2 and 3). However, such relatively weak trends might not be of biological importance and should thus be interpreted with caution. Additionally, *p*-values also should be evaluated cautiously as the large sample size increases the significance of the indices even with low r2 as seen in both indices (Table 1). These relatively low r<sup>2</sup> values might be caused by noise from other factors influencing the density of bird flocks, masking the disturbing effects of the landscape elements. Such factors might include flock size [11] as well as food and water distribution [12], which are both known to affect the density and spatial distribution of bird flocks. Especially BVAR might have been affected by variations in flock size as flocks of different size prioritise anti-predation behaviour and foraging beneficial behaviour differently. A study by Lazarus [11], who examined the influence of flock size on the vigilance of white-fronted geese (*Anser albifrons*), noted that the percentage of vigilant birds in a flock would decrease steeply at flock sizes above 200–300 individuals [11]. Hence, larger flocks might prioritise foraging beneficial behaviour while smaller flocks prioritise anti-predator behaviour through the dilution effect [7,11,13,14]. However, trends were still observed for both BSYM and BVAR that prove that both methods can be applied to wildlife behaviour with multiple random influences, as long as the factor of interest is measured across a correlated variable.

In our study linear relationships were assumed and therefore only linear regression has been utilized. However, BSYM and BVAR could also be applied in studies with a non-linear relationship between the independent value and the dependent value. In this case, the asymmetry on the left and right side of mean can vary and might be estimated by the first derivative, which is equivalent to the slope observed for a certain value of the dependent value. Therefore, the BSYM index can be considered the equivalent of the absolute values of the difference between right and left value of a trait, which was referred to as fluctuating asymmetry (FA1) [15] in a study by Palmer and Strobeck [15]. Whereas, BVAR indices can be considered the equivalent of the phenotypic variability of trait (Vp). For both indices, FA1 and Vp are often utilised as estimators of developmental instability (DI) [16].

#### **5. Conclusions**

We have tested both indices against distance to landscape elements as an independent variable; however, the indices could also be utilised for other possible variables, e.g., regression against time or temperature. The merit of BSYM and BVAR is that they do not need bilateral traits for the estimation of the instability, which in our case we have called behavioural instability as we are estimating a deviation from expected behaviour, which is a constant distance between every single individual due to a uniform distribution of the individuals. These indices have proven effective tools for assessing the effect of environmental stress on expected behaviour. Thus, the term behavioural instability has been proven useful for data other than bilateral data as discussed in Pertoldi et al. [1]. We expect that the application of these two indices BSYM and BVAR could have several applications in applied ecology as both indices can be regressed against environmental gradients or temporal series.

**Author Contributions:** M.B.-H., R.M.K. and C.P. conceptualised the project and did the methodology. C.P. provided the resources, while M.B.-H., R.M.K., J.H.F.C. and J.B.G. did the investigation and data curation. M.B.-H., R.M.K., D.B., J.B.G. and C.P. did the formal analysis and visualization. M.B.-H., R.M.K. and C.P. took the lead in writing the manuscript. All authors provided critical feedback and helped shape the research, analysis and manuscript. M.B.-H. and R.M.K. were responsible for the project administration as well as, research activities and planning, who were supervised by D.B., B.L., L.F.J. and C.P.

**Funding:** Mads Bech-Hansen, Rune M. Kallehauge and Cino Pertoldi were financed by the Aalborg Zoo Conservation Foundation (AZCF: grant number 2018-4).

**Acknowledgments:** We would like to thank the Lead Guest Editor and the anonymous reviewers for invaluable suggestions and comments on the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Evolutionary Game Research on Symmetry of Workers' Behavior in Coal Mine Enterprises**

#### **Kai Yu 1, Lujie Zhou 1,2,\*, Qinggui Cao 1,\* and Zhen Li <sup>1</sup>**


Received: 4 January 2019; Accepted: 28 January 2019; Published: 31 January 2019

**Abstract:** Statistics show that humans' unsafe behaviors are the main cause of accidents. Because of the asymmetry of game benefits between managers and coal miners, the stability of workers' behaviors is affected and unsafe behaviors are produced. In this paper, the symmetry of the behavior benefits of coal mine workers is studied, using game theory. In order to observe the dynamic game evolution process of behavioral stability, the paper establishes a system dynamics (SD) model and simulates it. The SD simulation results show that with the continuation of the game, when the benefits for safety managers and workers are asymmetric and the safety manager's safety inspection benefits are less than the non-inspection benefits, the manager may not conduct safety inspections, which poses a great hidden danger to safety production. Through dynamic incentives to regulate the symmetry of income of coal mine safety managers and coal mine workers, the purpose of enhancing the stability of safety behavior is achieved. The research results of the paper have been successfully applied to coal mine enterprises.

**Keywords:** benefits symmetry; behavioral stability; evolutionary game; system dynamics; dynamic incentive

#### **1. Introduction**

A large number of statistics show that in some large-scale accidental disasters, unsafe human behavior is the main cause of the accident [1–3]. More than 80% of coal mine accidents in China are directly or indirectly caused by unsafe human behaviors [4]. In safety production, the game income of coal mine managers and workers is asymmetric, which leads to unstable behavior on both sides. Therefore, it is very important to study the symmetry of game benefits and behavior stability between managers and workers in coal mine enterprises.

In terms of coal mine accidents, Chen et al. studied the development trend of coal mine accidents in China, and analyzed the human factors, indicating that the unsafe behavior of human beings is the main cause of the accidents [5]. Nouri assessed the safety of the Iranian Gas Company and concluded that workers' unsafe behavior is the main cause of accidents [6]. Ianole studied applied behavioral economics and its trends [7]. Experts have shown from a wide range of studies the importance of unsafe human behavior to safety production.

If safety behaviors are studied in depth, we must first understand their specific definition. In the late 1970s, Heinrich proposed that people's unsafe behavior and material insecurity led to accidents [8]. Based on Heinrich, Reason proposed a model of human factors in the 1990s and established an accident

model, which was named "cheese". The theories presented that the potential risks in the system are deeply defended at all levels and ultimately lead to accidents [9].

Based on the definition of unsafe behavior, behavioral symmetry [10] and behavioral stability [11] were studied. In the construction industry, safety behaviors are studied from the aspects of safety cognition and safety climate. For example, Fang et al. established a cognitive model of unsafe behavior of construction workers [12]. Lyu et al. studied delay behavior from the aspects of safety climate and safety results [13]. Safety leadership has a great influence on the evolution of safety behaviors. For example, Beatriz et al. studied the influence of safety leadership on workers' safety behaviors in the light of working conditions [14], and Shen studied the transfer mechanism of safety behaviors [15].

In summary, most of the studies focus on the factors, formation mechanisms, and models of safety behaviors, while few focus on the dynamic evolution process of safety behaviors. Therefore, the paper uses dynamic evolutionary game theory to analyze the evolution process of safety behaviors.

Game theory has been widely used in many fields of scientific research. It discusses open research topics of importance to economics and the broader social sciences [16]. Game theory has been used to study the improvement of resources and environment. Zhao et al. analyzed environmental protection materials and environmental risks [17]. Using game theory, Feng et al. put forward the method of optimizing resource allocation [18]. In social science, Lu et al. established a multi-party evolutionary game model [19]. Sun et al. focused on the game between employer and worker mobility [20]. Wang et al. studied the interaction among government, enterprises, and workers [21].

Game theory has made outstanding contributions in social science and other research. Combining game theory with safety behaviors of coal miners provides a new research direction for studying the evolution process of them. If the game theory is used alone, the evolutionary results of safety behaviors in equilibrium state can be obtained, but the evolutionary state of safety behaviors in each time period cannot be observed better. Therefore, the paper combines game theory with system dynamics, gives full play to system dynamics, and deeply studies the changing state of workers' safety behaviors in each time period.

Professor Forrester (Jay W. Forrester) of the Massachusetts Institute of Technology (MIT) pioneered System Dynamics (SD) in 1956, which is combined qualitative analysis with quantitative analysis to study the functions of complex systems and the interaction of behaviors through model simulations [22,23]. In the construction industry, system dynamics is used to analyze the safety behaviors of workers, such as the feedback mechanism of workers' safety attitudes and safety behaviors [24], the influencing factors of workers' unsafe behaviors [25], and the influence of workers' work interference and family conflict on safety behaviors [26]. In enterprise management, system dynamics is used to simulate the relationship between safety investment and coal mine accidents [27], and optimize the risk management of chemical enterprises [28]. In other research fields, system dynamics is also used as a research tool to analyze the dependence between safety factors [29], and the safety psychological process of railway workers [30].

The behavior benefits matrix of the safety management department workers group is constructed to solve the problem of asymmetry of benefits and unstable behaviors between the safety management department and multiple workers groups. Based on the matrix, the evolutionary game analysis is carried out on the behavior stability of coal mine workers. The evolutionary game method is combined with system dynamics (SD), and Vensim is used to build the evolutionary SD model of safety behaviors. The dynamic evolution rules of the safety management department and the two coal mine workers groups are deeply analyzed through simulation by SD.

#### **2. Evolutionary Game Analysis**

Coal mining enterprises usually include frontline units and safety management departments. Workers are the main producers of behavior, which determine whether coal mining enterprises can operate safely. As an external factor, managers play a direct role in supervising the behaviors of workers, and will assume corresponding responsibilities.

Under the inspection of safety managers, workers may follow safety instructions to take actions, or may not comply with instructions to choose risky operations [31–33]. Workers are affected by safety managers, but also affect the decision-making behavior of the managers. Therefore, it can be considered that there is a game relationship between the inspection of the safety management department and the behavior choice of coal mine workers, but in the actual production of coal mine enterprises, the safety management department is facing more than a group of workers. Therefore, the paper analyses the game between a safety management department and two groups of coal mine workers.

The relevant assumptions of the game are as follows:

(1) Limited rational game group: coal mine workers Group 1, coal mine workers Group 2, coal mine safety management department. The strategies of coal mine workers groups include "safe behavior" and "unsafe behavior". Safety management departments have "inspection" and "no inspection" strategies.

(2) It is assumed that managers and coal miners are rational economic persons who determine their behavior based on cost-benefit analysis. Coal mine workers do not consciously abide by safety operating rules at all times. Similarly, managers do not always carry out safety checks, but once they do, they will strictly abide by the rules and regulations. That is to say, if workers' unsafe behaviors are found, they will be punished according to safety management regulations [34–36].

(3) Coal mine safety management departments and coal miners have limited rational characteristics. That is, there will inevitably be errors in logical reasoning or decision-making judgment and self-interest considerations. The logical reasoning or decision-making judgment is illustrated by examples. Managers would spend a lot of time checking Group 1 if they thought that Group 1 often had unsafe behavior, while Group 1 thought that the manager would not carry out safety checks immediately after checking, so the safety consciousness of Group 1 began to slacken. Consideration of self-interest is illustrated by examples. Group 2 believes that inspectors will inspect Group 1, according to the usual practice. Unsafe behaviors occur repeatedly in a short period of time, while members of groups with high safety awareness choose conformity behaviors under group pressure. From a group perspective, this may move the coal miners far away from the optimal benefit decision-making choice.

(4) The *r* is the reward for safety behaviors of coal mine workers. The *d* is the cost of taking safe actions, such as the extra physical and time required to take safety behaviors. If unsafe behavior is taken, there is no need to pay the corresponding cost. On the contrary, this part of the benefit will be obtained. The *e* is the additional benefits (such as psychological benefits, economic benefits, time benefits, etc.) for the coal miners who take unsafe behaviors. If the safety management department finds unsafe behaviors, it will take measures to correct unsafe behavior, such as education and training, and impose fines. The *F* is a penalty for unsafe behaviors (mainly economic penalties). The *H* is the loss that the group itself needs to bear, when the group has unsafe behaviors. The income obtained by the safety management department from inspecting the behaviors of coal miners mainly comes from the punishment (*F*) for the unsafe behaviors of coal miners. The *C* is the cost of safety inspection. The *D* is the loss sustained by the coal mine safety management department when the coal mine workers take unsafe behaviors. The probability of workers Group 1 adopting unsafe behaviors is *x1*, and that of Group 2 is *x*2. The probability of safety management inspection is *a1*.

The benefits matrix of Group 1 and Group 2 is shown in Table 1.


**Table 1.** The benefits matrix of Group 1 and Group 2.

*Symmetry* **2019**, *11*, 156

Group 1 chooses the benefits of safe behaviors *v*<sup>11</sup> and the benefits of choosing unsafe behaviors *v*<sup>12</sup> as Equations (1) and (2).

$$\upsilon\_{11} = (r\_1 - d\_1) \times a\_1 \times x\_2 + (1 - a\_1) \times x\_2 \times (r\_1 - d\_1) + (r\_1 - d\_1) \times a\_1 \times (1 - x\_2) + (1 - a\_1) \times (1 - x\_2) \times (r\_1 - d\_1) \tag{1}$$

$$\begin{aligned} \mathbf{v}\_{12} &= a\_1 \times \mathbf{x}\_2 \times \left( d\_1 - H\_1 - F\_1 - a\_1 \times F\_2 \right) + a\_1 \times \left( 1 - \mathbf{x}\_2 \right) \times \left( d\_1 - H\_1 - F\_1 - a\_1 \times F\_2 \right) \\ &+ \left( 1 - a\_1 \right) \times \mathbf{x}\_2 \times \left( d\_1 - H\_1 - F\_1 - a\_1 \times F\_2 \right) + \left( 1 - a\_1 \right) \times \left( 1 - \mathbf{x}\_2 \right) \times \left( d\_1 - H\_1 - F\_1 - a\_1 \times F\_2 \right) \end{aligned} \tag{2}$$

Therefore, Group 1 chooses the average expected benefits of safe behaviors and unsafe behaviors as Equation (3).

$$
\upsilon\_1 = \mathbf{x}\_1 \times \upsilon\_{11} + (1 - \mathbf{x}\_1) \times \upsilon\_{12} \tag{3}
$$

The d*x*1/d*t* represents the rate of change of the proportion of safe behaviors of Group 1 with time, and the replication dynamics of Group 1's selected safe behaviors can be obtained.

$$\mathbf{G(x\_1) = dx\_1 / dt} = \mathbf{x\_1 \times (v\_{11} - v\_1) = x\_1 \times (1 - x\_1) \times (v\_{11} - v\_{12})} \tag{4}$$

Group 2 chooses the benefits of safe behaviors *v*<sup>21</sup> and the benefits of choosing unsafe behaviors *v*<sup>22</sup> as Equations (5) and (6).

$$\mathbf{v}\_{21} = (\mathbf{r}\_2 - d\mathbf{j}) \times a\mathbf{i} \times \mathbf{x}\_1 + (\mathbf{1} - a\mathbf{j}) \times \mathbf{x}\_1 \times (\mathbf{r}\_2 - d\mathbf{j}) + (\mathbf{r}\_2 - d\mathbf{j}) \times a\mathbf{i} \times (\mathbf{1} - x\mathbf{j}) + (\mathbf{1} - a\mathbf{j}) \times (\mathbf{1} - x\mathbf{j}) \times (\mathbf{r}\_2 - d\mathbf{j}) \tag{5}$$

$$\begin{aligned} \mathbf{v\_{22}} &= a\_1 \times \mathbf{x}\_1 \times (d\_2 - H\_2 - F\_2 - a\_1 \times F\_1) + a\_1 \times (1 - \mathbf{x}\_1) \times (d\_2 - H\_2 - F\_2 - a\_1 \times F\_1) \\ &+ (1 - a\_1) \times \mathbf{x}\_1 \times (d\_2 - H\_2 - F\_2 - a\_1 \times F\_1) + (1 - a\_1) \times (1 - \mathbf{x}\_1) \times (d\_2 - H\_2 - F\_2 - a\_1 \times F\_1) \end{aligned} \tag{6}$$

Therefore, Group 2 chooses the average expected benefits of safe behaviors and unsafe behaviors as Equation (7).

$$
\upsilon\_2 = \upsilon\_2 \times \upsilon\_{21} + (1 - \upsilon\_2) \times \upsilon\_{22} \tag{7}
$$

The d*x*2/d*t* represents the rate of change of the proportion of safe behaviors of Group 2 with time, and the replication dynamics of Group 2 selected safe behaviors can be obtained.

$$\mathbf{G(x\_2) = dx\_2 / dt} = \mathbf{x\_2 \times (v\_{21} - v\_2) = x\_2 \times (1 - x\_2) \times (v\_{21} - v\_{22})} \tag{8}$$

The safety management department randomly competes with any of the two groups, and the benefits matrix is shown in Table 2.


**Table 2.** Benefits management matrix.

The benefits of safety managers' inspection and non-inspection are shown in Equations (9) and (10).

$$
\mu\_1 = \mathbf{x}\_2 \times (\mathbf{F}\_2 - D\_2) + \mathbf{x}\_1 \times (\mathbf{F}\_1 - D\_1) - \mathbf{C} \tag{9}
$$

$$
\mu\_2 = -\mathbf{x}\_2 \times D\_2 - \mathbf{x}\_1 \times D\_1 \tag{10}
$$

Therefore, the average expected benefits of safety managers' inspection and non-inspection are shown in Equation (11).

$$
u = a\_1 \times \mu\_1 + (1 - a\_1) \times \mu\_2 \tag{11}$$

The d*u*/d*t* is used to represent the change rate of the proportion of safety inspection with time, and the replication dynamics of safety inspection can be obtained.

$$\text{G(u)} = \text{d}u/\text{d}t = a\_1 \times (\text{u}\_1 - \text{u}) = a\_1 \times (1 - a\_1) \times (\text{u}\_1 - \text{u}\_2) \tag{12}$$

#### **3. Model Simulation Analysis**

#### *3.1. System Dynamics (SD) Simulation Model*

Safety management departments in coal mine enterprises face more than one workers group. Based on game theory, the paper establishes a system dynamics model of game between managers and multiple workers groups, as shown in Figure 1.

**Figure 1.** System dynamics (SD) model of evolutionary game.

The SD model of the evolutionary game has six state variables, three rate variables, 12 auxiliary variables, and 11 environmental variables. The system dynamics equation in the SD model is determined by the replication dynamic equation in the evolutionary game model.

#### *3.2. Strategies Simulation Analysis*

According to the above variables settings, the system dynamics model is analyzed by Vensim software (PLE, Ventana Systems. Inc, Harvard, MA., USA). The simulation results are shown in Figures 2–4.

**Figure 2.** The probability of manager safety inspection under general state and dynamic penalty.

**Figure 3.** The probability of Group 1 choosing unsafe behaviors under general state and dynamic penalty.

**Figure 4.** The probability of Group 2 choosing unsafe behaviors under general state and dynamic penalty.

In Figures 2–4, curves 1, 2, 4, and 5 show the dynamic evolution of strategies of each player when managers choose dynamic penalties. Curves 3 and 6 are the dynamic evolution of each player's strategy, when the probability of safety inspection is 50% and the unsafe actor is punished generally.

Curve 3 is the dynamic evolution of each player's strategy, when the level of unsafe behavior of both groups is 50%, the level of unsafe behavior of another group is 70%. Curve 6 is the dynamic evolution of the strategies of each player when the level of unsafe behaviors of two groups is 50%. Contrast curves 3 and 6 show that the level of unsafe behaviors of coal miners has an effect on the safety inspection of managers, while the level of unsafe behaviors of group 2 has little effect on the level of unsafe behaviors of group 1. Curves 3 and 6 in Figure 2 show that with the continuation of the game, asymmetry of the benefits between safety managers and workers began to appear. The level of group unsafe behaviors decreases, and the source of benefits for safety managers is 0. The safety manager chooses the safety probability not to rise but to fall. This will cause a great potential safety hazard to coal mining enterprises, which is not in line with the actual situation.

Curves 1 and 4, 2 and 5, and 3 and 6 show that when the punishment is positively correlated with the level of unsafe behavior, the greater the correlation coefficient is, the faster the level of unsafe behavior decreases, and the shorter the time it takes for the final group to choose safe behaviors. When the level of unsafe behaviors of two groups of coal miners is different, the two sides also influence each other, and the effect is obvious. This is because the level of safety behaviors of both sides has an impact on the behaviors of safety managers, and the punishment of the groups is also related to the level of unsafe behaviors of the other group.

Part of the benefit from coal mine management comes from fines for unsafe behaviors. Therefore, in order to maximize profits, safety management departments tend to choose safety inspection as their initial strategy. The higher the level of group unsafe behavior is, the more likely managers will choose safety management. With the improvement of the level of group safety behaviors, the benefit of the coal mine safety management department decreases. When the benefit of inspection is less than that of non-inspection, the probability of inspection decreases and finally evolves into a hidden danger state.

Therefore, only considering dynamic punishment can reduce the level of group unsafe behavior to a certain extent, but it can't motivate safety managers to carry out safety inspection continuously. This should consider how to adjust the benefit symmetry of managers and workers. In order to stimulate safety managers to conduct continuous safety inspections, dynamic incentives are considered here. The magnitude of dynamic incentive coefficient is negatively correlated with the level of group unsafe behaviors, that is, the greater the level of group safe behaviors, the greater the incentive of coal mine enterprises to the coal mine safety management department, and the greater the incentive to the group. The SD model with dynamic incentives is shown in Figure 5, and the simulation results are shown in Figures 6–8.

**Figure 5.** SD model with dynamic incentives.

**Figure 6.** The probability of managers choosing safety inspection with dynamic incentives.

**Figure 7.** The probability of Group 1 choosing unsafe behaviors with dynamic incentives.

**Figure 8.** The probability of Group 2 choosing unsafe behaviors with dynamic incentives.

From the curves 2 and 4 of Figure 6, it can be seen that introducing dynamic incentives can improve the symmetry of game benefits between managers and workers. At the same time, it can also motivate safety managers to carry out safety inspections continuously and ensure the stability of the safe operation of coal mine workers. It can also be seen that when the level of unsafe behavior of two groups is not the same, the lower the level of unsafe behavior of the groups, the more quickly safety managers can make a safety inspection.

It can be seen from Figures 7 and 8 that after with the dynamic incentives, part of the benefits of coal mine workers depend on the level of safety behaviors of the groups. Therefore, the order of the decline rate of unsafe behavior level is: Dynamic incentive > Dynamic punishment > General punishment.

It is concluded from the simulation analysis that the dynamic incentives have a good effect on improving the benefit symmetry and behavior stability of managers and workers, and can promote the continuous supervision of safety managers and reduce the level of unsafe behavior of coal miners. Therefore, in the process of safety production, coal mine enterprises should give reasonable rewards and punishments to safety managers and workers to regulate the benefit symmetry of managers and workers and encourage them to choose safety behaviors, so as to ensure the smooth progress of safety production.

#### **4. Example Application**

Through the game analysis and SD simulation analysis of the safety management department and workers groups, it can be seen that dynamic incentives play a positive role in the stability of safety behaviors. Under the guidance of the theory, the Tangkou Coal Industry has been investigated on the spot. The existing violation treatment system has been optimized, the main content of which is to introduce a dynamic incentive system. Figure 9 shows the trend of violations of regulations in the Machinery and Electricity Team of Tangkou Coal Industry from January 2017 to November 2018. The team began to implement the dynamic incentive system in November 2017.

**Figure 9.** The trend of violations of regulations in the Machinery and Electricity Team from January 2017 to November 2018.

As can be seen from Figure 9, before the implementation of the dynamic incentive system (January 2017–October 2017), the number of violations of the team has been between 170–220. However, starting from November 2017, after the implementation of the dynamic incentive system, the number of violations of the team dropped rapidly. The number of violations dropped to below 70 in August 2018. The number of violations was stable at around 60 in September 2018, October 2018, and November 2018.

This paper compares and analyzes the violations of the team in the same month in 2017 and 2018, in order to more clearly analyze the changes in violations before and after the implementation of the dynamic incentive system, as shown in Figure 10. It can be clearly seen from Figure 10 that after the dynamic incentive measures were adopted, the number of violations of the team decreased significantly during the same period.

**Figure 10.** Comparison of the violations of the team in 2017 and 2018.

This shows that after the dynamic incentive measures, the safety management departments are motivated to be more active in the regulation of violations. At the same time, under the dynamic incentive measures, workers can also get extra benefits from safety production in the process of work, and their safety awareness is also improved.

It can be seen from Figures 9 and 10 that although the number of violations of the team has been significantly reduced, it is still around 60. This shows that there is still room for decline in the team's violations. Combined with the actual situation of Tangkou Coal Industry, if the dynamic incentive system can continue to be properly adjusted, the number of violations of the team may continue to decline. This is also the focus of the next research work in the paper.

#### **5. Conclusions**

The paper studies the symmetry of game benefits and behavior stability between coal mine safety managers and two groups of coal mine workers. Through the SD model, the evolutionary game process between safety management department and two groups of coal miners is analyzed in detail. The main conclusions are as follows:

(1) The symmetry of game benefits between safety managers and workers is analyzed, and their benefit matrix and game model is established.

(2) On the basis of the benefit matrix and game model, the SD model is constructed. The evolutionary process of game between the safety management department and two workers groups is analyzed in detail. The research shows that general punishment and dynamic punishment have a certain effect on reducing the safety behavior of coal miners, but at the same time, they also have some hidden dangers. Because the main source of benefit of coal mine safety managers is the fines for unsafe behavior, when the benefit symmetry is broken and the benefits of the safety management department's safety inspection is less than that of non-inspection, the safety management department ultimately does not carry out the safety inspection, which can be a great hidden danger to coal mine safety production.

(3) The dynamic incentive system is added to the SD model. It is concluded from the simulation analysis of the SD model that dynamic incentives can effectively regulate benefit symmetry and behavior stability of managers and workers, and promote continuous inspection by safety managers and reduce the level of unsafe behavior of coal mine workers.

(4) The research results of this paper have been applied to coal mine enterprises. By introducing dynamic incentive system, the existing reward and punishment system of enterprises has been adjusted and optimized. It effectively reduces the number of violations in coal mine enterprises.

The simulation analysis and practical application in this paper show that the establishment of a reasonable dynamic incentive mechanism can effectively improve the safety level of coal mine enterprises in the process of safety production.

**Author Contributions:** Conceptualization, K.Y. and Z.L.; methodology, K.Y. and L.Z.; validation, K.Y., L.Z., and Q.C.; investigation, K.Y., L.Z., and Q.C.; data curation, L.Z.; writing—original draft preparation, K.Y.; writing—review and editing, L.Z. and Q.C.; visualization, K.Y. and Z.L.; supervision, L.Z. and Q.C.; project administration, L.Z. and Q.C.; funding acquisition, L.Z. and Q.C.

**Funding:** This research was funded by National Natural Science Foundation of China, grant number 51474138; 51574157; 51804180.

**Acknowledgments:** The authors would like to thank the National Natural Science Foundation of China and the authors of the references.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Route Choice Behavior: Understanding the Impact of Asymmetric Preference on Travelers' Decision Making**

#### **Kai Liu \* and Yuan Xu \***

School of Engineering and Management, Nanjing University, Nanjing 210093, China

**\*** Correspondence: dg1415007@smail.nju.edu.cn (K.L.); xy\_0011@163.com (Y.X.); Tel.: +86-159-520-35866 (K.L.); +86-180-129-68631 (Y.X.)

Received: 10 December 2018; Accepted: 4 January 2019; Published: 8 January 2019

**Abstract:** This paper investigated the impact of asymmetric preference on travelers' route choices. Firstly, a status quo-dependent route choice mode was developed to describe travelers' route choices. Then, based on that model, a route choice experiment was conducted, and during the experiment, participants were requested to choose a route from two arbitrary non-dominated routes. Finally, according to the observation data, data analysis and model parameter estimation were conducted. The results show that participants used different measures to trade off travel cost and travel time. Additionally, there was a gap between most participants' willingness to pay (WTP) and willingness to accept (WTA). Moreover, participants' WTP greater than their own WTA was the key reason resulting in the inertial route choices. The empirical results in this paper can help the traffic manager to understand travelers' inertial route choice behavior from a different perspective.

**Keywords:** asymmetric preference; inertial route choice; willingness to pay; willingness to accept; experimental study

#### **1. Introduction**

Travelers' route choice behaviors have been researched for many years. Traditionally, it was assumed that all travelers would choose the routes with the shortest travel time [1]. However, in reality, this assumption is very restrictive, because it ignores travelers' cognitive limitations and intrinsic preferences [2–4]. Moreover, many empirical studies have shown that travelers do not always choose the shortest routes [5–7]. Simon [8] indicates that people are boundedly rational, which means sometimes they prefer a satisfactory choice to an optimal one. Based on the concept of bounded rationality, Mahmassani and Chang [9] introduced the concept of the indifference band to study travelers' boundedly rational departure time choices. Lou et al. [10] further encapsulated the indifference band into route choice modeling and proposed a boundedly rational route choice model, which addressed the limitation on travelers' knowledge of traffic conditions and their capability of finding the best available routes. In their boundedly rational route choice model, travelers would not necessarily switch to the shortest route when the time (or cost) difference between the current route and the shortest one was lower than an inertia threshold. Zhao and Huang [11] applied the concept of aspiration level to investigate travelers' boundedly rational route choice behavior and found that travelers can only obtain the actual travel times of the chosen routes and are enabled to recognize the best routes.

The boundedly rational behavior in travel choice describes the facts that travelers will not always choose the routes with the shortest travel time. There are various causes raising such behavior, one of which is inertial choice. In fact, the term "inertia" was first mentioned in Newtonian physics, and it has been widely used to describe similar characteristics of human behavior. Inertial choice can be understood as the tendency to repeat previous choices [12]. In the domain of behavioral decision, inertia was always employed to characterize the tendency for decision makers to choose options that maintain the status quo regardless of its outcome in a subsequent decision scenario, though the concise definitions were slightly different among the literature [13–19].

In the area of travel behavior research, travelers' inertia was defined as the tendency to stick to the alternative that they had previously chosen [15,20–24]. Many empirical studies provided support for the inertial route choice behavior [7,25,26]. For example, Chorus and Dellaert [15] thought that the wish to save cognitive resources could lead travelers to inertial choices; Site and Filippi [27] explained that phenomenon in terms of "loss aversion"—the disadvantages of a move from the status quo are valued more heavily than the advantages. Besides saving cognition and avoiding loss aversion, other behavioral mechanisms which generate inertial behavior include the asymmetric preference [4,22], prevailing choice set [28], habitual behavior [29,30], familiarity, prior decision [21], risk aversion [31], endowment effect [32], and so on.

The asymmetric preference describes the phenomenon that travelers might feel different when they made a tradeoff between travel cost and travel time depending on whether they spend money to save time (willingness to pay, WTP) or they bear more time to receive money (willingness to accept, WTA) [33]. While in the conventional study of the value of time, the substitutions between money and time are assumed to be constant. Various theoretical and empirical studies have been developed to examine the value of WTP and WTA as well as to estimate the gap between those two values [34–38]. De Borger and Fosgerau [35] adopted the status quo as the reference state to distinguish travelers' judgment between gain and loss, and further examined the trade-off between money and time. Their empirical study confirmed the gap between WTP and WTA. However, De Borger and Fosgerau [35] used "wrong choice" to describe travelers' route choices in their research. In fact, travelers' route choices were not "wrong", their behavior was just coinciding with the inertial choice.

Therefore, based on the research of De Borger and Fosgerau [35], Xu et al. [4] proposed a status quo-dependent route choice model to handle the so called "wrong choice" behaviors from the view of inertial choice. In their proposed route choice model, travelers were assumed to "compare the travel cost to their status quo (travel cost of the currently used path) in deciding whether to switch to another alternative, and the underlying value of time is adaptive in the sense that it varies across different route choice contexts". Xu et al. [4] found that the status quo-dependent route choice model can explain the route choice inertia resulting from asymmetric preference. Moreover, the inertia is path-specific and can incorporate the scaling effect of travel cost on travelers' route choices.

However, there was a lack of empirical study in Xu et al.'s [4] research, thus they cannot clearly show that when travelers' WTP is smaller than their WTA, or when travelers' WTP is greater than their WTA, how travelers will make their route choices. In addition, how will the scaling effect of travel cost affect travelers' asymmetric preference? This paper is an extended work of Xu et al.'s research aiming to provide more detailed support as an empirical study. A route choice experiment was conducted to investigate the impact of asymmetric preference on travelers' inertial route choice. In the experiment, three different travel scenarios (short, middle, and long travel) were designed to collect participants' route choice data. Then, based on the experimental data, the value of participants' WTP and WTA were calculated and compared. According to the value of WTP and WTA, participants' route choice behaviors were analyzed. Moreover, the parameters in Xu's route choice model were also estimated.

The remainder of this paper is organized as follows. Section 2 describes the status quo-dependent route choice model on how to handle travelers' route choice behaviors that result from asymmetric preference. Section 3 describes the route choice experiment as well as the data collection procedure. Section 4 presents the calculation and analysis results of participants' WTP and WTA and route choices in the experiment. Section 5 presents the model parameters estimation. Finally, Section 6 concludes the paper.

#### **2. The Status Quo-Dependent Route Choice Model**

For the convenience of readers, the symbols used in this section are summarized as follows:


Considering two arbitrary non-dominated routes *r* and *j*, *r*, *j* ∈ *Rw*, *w* ∈ *W*, without loss of generality, it is assumed that *t w <sup>r</sup>* ≤ *t w <sup>j</sup>* and *<sup>τ</sup><sup>w</sup> <sup>r</sup>* ≥ *<sup>τ</sup><sup>w</sup> <sup>j</sup>* . For the travelers on route *r*, the status quo-dependent travel utility by changing to path *j* is calculated as Equation (1).

$$
\Delta I \left( c\_r^w, c\_j^w \right) = -\lambda\_1 \eta\_1 \left( t\_j^w - t\_r^w \right) + \eta\_2 \left( \tau\_r^w - \tau\_j^w \right) \tag{1}
$$

where *U cw <sup>r</sup>* , *c<sup>w</sup> j* is the status quo-dependent travel utility for the travelers on route *r*, with *c<sup>w</sup> r* representing the status quo and *c<sup>w</sup> <sup>j</sup>* the alternative. Parameter *λ*<sup>1</sup> is the travelers' loss aversion coefficient of the time cost. *η*<sup>1</sup> and *η*<sup>2</sup> are the coefficients that unify time and monetary factors with travel utility. Thus, travelers on route *r* would suffer a smaller travel time, but they would obtain a larger cost in the aspect of monetary factors, and travelers would stick to route *r* if and only if the inequality *U* - *cw <sup>r</sup>* , *c<sup>w</sup> j* ≤ 0 holds. Otherwise, the travelers on route *r* will change to route *j*. In addition, *U* - *cw <sup>r</sup>* , *c<sup>w</sup> j* ≤ 0 holds if and only if

$$-\lambda\_1 \eta\_1 \left( t\_j^w - t\_r^w \right) + \eta\_2 \left( \tau\_r^w - \tau\_j^w \right) \le 0 \Rightarrow \tau\_r^w - \tau\_j^w \le \frac{\lambda\_1 \eta\_1}{\eta\_2} \left( t\_j^w - t\_r^w \right). \tag{2}$$

As indicated by Equations (1) and (2), travelers on route *r* would not change to route *j* if the potential improvement of monetary cost is not great enough (i.e., less than or equal to the right-hand-side value), and the coefficient *λ*1*η*1/*η*<sup>2</sup> in Equation (2) is equal to the travelers' willingness to accept (WTA).

On the other hand, for the travelers on route *j* the status quo-dependent travel utility by changing to route *r*, is calculated as Equation (3).

$$
\Delta I\left(c\_{\dot{\gamma}}^{w}, c\_{r}^{w}\right) = \eta\_{1}\left(t\_{\dot{\gamma}}^{w} - t\_{r}^{w}\right) - \lambda\_{2}\eta\_{2}\left(\pi\_{r}^{w} - \pi\_{\dot{\gamma}}^{w}\right) \tag{3}
$$

where *U cw <sup>j</sup>* , *<sup>c</sup><sup>w</sup> r* is the status quo-dependent travel utility for the travelers on route *j*, with *c<sup>w</sup> j* representing the status quo and *c<sup>w</sup> <sup>r</sup>* the alternative. Parameter *λ*<sup>2</sup> is the travelers' loss aversion coefficient of the monetary cost. Travelers on route *j* would suffer a smaller monetary cost, but they would get a larger cost in the aspect of travel time and travelers would stick to route *j* if and only if the inequality *U* - *cw <sup>j</sup>* , *<sup>c</sup><sup>w</sup> r* ≤ 0 holds. Otherwise, travelers on route *j* will change to route *r*. In addition, *U* - *cw <sup>j</sup>* , *<sup>c</sup><sup>w</sup> r* ≤ 0 holds if and only if

$$
\eta\_1 \left( t\_{\dot{j}}^w - t\_r^w \right) - \lambda\_2 \eta\_2 \left( \tau\_r^w - \tau\_{\dot{j}}^w \right) \le 0 \\
\Rightarrow \tau\_r^w - \tau\_{\dot{j}}^w \ge \frac{\eta\_1}{\lambda\_2 \eta\_2} \left( t\_{\dot{j}}^w - t\_r^w \right). \tag{4}
$$

As indicated by Equations (3) and (4), travelers on route *j* would not change to route *r* if the increase in monetary cost is greater than or equal to the right-hand-side value, and the coefficient

*η*1/*λ*2*η*<sup>2</sup> in Equation (4) is equal to the travelers' willingness to pay (WTP). Then, based on the above analysis, a route choice experiment will be conducted in the following section to verify travelers' inertial route choices that resulted from asymmetric preference.

#### **3. Experiment Design**

The experiment was carried out in the laboratory of the School of Engineering and Management, Nanjing University, in November 2016. We used the z-Tree as our experimental platform [39] to deal with the inputs of all participants' route choices. Next, we will introduce in detail the participants, experimental design, and procedure in the experiment.

#### *3.1. Participants*

In total, 30 participants—12 females and 18 males—were recruited from Nanjing University, and all of them had never participated in a similar experiment. Participants were required to make route choice decisions for a total of 120 test rounds, according to the experimental settings and procedures. During the experiment, participants were not allowed to have any interaction with others.

#### *3.2. Experimental Design*

The experimental design included three travel scenarios reflecting three different levels of travel time—short, middle, and long. In each of the scenarios, there were two non-dominated routes, e.g., Route 1 and Route 2, and the route's status quo is characterized by (*tr*, *τr*). During the experiment, participants needed to choose one's routes as their travel route, and in each scenario, participants would make route choices for 40 test rounds. Since we focused on travelers' route choice behaviors in this study, thus, the interaction between travelers and the collective effect of travelers' behavior on congestion are not considered in the experiment design. Tables 1 and 2 show the design of travel time and travel cost of two routes in the experiment.

**Table 1.** Investigating participants' willingness to pay (travel time \*, monetary cost \*\*).


\* The unit of travel time is minute. \*\* The unit of monetary cost is Yuan, which is the Chinese currency unit.


**Table 2.** Investigating participants' willingness to accept (travel time, monetary cost).

Table 1 is used to investigate participants' WTP, and Table 2 is used to investigate participants' WTA. *τ<sup>i</sup>* <sup>1</sup>, *<sup>τ</sup><sup>i</sup>* <sup>2</sup>, and *<sup>τ</sup>*ˆ*<sup>i</sup>* <sup>2</sup>(*i* = 1, 2, 3) represent travel money costs in the three different travel scenarios, and to note that the values of *τ<sup>i</sup>* <sup>1</sup>, *<sup>τ</sup><sup>i</sup>* <sup>2</sup>, and *<sup>τ</sup>*ˆ*<sup>i</sup>* <sup>2</sup> are calculated based on the participants inputting answers about their WTP at the initial stage of the experiment. Figure 1 shows the experimental design steps. As shown in Figure 1, for a participant, at the initial stage of the experiment they needed to input their answers for "if the travel time can save 10 min, how much you willing to pay", "if the travel time can save 25 min, how much you willing to pay", and "if the travel time can save 45 min, how much you willing to pay". Then, based on the three inputted answers, we can calculate the values of *τ<sup>i</sup>* <sup>1</sup>, *<sup>τ</sup><sup>i</sup>* 2, and *τ*ˆ*<sup>i</sup>* <sup>2</sup>(*i* = 1, 2, 3) in Tables 1 and 2. Additionally, from participants' own points of view, we can get the (25, 3) is equal to 10, *τ*<sup>1</sup> 2 , 24, *τ*<sup>1</sup> 1 , and 12, *τ*ˆ<sup>1</sup> 2 (similarly, we can also get the (50, 10) is equal to 30, *τ*<sup>2</sup> 2 , 45, *τ*<sup>2</sup> 1 , and 25, *τ*ˆ<sup>2</sup> 2 ; (90, 20) is equal to 60, *τ*<sup>3</sup> 2 , *c*<sup>1</sup> = 85, *τ*<sup>3</sup> 1 , and 65, *τ*ˆ<sup>3</sup> 2 ). Note that in Tables 1 and 2 the parameters *α<sup>i</sup> <sup>n</sup>* and *β<sup>i</sup> <sup>k</sup>*(*i* = 1, 2, 3) are constants, which are appended to the monetary

cost, and the *n* and *k* index for the round of experiment. Moreover, as the experiment progressed, there were *α<sup>i</sup>* <sup>1</sup> > *<sup>α</sup><sup>i</sup>* <sup>2</sup> > ··· > *<sup>α</sup><sup>i</sup> <sup>t</sup>* = <sup>0</sup> > ··· > *<sup>α</sup><sup>i</sup> <sup>T</sup>* and *<sup>β</sup><sup>i</sup>* <sup>1</sup> > *<sup>β</sup><sup>i</sup>* <sup>2</sup> > ··· > *<sup>β</sup><sup>i</sup> <sup>t</sup>* = <sup>0</sup> > ··· > *<sup>β</sup><sup>i</sup> <sup>T</sup>*. This design ensures that participants will choose Route 1 in the inertial test rounds when investigating their WTP, and they will choose Route 2 in the inertial test rounds when investigating their WTA.

**Figure 1.** The design of experiment. WTP: willingness to pay; WTA: willingness to accept.

#### *3.3. Procedure*

Participants were invited to participate in the experiment at the University laboratory, and each participant was seated separately in front of a computer screen displaying the simulation. In order to encourage realistic behavior, participants were asked to imagine that they were commuters and driving to the workplace. The monthly salary was three times of their school living expenses. They were also informed that from the home to the workplace there were two routes that could be chosen. One route's travel time was larger than another route's travel time, but the travel money cost was smaller. For example, Route 1 was a congested road but toll-free, thus the monetary travel cost on Route 1 only consisted of an emission fee, and vehicle operating cost; Route 2 was a toll road but not congested, thus the monetary travel cost on Route 2 consisted of a congestion toll, emission fee, and vehicle operating cost. Note that in the cases of WTA, although a constant *β<sup>i</sup> <sup>k</sup>* was appended to the Route 1 s monetary cost, Route 1 was still cheaper than Route 2 in the experiment. Figure 2 gives examples of actual choice situations in the experiment. No other information was provided before beginning the experiment.

**Figure 2.** Snapshot of the route choice simulation window (Scenario 1).

The task included making a series of 120 consecutive route choices. Each test round represented a working day and participants were asked to make daily route choices while driving from home to workplace. Moreover, before the experiment, participants were requested to complete a post-task questionnaire that enquired about their loss aversion of the time cost and the monetary cost. The experimental steps are specified as follows:

(1) At the initial time, each participant needed to input the answers about their WTP, i.e., how much you are willing to pay if travel time can save 10 min; how much you are willing to pay if travel time can save 25 min; and how much you are willing to pay if travel time can save 45 min.


#### **4. Results and Discussion**

#### *4.1. The Values of Participants' WTP and WTA*

Based on the experimental data, the values of participants' WTP and WTA can be calculated. Note that the WTP refers to the value of a unit time saving in the case of loss in monetary cost, and WTA refers to the monetary compensation a traveler requires to agree to a unit travel time increase. Thus, in the experiment, the participants' WTP and WTA can be calculated as shown in Figure 3. It is worth mentioning that, for some participants, their initial stated WTP was different from the experimental WTP, in other words, this means that the result of the report is not equal to the result of the action. Considering the action result can better reflect participants' psychological index of decision-making, therefore, the experimental WTP and WTA were used in this paper. Figure 4 shows the calculation results. It clearly shows that for most participants, their WTP was smaller than WTA in all of three travel scenarios. This means that for these participants, the monetary compensation for the value of unit travel time increasing is larger than the monetary payment for the same value of unit travel time saving. Moreover, the test of null hypothesis shown in Table 3 suggests that the value of WTP is not equal to the value of WTA.

**Figure 3.** The calculation method of WTP and WTA.

**Figure 4.** *Cont.*

**Figure 4.** The value of WTP and WTA for all participants in three travel scenarios.



However, regardless of WTP < WTA or WTP > WTA, this result confirms that participants use different measures to distinguish the value of time, rather than constant substitution between money and travel time [33]. Thus, in the experiment, the asymmetric preference exists in participants' route choices. In addition, for the average values of the WTP, there is WTP1 < WTP2 < WTP3 (*i* = 1, 2, 3 represents Scenario 1, Scenario 2, and Scenario 3). This result indicates that with the increase of the travel time, participants were willing to pay more to save travel time. Then, in order to explain the impact of asymmetric preference on participants' inertial route choices, in the following subsection, participants' route choice behaviors will be analyzed.

#### *4.2. Influence of Asymmetric Preference on Participants' Route Choices*

Since in the experiment, many participants' WTP < WTA or WTP > WTA, the impact of asymmetric preference on participants' route choices will be discussed with the two cases, i.e., WTP < WTA and WTP > WTA.

#### *Case 1:* WTP < WTA

In this subsection, Participant 2 will be taken as an example, and Figures 5–7 show Participant 2 s route choices in the experiment under three different travel scenarios. As shown in Figure 5, in Scenario 1, Participant 2 s WTP equals 0.2 and WTA equals 0.33. Then, assuming that if Participant 2 s WTP equals WTA (i.e., WTP = WTA = 0.33), they should switch to Route 2 when its travel monetary cost is 8. However, in fact, Participant 2 continues to choose Route 1 until Route 2 s travel monetary cost is 6. On the other hand, assuming that if Participant 2 s WTA equals WTP (i.e., WTP = WTA = 0.2), they should switch to Route 1 when its travel monetary cost is 5.5. However, in fact, Participant 2 continues to choose Route 2 until Route 1 s travel monetary cost is 4. Moreover, for Participant 2 (even those other participants whose WTP was smaller than WTA), similar phenomena can be found in the other two scenarios. Therefore, when the participants' WTP was smaller than WTA, the asymmetric preference can make them stick to the current routes until one alternative route's monetary cost is enough to tempt them to choose. This analysis result verifies that asymmetric preference (i.e., WTP < WTA) can lead travelers to make inertial choices (e.g., stick to the current route).

**Figure 6.** Scenario 2—middle travel time.

**Figure 7.** Scenario 3—long travel time.

*Case 2:* WTP > WTA

In this subsection, Participant 10 (Scenario 1), Participant 22 (Scenario 2), and Participant 21 (Scenario 3) will be taken as examples, and Figures 8–10 show participants' route choices in the experiment under three different travel scenarios. As shown in Figure 8, in Scenario 1, Participant 10's WTP equals 0.47 and WTA equals 0.33. Then, assuming that if Participant 10's WTP equals WTA (i.e., WTP = WTA = 0.33), they should switch to Route 2 when its travel monetary cost is 8. However, in fact, Participant 10 switches to Route 2 when its travel monetary cost is 11. On the other hand, assuming that if Participant 10's WTA equals WTP (i.e., WTA = WTP = 0.47), they should switch to Route 1 when its travel monetary cost is 2.5. However, in fact Participant 10 switches to Route 1 when its travel monetary cost is 4.2. Moreover, for another two participants, similar phenomena can be found in the other two scenarios. Therefore, when participants' WTP was larger than WTA, the asymmetric preference can make them tend to choose a route with the shortest travel time.

**Figure 8.** Scenario 1—short travel time (Participant 10).

**Figure 9.** Scenario 2—middle travel time (Participant 22).

**Figure 10.** Scenario 3—long travel time (Participant 21).

The above analysis results suggest that when participants' WTP is smaller than WTA, the asymmetric preference can make them stick to the current routes. Next, based on the observation data, the parameters *λ*1, *λ*2, *η*1, and *η*<sup>2</sup> in Xu's [4] status quo-dependent route choice model will be estimated.

#### **5. Model Parameters Estimation**

#### *5.1. Estimation of λ*<sup>1</sup> *and λ*<sup>2</sup>

The cumulative prospect theory (CPT), which was proposed by Tversky and Kahneman [40], is a valid utility measurement system. Many researchers used CPT to study individual's risk attitude and choice behavior under uncertainty over the past years. Xu et al. [41] developed a general travel decision-making rule utilizing CPT to investigate travelers' risk attitudes and route choice behaviors under uncertainty. Zhang et al. [42], based on the CPT, developed a day-to-day route-choice learning model with friends' travel information. Therefore, in this subsection, the cumulative prospect theory (CPT) was applied to measure participants' risk attitudes, i.e., loss aversion of travel time and travel monetary cost. Without loss of generality, let *x*<sup>0</sup> define a reference point in the outcome domain, and the value function can be defined as Equation (5), and the weighting function can be defined as Equation (6).

$$w(\mathbf{x}) = \begin{cases} (\mathbf{x} - \mathbf{x}\_0)^a, \mathbf{x} > \mathbf{x}\_0 \\\ -\lambda \cdot (\mathbf{x}\_0 - \mathbf{x})^\mathcal{G}, \mathbf{x} \le \mathbf{x}\_0 \end{cases} \tag{5}$$

$$w(p\_i) = p\_i^{\gamma} / \left[p\_i^{\gamma} + (1 - p\_i)^{\gamma}\right]^{1/\gamma} \tag{6}$$

where the parameters *α* ≤ 1, *β* ≤ 1 measure the degree of the diminishing sensitivity to change in both directions from the reference point *x*0, parameter *λ* ≥ 1 captures the degree of loss aversion, parameters 0 < *γ* < 1 reflect the level of distortion in the probability judgment. *pi* is the probability of an outcome *xi*. Figure 11 illustrates the value function and the weighting function. Suppose an alternative denoted by the pair (*x*; *p*) is composed of *m* + *n* + 1 possible outcomes *x*−*<sup>m</sup>* < ··· < *x*<sup>0</sup> < ··· < *xn* with probabilities *p*−*m*, ··· , *pn*, respectively. Then, the cumulative decision weights are defined as follows:

$$w\_i^+(p\_i) = w^+(p\_i + \dots + p\_n) - w^+(p\_{i+1} + \dots + p\_n), 0 \le i \le n - 1\tag{7}$$

$$
\pi\_{-j}^{-}(p\_{-j}) = w^{-} \left( p\_{-m} + \dots + p\_{-j} \right) - w^{-} \left( p\_{-m} + \dots + p\_{-j-1} \right), 1 - \mathbf{m} \le -j \le 0. \tag{8}
$$

Accordingly, the CPT value of (*x*; *p*) can be calculated by Equation (9).

$$MI(\mathbf{x}; p) = \sum\_{i=0}^{n-1} v(\mathbf{x}) \cdot \pi\_i^+(p\_i) + \sum\_{j=1-m}^0 v(\mathbf{x}) \cdot \pi\_{-j}^-(p\_{-j}) \tag{9}$$

**Figure 11.** Value function and weighting function.

Firstly, participants were requested to complete a post-task questionnaire. Then, based on the survey data from the questionnaire, the parameter *λ*<sup>1</sup> and *λ*<sup>2</sup> were estimated. The design of the questionnaire considered the following two scenes:

**Scene 1** (Investigate the loss aversion of the travel time): The travel times on both routes (Route A and Route B) are likely to increase because of congestion. Let (*t*1, *p*1%; *t*2, *p*2%) denote the distribution of travel time on each route. For a traveler, there is a probability of *p*1% to accomplish the trip with *t*<sup>1</sup> more minutes, a probability of *p*2% to accomplish the trip with *t*<sup>2</sup> more minutes, and a probability of 1 − *p*1% − *p*2% with no extra time.

**Scene 2** (Investigate the loss aversion of the travel monetary cost): The travel costs on both routes (Route R and Route S) are likely to increase because of the late penalty. Let (*c*1, *p*1%; *c*2, *p*2%) denote the distribution of travel cost on each route. For a traveler, there is a probability of *p*1% to accomplish the trip with *c*<sup>1</sup> more costs (Yuan, Chinese currency), a probability of *p*2% to accomplish the trip with *c*<sup>2</sup> more costs, and a probability of 1 − *p*1% − *p*2% with no extra costs.

The route choice preferences of participants between A and B are shown in Table 4. The route choice preferences of participants between R and S are shown in Table 5.


**Table 4.** Results of route choice in Scene 1.



Then, based on the survey data in Tables 4 and 5, the values of parameters *λ*<sup>1</sup> and *λ*<sup>2</sup> can be estimated using the least squares method and the cumulative squared residual can be given by following Equations (10) and (11), respectively:

$$f(\beta\_1, \lambda\_1, \delta) = \sum\_{k=1}^{k=5} (A\_k \% - P(A\_k > B\_k)) \tag{10}$$

$$f(\beta\_{2\prime}\lambda\_{2\prime}\delta) = \sum\_{k=6}^{k=12} \left( R\_k \% - P(R\_k > S\_k) \right)^2 \tag{11}$$

where *A*% is the proportion of choosing Route *A* in case *k*, *P*(*Ak* > *Bk*) is the probability that Route *A* is more valuable than Route *B* in case of *k*, and *P*(*Ak* > *Bk*) can be given by Equation (12). *R*% is the proportion of choosing Route *R* in case of *k*, *P*(*Rk* > *Sk*) is the probability that Route *R* is more valuable than Route *S* in case *k*, and *P*(*Rk* > *Sk*) can be given by Equation (12). Here, *U*(·) is the route's cumulative prospect value and it can be calculated by Equation (9).

$$P(A\_k > B\_k) = \frac{1}{1 + \exp(\mathcal{U}(B\_k) - \mathcal{U}(A\_k))}\tag{12}$$

$$P(R\_k > S\_k) = \frac{1}{1 + \exp\left(\mathcal{U}(S\_k) - \mathcal{U}(R\_k)\right)}\tag{13}$$

Note that in this paper, the value of parameter *δ* is set as 0.74, which was estimated by Wu and Gonzalez [43]. Participants' perception errors of the cumulative prospect value are assumed to follow the independent and identical Gumbel distributions. Then, based on these settings, the estimation values of *β*<sup>1</sup> is 0.35 and *λ*<sup>1</sup> is 1.41 in Scene 1; *β*<sup>2</sup> is 0.23 and *λ*<sup>2</sup> is 1.26 in Scene 2.

#### *5.2. Estimation of η*<sup>1</sup> *and η*<sup>2</sup>

Until now, we have estimated the values of parameters *λ*<sup>1</sup> and *λ*2, associated with the CPT-based utility measurement system. Specifically, *λ*<sup>1</sup> = 1.41 and *λ*<sup>2</sup> = 1.26. With them, based on the experimental data, we used the least squares method to estimate the values of parameters *η*<sup>1</sup> and *η*2. Data of three thousand and six hundred observations of 30 participants were used. Table 6 presents the estimation results.


**Table 6.** Model estimation results.

As shown in Table 6, the estimation values of all parameters are positive and significant. Based on the estimation values of *λ*1, *λ*2, *η*1, and *η*2, participants' WTP (i.e., *η*1/*λ*2*η*2) and WTA (i.e., *λ*1*η*1/*η*2) can be calculated. Table 5 also shows the results of the calculation. As shown in Table 6, there is *WTP*<sup>1</sup> < *WTP*<sup>2</sup> < *WTP*<sup>3</sup> and *WTA*<sup>1</sup> < *WTA*<sup>2</sup> < *WTA*<sup>3</sup> (*i* = 1, 2, 3 represents Scenario 1, Scenario 2, and Scenario 3). These results are consistent with the aforementioned data analysis. Moreover, the calculation results of *WTP* < *WTA* for three scenarios further explain that the asymmetric preference exists in participants' route choices.

#### **6. Conclusions**

In this research, a route choice experiment was conducted to investigate travelers' inertial route choice behaviors that resulted from asymmetric preference. Through data analysis, we can confirm that in the experiment participants used different measures to distinguish the value of time rather than constant substitution between money and travel time. Additionally, asymmetric preference significantly exists in participants' route choices. Especially, for most participants, their WTP was smaller than WTA, and this makes them stick to the current route until one alternative route's monetary cost is enough to tempt them to choose.

The value of this study lies in verifying travelers' inertial route choices that resulted from asymmetric preference and estimating the values of parameters in Xu's [4] status quo-dependent route choice model. Moreover, by considering the asymmetric preference, we can explain why some travelers stick to a toll lane, even when there is an alternative road that is slightly better. However, there are some limitations in our experimental design. Firstly, we used the questionnaire survey data of 30 participants to estimate the values of parameters *λ*<sup>1</sup> and *λ*2. However, a behavioral experiment is needed to collect individuals' route-choice data to estimate the values of parameters *λ*<sup>1</sup> and *λ*2. Secondly, as noted earlier, our experimental design was conducted with a lack of interaction between the participants' choices. Thus, there is meaningful reason for designing a behavioral experiment, and to incorporate the asymmetric preference into traffic equilibrium modeling.

**Author Contributions:** Model construction, experiment design, software, analysis of results and model parameters estimation: K.L.; introduction writing, model parameters estimation: Y.X.

**Funding:** This work was supported by a national natural science foundation of China [No. 71571097] and a humanities and social sciences project of Anhui, China [No. SK2018A0025].

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
