A method for assessing the accessibility of urban ecological space based on ecosystem services

Abstract

The present invention relates to the technical field of urban ecological space accessibility assessment, and discloses an urban ecological space accessibility assessment method based on ecosystem services, including obtaining thematic element data; calculating the ecosystem service factor of the urban ecological space according to the thematic element data; generating a residential irregular grid according to the thematic element data; calculating the reachable range of each residential irregular grid according to the selected travel mode combined with the residential irregular grid and the thematic element data; calculating the attenuation coefficient of the urban ecological space position within the reachable range of each residential irregular grid according to the thematic element data; calculating the urban ecological space accessibility of each residential irregular grid based on ecosystem services according to the ecosystem service factor of the urban ecological space, the attenuation coefficient of the urban ecological space position in each residential irregular grid, and the thematic element data. The present invention can improve the calculation accuracy of the urban ecological space accessibility measurement.

Description

Technical Field

[0001] The present invention relates to the technical field of urban ecological space accessibility assessment, and in particular to an urban ecological space accessibility assessment method based on ecosystem services. Background Art

[0002] Accessibility is defined as the size of the opportunity for each node in the transportation network to interact with each other. The accessibility of urban space includes time, distance, space and other costs. At present, the research content on the accessibility of urban ecological space is extensive and has many applications, but few studies have explored the accessibility of residential areas to the use of urban ecological space and the benefits of the ecosystem provided to humans. Urban terrain conditions have prominent restrictions on travel activities (walking, cycling, driving, etc.). Under the conditions of complex terrain in mountainous cities, the actual path to the destination is longer and the psychological barriers are greater, which hinders the accessibility of urban ecological space in mountainous environments.

[0003] In summary, the impact of urban topography on the accessibility of urban ecological space cannot be ignored. At present, there are three major limitations in most urban ecological space accessibility technologies: first, few studies use residential areas as the starting point to explore the accessibility of residents in residential areas to surrounding urban ecological spaces; second, almost no research considers the ecosystem service function of urban ecological space in the urban ecological space accessibility model; third, few studies introduce terrain factors into the accessibility model (or the attenuation coefficient of the model). These limitations make the expression of the accessibility model of urban ecological space deviate from the actual needs of urban residents for urban ecological space, and it is difficult to serve as a reference for urban ecological space planning or residents' happiness. Summary of the invention

[0004] In view of the above-mentioned deficiencies in the prior art, the present invention provides an urban ecological space accessibility assessment method based on ecosystem services.

[0005] In order to achieve the above-mentioned object of the invention, the technical solution adopted by the present invention is:

[0006] An urban ecological space accessibility assessment method based on ecosystem services includes the following steps:

[0007] Get the thematic element data;

[0008] Calculate the ecosystem service factor of the urban ecological space based on the thematic element data;

[0009] Generate irregular residential grids based on thematic element data;

[0010] Calculate the reach of each residential irregular grid based on the selected travel mode combined with the residential irregular grid and thematic element data;

[0011] Based on the thematic element data, the attenuation coefficient of the urban ecological space position within the reach of each irregular residential grid is calculated;

[0012] According to the ecosystem service factor of urban ecological space, the attenuation coefficient of the urban ecological space position within the reach of each irregular residential grid, and the thematic element data, the accessibility of urban ecological space based on ecosystem services for each irregular residential grid is calculated.

[0013] Furthermore, the thematic element data includes:

[0014] The scope of the area to be evaluated, residences, road networks, urban ecological space, permanent population, land use data, nature reserves, normalized difference vegetation index, average daily actual evapotranspiration, digital elevation model data, meteorological data, surface runoff, negative ion concentration, tree height and soil physical and chemical properties data.

[0015] Furthermore, based on the thematic element data, the ecosystem service factors of the urban ecological space are calculated, including:

[0016] Calculate the service functions of various ecosystems based on the data of thematic elements; ecosystem service functions include air purification function, water purification function, climate regulation function, negative ion release function, water conservation function, flood regulation function, soil conservation function, and biodiversity maintenance function;

[0017] Based on the service function quantities of various ecosystems, the ecosystem service factor of the urban ecological space is calculated.

[0018] Furthermore, according to the service function of each type of ecosystem, the ecosystem service factor of the urban ecological space is calculated, which is as follows:

[0019]

[0020] Among them, ESfactor(p) represents the ecosystem service factor at pixel P, M(p,ec) is the standardized value of the ecosystem service function at pixel p, W(ec) represents the weight of the ecosystem service function, and ecn represents the quantity of ecosystem service function.

[0021] Furthermore, based on the thematic element data, an irregular residential grid is generated, including:

[0022] Divide the minimum rectangular range that completely surrounds the area to be evaluated into square grid data of the first side length, and superimpose the spatial vector data of the residence, and calculate the residential area ratio within each square grid of the set side length in turn;

[0023] The residential area ratios in the square grids of the first side length are judged in turn, the square grids with a residential area ratio of 0% are deleted, and the square grids with a residential area ratio less than the first ratio threshold are further divided into square grids of the second side length and the corresponding residential area ratios are calculated;

[0024] The residential area ratios in the square grids of the second side length are judged in turn, the square grids with a residential area ratio of 0% are deleted, and the square grids with a residential area ratio less than the first ratio threshold are further divided into square grids of the third side length and the corresponding residential area ratios are calculated;

[0025] The residential area ratios in the square grids of the third side length are judged in turn, and the square grids with a residential area ratio of 0% are deleted to obtain the remaining square grids of the first side length, the second side length, and the third side length;

[0026] The remaining square grids of the first side length, the second side length and the third side length are merged with adjacent elements, and all non-adjacent irregular grids after fusion are the results of the residential irregular grid division.

[0027] Furthermore, according to the selected travel mode, the reachable range of each residential irregular grid is calculated by combining the residential irregular grid and thematic element data, including:

[0028] Select a travel mode, overlay the slope and road network raster data, and calculate the true path distance for each road network raster pixel;

[0029] Based on the residential irregular grid data, road network data and urban ecological space raster data, taking the geometric center of each residential irregular grid as the starting point, determine the maximum reachable range of each residential irregular grid along the road network under the currently selected travel mode;

[0030] When the maximum reachable range of the residential irregular grid does not exceed the range of the current residential irregular grid, the reachable range of the current residential irregular grid is set to the urban ecological space grid pixels of all spatial resolutions within the residential irregular grid.

[0031] Furthermore, the actual path distance of each road grid pixel is calculated as follows:

[0032]

[0033] Among them, RL(k) represents the true path distance at pixel k, AR represents the pixel area, and θ(k) represents the slope at pixel k.

[0034] Furthermore, based on the thematic element data, the attenuation coefficient of the urban ecological space position within the reach of each residential irregular grid is calculated, specifically:

[0035]

[0036] Among them, ESSJ(d,m) represents the attenuation coefficient of pixel m in the irregular residential grid d, RLw(d,m) represents the shortest true path distance from the geometric center of the irregular residential grid d to pixel m under travel mode w, RLmax(d) represents the maximum value of the shortest true path distance from the geometric center of the irregular residential grid d to any urban ecological space pixel within its reach under travel mode w, and α(d,m) represents the terrain factor.

[0037] Furthermore, according to the ecosystem service factor of the urban ecological space, the attenuation coefficient of the urban ecological space position within the reach of each irregular residential grid, and the thematic element data, the urban ecological space accessibility of each irregular residential grid based on ecosystem services is calculated, specifically:

[0038]

[0039] Among them, ACC(d) represents the accessibility of urban ecological space based on ecosystem services of residential irregular grid d, ESSJ(d,m) represents the attenuation coefficient of pixel m in residential irregular grid d, ESfactor(d,m) represents the ecosystem service factor of pixel m in residential irregular grid d, PO(d,n) represents the number of permanent residents in pixel n in residential irregular grid d, dpm represents the number of urban ecological space grid pixels with a set spatial resolution in residential irregular grid d, dpn represents the number of permanent residents pixels with a set spatial resolution in residential irregular grid d, and n represents the permanent residents pixels with a set spatial resolution within the reachable range of residential irregular grid d.

[0040] The present invention has the following beneficial effects:

[0041] The present invention solves the limitations of current urban ecological space accessibility technology by dividing irregular residential grids, considering ecosystem service factors and travel modes, and introducing terrain factors into the attenuation coefficient. It can provide a decision-making basis for urban planning, urban renewal and other work. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Figure 1 This is a flowchart of an urban ecological space accessibility assessment method based on ecosystem services. DETAILED DESCRIPTION

[0043] The specific implementation modes of the present invention are described below so that those skilled in the art can understand the present invention. However, it should be clear that the present invention is not limited to the scope of the specific implementation modes. For those of ordinary skill in the art, as long as various changes are within the spirit and scope of the present invention as defined and determined by the attached claims, these changes are obvious, and all inventions and creations utilizing the concept of the present invention are protected.

[0044] The technical problems addressed by the present invention include:

[0045] 1. The accessibility model does not take into account the ecosystem service functions that the urban ecological space itself can provide.

[0046] Most current methods only consider the ratio of urban ecological space area to population (supply and demand). However, the ecosystem services (benefits to humans) that can be generated by urban ecological space of the same area are closely related to its cover type and quality. This will also make urban ecological space of the same area in different spaces have very different ecosystem service function quantities. Therefore, it is obviously unreasonable not to consider ecosystem services in the process of calculating the accessibility of urban ecological space.

[0047] 2. The division of demand grids does not reflect the spatial distribution characteristics of communities / parks / residential areas.

[0048] At present, most of the demand grids in the accessibility calculation process are based on the same size and shape (such as hexagons, squares, etc. of the same size). Often a community / park / residential area will be unevenly distributed in multiple grids, but the accessibility of the surrounding urban ecological space to the residents in the same community / park / residential area is the same / consistent. Therefore, it is more reasonable to divide the same community / park / residential area into the same irregular grid as much as possible.

[0049] 3. Insufficient accuracy of attenuation coefficient in accessibility model

[0050] At present, most existing technologies only consider the distance of the road / straight-line distance as a parameter when calculating the attenuation coefficient in the accessibility model, and do not consider the terrain factors on the road distance / straight-line distance. In fact, the terrain on the road has a greater impact on the accessibility of the urban ecological space of residents, that is, the longer the distance, the more complex the road terrain (the greater the terrain undulation), the more difficult the road is to walk, the lower the willingness of residents to go to the urban ecological space at the end of the road, and the lower the accessibility.

[0051] The solution adopted by the present invention includes:

[0052] 1. The present invention considers the ecosystem service factor in the process of calculating the accessibility of urban ecological space. Different areas to be evaluated select some or all ecosystem services and set weights according to actual conditions. It can reflect the benefits of the ecosystem to humans in the process of calculating the accessibility calculation index, which is more scientific and reasonable. It can solve the problem that the current urban ecological space accessibility calculation method does not consider the benefits of the urban ecological space itself to humans.

[0053] 2. The present invention targets the residential space distribution characteristics of the area to be evaluated, uses squares of different sizes to divide the grid and then fuses adjacent elements to obtain an irregular residential grid. This can greatly solve the problem that the grid division in the prior art does not reflect the spatial distribution characteristics of the community / park / residential area, and improve the calculation accuracy of the urban ecological space accessibility measurement.

[0054] 3. The present invention combines path distance and terrain factor (terrain undulation) in the calculation of the attenuation coefficient, reflecting that the longer the road distance and the more complex the road terrain (the greater the terrain undulation), the less beneficial the urban ecological space at the end of the road is to residents, which can solve the problem of insufficient accuracy of the attenuation coefficient in the prior art.

[0055] like Figure 1 As shown, an urban ecological space accessibility assessment method based on ecosystem services provided by an embodiment of the present invention includes the following steps S1 to S6:

[0056] S1. Obtain thematic element data;

[0057] In an optional embodiment of the present invention, the thematic element data acquired in step S1 includes:

[0058] The scope of the area to be evaluated, residences, road networks, urban ecological space, permanent population, land use data, nature reserves, normalized difference vegetation index, average daily actual evapotranspiration, digital elevation model data, meteorological data, surface runoff, negative ion concentration, tree height and soil physical and chemical properties data.

[0059] The subject element data of the area to be evaluated, the format description and the source description obtained in this embodiment are shown in Table 1:

[0060] Table 1

[0061]

[0062] This embodiment preprocesses various types of data as follows, which serves as the basis for subsequent steps:

[0063] 1. The data of urban ecological space, road network, nature reserve, etc. are first clipped to the scope of the area to be evaluated, and then converted from the spatial vector format and resampled to raster data. The spatial resolution of the newly generated raster is 10 m, and the spatial coordinates are consistent with the raster data after the data format conversion of the area to be evaluated;

[0064] 2. The grid data of permanent population, land use data, average daily actual evapotranspiration, surface runoff, digital elevation model data, normalized difference vegetation index, etc. are first clipped to the scope of the area to be evaluated, and then resampled to generate new raster data with a spatial resolution of 10 m. The spatial coordinates are consistent with the raster data after the data format conversion of the scope of the area to be evaluated;

[0065] 3. Calculate the slope, terrain relief and flow direction based on the generated digital elevation model data with a spatial resolution of 10 m;

[0066] 4. According to the daily rainfall in the meteorological data, the annual rainfall (the cumulative daily rainfall throughout the year) and the annual rainfall during heavy rain (the cumulative daily rainfall throughout the year is greater than or equal to 50 mm) of each point are counted. Then, the annual rainfall, annual rainfall during heavy rain, and soil physical and chemical properties data are spatially interpolated and resampled to a raster format. The spatial resolution of the newly generated raster is 10 m, and the range is the range of the area to be evaluated. The spatial coordinates are consistent with the raster data after the data format conversion of the area to be evaluated.

[0067] 5. According to meteorological data, count the number of days with a maximum daily temperature greater than or equal to 26 degrees Celsius and the number of days with a relative humidity less than or equal to 45% in the area to be evaluated; according to land use data, extract the total area of ​​water bodies in the area to be evaluated; according to the daily actual evapotranspiration raster data with a spatial resolution of 10 m, calculate the annual actual evapotranspiration raster with a spatial resolution of 10 m;

[0068] The above-mentioned operations such as format conversion, resampling, area extraction, clipping, spatial overlay analysis, calculation of flow direction, slope, terrain undulation, etc. can all be achieved through GIS software, such as ArcGis, QGIS, MAPGIS and other software.

[0069] S2. Calculate the ecosystem service factor of the urban ecological space based on the thematic element data;

[0070] In an optional embodiment of the present invention, step S2 calculates the ecosystem service factor of the urban ecological space according to the thematic element data, including:

[0071] Calculate the service functions of various ecosystems based on the data of thematic elements; ecosystem service functions include air purification function, water purification function, climate regulation function, negative ion release function, water conservation function, flood regulation function, soil conservation function, and biodiversity maintenance function;

[0072] According to the service function of various ecosystems, the ecosystem service factor of urban ecological space is calculated, which is as follows:

[0073]

[0074] Among them, ESfactor(p) represents the ecosystem service factor at pixel P, M(p,ec) is the standardized value of the ecosystem service function at pixel p, W(ec) represents the weight of the ecosystem service function, and ecn represents the quantity of ecosystem service function.

[0075] This embodiment first calculates the functional amount of various ecosystem services, and then calculates the ecosystem service factor according to the weights of various ecosystem services. Ecosystem services include: air purification, water purification, climate regulation, release of negative ions, water conservation, flood regulation, soil conservation, maintenance of biodiversity, etc. All or part of the ecosystem services can be selected according to the actual situation of the area to be evaluated, and the sum of the weights of the ecosystem services involved in the calculation is 1.

[0076] Calculate the air purification function:

[0077]

[0078] In the formula, p represents any 10 m spatial resolution pixel position in the urban ecological space of the area to be evaluated; ESA(p) represents the air purification function at pixel p, in t / year; i represents any air pollutant among sulfur dioxide, nitrogen oxides and particulate matter; QA(p,i) represents the purification amount of air pollutant i by the ecosystem at pixel p, in t / year; AR is the pixel area, in km 2 , that is, 10 m × 10 m = 100 m 2 = 0.0001 km 2 ; QAd(p,i) represents the purification amount per unit area of ​​air pollutant i by the ecosystem at pixel p, in units of t / (km 2 .Year).

[0079] The purification capacity of air pollutants per unit area of ​​the above-mentioned ecosystems can be referred to literature data or field monitoring data.

[0080] Calculate the water purification function:

[0081]

[0082] In the formula, p represents any 10 m spatial resolution pixel position in the urban ecological space of the area to be evaluated; ESW(p) represents the water purification function at pixel p, in t / year. If pixel p is not a water body, ESW(p) is equal to 0; i represents any water pollutant among chemical oxygen demand, total nitrogen and total phosphorus; WP(p,i) represents the total purification capacity of water pollutant i by the ecosystem at pixel p, in t / year; WPd(p,i) represents the purification capacity per unit area of ​​wetland for water pollutant i, in t / (km 2 .year); AR is the pixel area, unit is km 2 , that is, 10 m × 10 m = 100 m 2 = 0.0001 km 2 .

[0083] The above-mentioned unit area wetland purification capacity of water pollutants per unit area can be referred to literature data or field monitoring data.

[0084] Calculate the climate regulation capacity:

[0085]

[0086] Where p represents any 10 m spatial resolution pixel position of the green infrastructure in the area to be evaluated; ESC(p) represents the climate regulation function at pixel p, in kWh / year; Ept(p) represents the cooling function at pixel p, in kWh / year. If there is no vegetation or water at pixel p, Ept(p) is equal to 0; Ewe(p) represents the humidification function at pixel p, in kWh / year. If there is no water at pixel p, Ewe(p) is equal to 0; DET(p) represents the average daily actual evapotranspiration per unit area at pixel p, in mm / day; AR is the pixel area, in m 2 , that is, 10 m × 10 m = 100 m 2 ; D1 is the number of days when the air conditioner is open in the area to be evaluated, that is, the number of days when the maximum temperature is greater than or equal to 26°C, in days / year; ρ is the density of water, which is 1000 kg / m 3 ; q is the latent heat of volatilization at a standard atmospheric pressure and 100°C, with a value of 2257.2 kJ / kg; JW is the conversion coefficient from kJ to kilowatt-hour, with a value of 1 / 3600 kWh / kJ; r is the energy efficiency ratio of household air conditioners, dimensionless, with a value of 3; D2 is the number of days when the relative humidity in the area to be evaluated is less than or equal to 45%, in days / year; y is the power consumption of the humidifier to convert water into steam, with a value of 120 kWh / m 3 .

[0087] Calculate the amount of negative ion release:

[0088]

[0089] Where p represents any 10 m spatial resolution pixel position of the green infrastructure in the area to be evaluated; ESG(p) represents the negative ion release function at pixel p, in units of per year. If pixel p is not forest land, ESG(p) is equal to 0; QGd(p) represents the negative ion concentration at pixel p, in units of per cm 3 ; AR is the pixel area, unit is m 2 , that is, 10 m × 10 m = 100 m 2 ; H is the average tree height of the forest in the area to be evaluated, in meters; L is the negative ion lifetime, in minutes, which can be taken as 1 minute when there is no monitoring value.

[0090] Calculate water conservation function:

[0091]

[0092] Where p represents any 10 m spatial resolution pixel position of the green infrastructure in the area to be evaluated; ESR(p) represents the water conservation function at pixel p, with the unit of m 3 / year; Pre(p) is the annual rainfall at pixel p, in mm / year; R(p) is the surface runoff at pixel p, in mm / year; AET(p) is the actual annual evapotranspiration per unit area at pixel p, in mm / year; AR is the pixel area, in m 2 , that is, 10 m × 10 m = 100 m 2 .

[0093] Calculate the flood storage function:

[0094]

[0095] Where p represents any 10 m spatial resolution pixel position of the green infrastructure in the area to be evaluated; ESF(p) represents the flood storage function at pixel p, with the unit of m 3 / year; Ev(p) is the flood storage function of vegetation at pixel p, in m 3 / year, if there is no vegetation at pixel p, Ev(p) is equal to 0; Em(p) is the flood regulation function of the swamp at pixel p, the unit is m 3 / year, if pixel p is not a swamp, then Em(p) is equal to 0; Prb(p) represents the annual rainfall of the rainstorm at pixel p, in mm / year; pa and pb are constant coefficients, which are determined according to the land use type; AR is the pixel area, in m 2 , that is, 10 m × 10 m = 100 m 2 ; ph is the water storage capacity per unit area of ​​the swamp, and the reference value is 2.47 m3 / m 2 ; pi is the surface water height of swamp soil, and the reference value is 0.3 m.

[0096] The above pa and pb constant coefficients are shown in Table 2:

[0097] Table 2

[0098]

[0099] Calculate soil conservation capacity:

[0100]

[0101] Where p represents any 10 m spatial resolution pixel position of the green infrastructure in the area to be evaluated; ESS(p) represents the soil conservation function at pixel p, in units of t / year; R(p) represents the rainfall erosivity factor at pixel p, in units of MJ.mm / hm 2 .h.year; K(p) represents the soil erodibility factor at pixel p, in t.hm 2 .h / hm 2 .MJ.mm; LS(p) represents the slope length and slope factor at pixel p, dimensionless; C(p) represents the vegetation cover and management factor at pixel p, dimensionless; AR is the pixel area, unit is hm 2 , that is, 10 m × 10 m = 100 m 2 = 0.01 hm 2 .

[0102] The above rainfall erosivity factor is calculated as follows:

[0103]

[0104] Where: p represents any 10 m spatial resolution pixel position of the green infrastructure in the area to be evaluated; R(p) represents the rainfall erosivity factor at pixel p, in units of MJ.mm / hm 2 .h.year; Pre(p) represents the annual rainfall at pixel p, in mm / year.

[0105] The above soil erodibility factor is calculated as follows:

[0106]

[0107] Where: p represents any 10 m spatial resolution pixel position of the green infrastructure in the area to be assessed; K(p) represents the soil erodibility factor at pixel p, in t.hm 2 .h / hm 2.MJ.mm; KE(p) is the soil erodibility factor calculated by the erosion-productivity evaluation model at pixel p, in t.hm 2 .h / hm 2 .MJ.mm; ms(p) is the percentage of sand particles (0.05~2 mm) at pixel p, unit is %; mi(p) is the percentage of silt particles (0.002~0.05 mm) at pixel p, unit is %; mc(p) is the percentage of clay particles (<0.002 mm) at pixel p, unit is %; orgC(p) is the percentage of organic carbon at pixel p, unit is %.

[0108] The above slope length and slope factor are calculated as follows:

[0109]

[0110] Where: p represents any 10 m spatial resolution pixel position of the green infrastructure in the area to be evaluated; LS(p) represents the slope length factor at pixel p, dimensionless; L(p) is the slope length factor, dimensionless; S(p) is the slope factor, dimensionless; λ=λ coefficient×pixel resolution (m), the λ coefficient is determined according to the flow direction, if the direction value is 1, 4, 16, 64, the λ coefficient is 1, if the direction value is 2, 8, 32, 128, the λ coefficient is , the pixel resolution is 10 m; m(p) and β(p) are both dimensionless slope length factor parameters at pixel p; θ(p) is the slope at pixel p, in degrees.

[0111] The above vegetation cover and management factors are calculated as follows:

[0112]

[0113] Where: p represents any 10 m spatial resolution pixel position of the green infrastructure in the area to be assessed; C(p) represents the vegetation coverage and management factor at pixel p, dimensionless; FVC(p) represents the vegetation coverage at pixel p, unit is %; NDVI(p) represents the normalized vegetation index at pixel p, dimensionless; NDVIsoil is the value closest to the 5% cumulative percentage of the normalized vegetation index in the area to be assessed (or the minimum NDVI value in the area to be assessed); NDVIveg is the value closest to the 95% cumulative percentage of the normalized vegetation index in the area to be assessed (or the maximum NDVI value in the area to be assessed).

[0114] Calculate the amount of biodiversity maintenance function:

[0115]

[0116] Where p represents any 10 m spatial resolution pixel position of the green infrastructure in the area to be assessed; ESB(p) is the biodiversity function at pixel p, in m 2 / year, if pixel p is not a nature reserve, ESB(p) is equal to 0; AR is the pixel area, unit is m 2 , that is, 10 m × 10 m = 100 m 2 ;

[0117] Calculate the ecosystem service factor value:

[0118]

[0119] Where p represents any 10 m spatial resolution pixel position of the green infrastructure in the area to be evaluated; ESfactor(p) represents the ecosystem service factor at pixel P, which is dimensionless and ranges from 0 to 100. The larger the value, the greater the benefit of the ecosystem at pixel p to humans; ecn is the number of ecosystem service types, including air purification, water purification, climate regulation, negative ion release, water conservation, flood storage, soil conservation, and biodiversity maintenance. Some or all of them can be selected according to the actual situation of the area to be evaluated; ec represents any ecosystem service in ecn; M(p,ec) is the standardized value of the ecosystem service ec at pixel p, which is dimensionless and ranges from 0 to 100; W(ec) is the weight of the ecosystem service ec, and the sum of the weights of the ecn ecosystem service types is 1; ES(p,ec) represents the functional value of the ecosystem service ec at pixel p; ESminec represents the minimum functional value of the ecosystem service ec of the urban ecological space in the area to be evaluated; ESmaxec represents the maximum functional value of the ecosystem service ec of the urban ecological space in the area to be evaluated; pn represents the total number of 10 m spatial resolution pixels of the green infrastructure in the area to be evaluated.

[0120] The weights W(ec) of the above ecosystem services can be determined based on a questionnaire survey of residents in the area to be evaluated. If all ecosystem services such as air purification, water purification, climate regulation, negative ion release, water conservation, flood regulation, soil conservation, and biodiversity are selected, the weights are shown in Table 3:

[0121] Table 3

[0122]

[0123] S3, generating irregular residential grids based on thematic element data;

[0124] In an optional embodiment of the present invention, step S3 generates an irregular residential grid according to the thematic element data, including:

[0125] Divide the minimum rectangular range that completely surrounds the area to be evaluated into square grid data of the first side length, and superimpose the spatial vector data of the residence, and calculate the residential area ratio within each square grid of the set side length in turn;

[0126] The residential area ratios in the square grids of the first side length are judged in turn, the square grids with a residential area ratio of 0% are deleted, and the square grids with a residential area ratio less than the first ratio threshold are further divided into square grids of the second side length and the corresponding residential area ratios are calculated;

[0127] The residential area ratios in the square grids of the second side length are judged in turn, the square grids with a residential area ratio of 0% are deleted, and the square grids with a residential area ratio less than the first ratio threshold are further divided into square grids of the third side length and the corresponding residential area ratios are calculated;

[0128] The residential area ratios in the square grids of the third side length are judged in turn, and the square grids with a residential area ratio of 0% are deleted to obtain the remaining square grids of the first side length, the second side length, and the third side length;

[0129] The remaining square grids of the first side length, the second side length and the third side length are merged with adjacent elements, and all non-adjacent irregular grids after fusion are the results of the residential irregular grid division.

[0130] In this embodiment, the minimum rectangular range that completely surrounds the area to be evaluated is first divided into square grid data with a side length of 500 m, and the spatial vector data of the residential buildings is superimposed, and the residential area ratio in each square grid with a side length of 500 m is calculated one by one:

[0131]

[0132] Where i represents a square grid with a side length of 500 m in the area to be assessed, and the value of i ranges from 1 to n5; n5 is the total number of 500 m square grids with residential buildings in the area to be assessed; BP500(i) represents the proportion of residential area in grid i, in units of %; B500(i) is the area of ​​residential buildings in grid i, in units of m 2 ; 500×500 represents the area of ​​grid i, in m 2 , that is, 500 m × 500 m = 250000 m 2 .

[0133] This embodiment then determines the residential area ratio of each square grid with a side length of 500 m in the area to be evaluated one by one. If the residential area ratio BP500(i) in the grid is 0%, the grid is deleted. If the building area ratio BP500(i) of the grid is less than 70% and greater than 0%, the grid is further divided into 4 square grids with a side length of 250 m and the building area ratios corresponding to the 4 grids are calculated:

[0134]

[0135] Where j represents a square grid with a side length of 250 m, and j ranges from 1 to 4; BP250(j) represents the proportion of residential area in grid j, in units of %; B250(j) represents the residential area in grid j, in units of m 2 ; 250×250 represents the area of ​​grid j, in m 2 , that is, 250 m × 250 m = 62500 m 2 .

[0136] This embodiment then determines the residential area ratio of the square grids with a side length of 250 m in the area to be evaluated one by one. If the building area ratio in the grid is 0%, the grid is deleted. If the residential area ratio in the grid is less than 70% and greater than 0%, the grid is divided into 25 square grids with a side length of 50 m and the residential area ratio corresponding to the 25 grids is calculated:

[0137]

[0138] Where k represents a square grid with a side length of 50 m, and the value of k ranges from 1 to 25; BP50(k) represents the proportion of residential area in grid k, in units of %; B50(k) represents the residential area in grid k, in units of m 2 ; 50×50 represents the area of ​​grid k, in m 2 , that is, 50 m × 50 m = 2500 m 2 .

[0139] Finally, this embodiment determines the residential area ratio of each square grid with a side length of 50 m in the area to be evaluated one by one, and if the building area ratio in the grid is 0%, the grid is deleted.

[0140] All the remaining 500 m, 250 m or 50 m square grids are merged with adjacent elements, and all non-adjacent irregular grids obtained after fusion are the irregular grid division results of the residential areas to be evaluated.

[0141] S4, calculating the reachable range of each residential irregular grid based on the selected travel mode combined with the residential irregular grid and thematic element data;

[0142] In an optional embodiment of the present invention, step S4 calculates the reachable range of each irregular residential grid according to the selected travel mode combined with the irregular residential grid and thematic element data, including:

[0143] Select a travel mode, overlay the slope and road network raster data, and calculate the true path distance for each road network raster pixel;

[0144] Based on the residential irregular grid data, road network data and urban ecological space raster data, taking the geometric center of each residential irregular grid as the starting point, determine the maximum reachable range of each residential irregular grid along the road network under the currently selected travel mode;

[0145] When the maximum reachable range of the residential irregular grid does not exceed the range of the current residential irregular grid, the reachable range of the current residential irregular grid is set to the urban ecological space grid pixels of all spatial resolutions within the residential irregular grid.

[0146] In this embodiment, a travel mode (denoted as w) is first selected, such as walking, non-motorized vehicles (bicycles, electric bicycles, motorcycles) or motor vehicles. The grid data with a spatial resolution of 10 m after slope and road network preprocessing is superimposed to calculate the actual path distance of each road network grid pixel:

[0147]

[0148] Where k is any 10 m spatial resolution pixel position in the road network data of the area to be evaluated; RL(k) represents the real path distance at pixel k, in m; AR is the pixel area, in m 2 , that is, 10 m × 10 m = 100 m 2 ; θ(k) is the slope at pixel k, in degrees.

[0149] This embodiment then superimposes the irregular residential grid data, road network data, and the 10 m spatial resolution raster data after urban ecological space preprocessing, and takes the geometric center of each irregular residential grid as the starting point to determine the maximum reachable range of each grid along the road network under the currently selected travel mode. Within the reachable range, starting from the geometric center of the grid to any location within the reachable range does not exceed the maximum travel time for leisure and recreation under the currently selected travel mode:

[0150]

[0151] Where d represents any residential irregular grid in the area to be evaluated; AG(d) represents the reachable range of the residential irregular grid d; m represents any urban ecological space grid pixel with a spatial resolution of 10 m within the reachable range AG(d) of the residential irregular grid d; TO(d,m) represents the minimum time from the geometric center of the residential irregular grid d to pixel m, in minutes; Twmax represents the maximum time for leisure and recreational travel when travel mode w is selected, in minutes; n represents the type of urban road, 1 to 4 are expressway, main road, secondary road and branch road respectively; RLwn(d,m) represents the length of urban road type n in the shortest true path distance from the geometric center of the residential irregular grid d to pixel m, in km; wsn represents the average speed of urban road type n when travel mode w is selected, in km / h.

[0152] When the above reachable range does not exceed the range of the current grid, the reachable range of the grid is set to all urban ecological space grid pixels with a spatial resolution of 10 m in the grid. The above RLwn(d,m), that is, the length of urban road type n in the shortest real path distance from the geometric center of the residential irregular grid d to pixel m, can be implemented in GIS software (such as ArcGIS, QGIS, MAPGIS, etc.), using the road network data of the real path distance to extract the shortest real path distance from the geometric center of the residential irregular grid d to pixel m, and then extract the length of urban road type n on this path according to the road network type attribute.

[0153] The above travel modes and their corresponding maximum travel time for leisure and recreation and average speed on different types of urban roads are shown in Table 4 (can also be obtained through relevant departments or field surveys based on the actual situation of the area to be evaluated):

[0154] Table 4

[0155]

[0156] S5. Based on the thematic element data, calculate the attenuation coefficient of the urban ecological space location within the reach of each irregular residential grid;

[0157] In an optional embodiment of the present invention, step S5 calculates the attenuation coefficient of the urban ecological space position within the reachable range of each irregular residential grid according to the thematic element data, specifically:

[0158]

[0159] Wherein, d represents any residential irregular grid in the area to be evaluated; m represents any urban ecological space grid pixel with a spatial resolution of 10 m within the reachable range of the residential irregular grid d; ESSJ(d,m) represents the attenuation coefficient of pixel m in the residential irregular grid d, and its value range is between 0 and 1. The smaller the value, the smaller the distance from the resident to pixel m; RLw(d,m) represents the shortest true path distance from the geometric center of the residential irregular grid d to pixel m under travel mode w, in km; RLmax(d) represents the maximum value of the shortest true path distance from the geometric center of the residential irregular grid d to any urban ecological space pixel within its reach under travel mode w, in km; α(d,m) represents the terrain factor, that is, the terrain undulation on the shortest true path from the geometric center of the residential irregular grid d to pixel m, that is, the maximum elevation value on the path minus the minimum elevation value on the path, in m.

[0160] S6. Calculate the accessibility of urban ecological space based on ecosystem services for each irregular residential grid based on the ecosystem service factor of the urban ecological space, the attenuation coefficient of the urban ecological space location within the reach of each irregular residential grid, and the thematic element data.

[0161] In an optional embodiment of the present invention, step S6 calculates the urban ecological space accessibility of each irregular residential grid based on ecosystem services according to the ecosystem service factor of the urban ecological space, the attenuation coefficient of the urban ecological space position within the reachable range of each irregular residential grid, and the thematic element data, specifically:

[0162]

[0163] In the formula, d represents any irregular residential grid divided into the area to be evaluated; ACC(d) represents the urban ecological space accessibility based on ecosystem services of grid d. The larger the value, the better the urban ecological space accessibility based on ecosystem services of the grid, and the greater the ecosystem benefits that residents of the grid can enjoy when traveling; m represents any 10 m spatial resolution urban ecological space grid pixel within the reachable range of the irregular residential grid d; dpm is the number of 10 m spatial resolution urban ecological space grid pixels in the irregular residential grid d; ESSJ(d,m) represents the ecosystem service attenuation coefficient of pixel m in the irregular residential grid d; ESfactor(d,m) is the ecosystem service factor of pixel m in the irregular residential grid d, dimensionless, and the data is extracted through the ecosystem service factor ESfactor(p) at pixel P in the area to be evaluated; n represents any 10 m spatial resolution permanent population pixel within the reachable range of the irregular residential grid d; dpn is the number of 10 m spatial resolution urban ecological space grid pixels within the irregular residential grid d. m is the number of pixels with permanent population at spatial resolution; PO(d,n) is the number of permanent population of pixel n in grid d, in persons.

[0164] The present invention considers the benefits of various ecosystems to humans when calculating the accessibility of urban ecological space. By dividing the irregular residential grid and considering the terrain factors in the attenuation coefficient calculation process, the accessibility of urban ecological space based on ecosystem services is obtained, which can solve the following problems:

[0165] 1. This patent considers the ecosystem service factor in the process of calculating the accessibility of urban ecological space, solving the problem that the current urban ecological space accessibility technology does not reflect the benefits of the urban ecological space itself to humans. At present, most urban ecological space accessibility technologies only consider the area and population ratio (supply and demand) of the urban ecological space, and do not take into account the ecosystem service functions (benefits to humans) that the urban ecological space can generate. This patent considers the ecosystem service factor in the process of calculating the accessibility of urban ecological space, selects part or all of the ecosystem service function quantities according to the actual situation in different areas to be evaluated, and calculates the ecosystem service factor, which can reflect the benefits of the ecosystem to humans in the process of calculating the accessibility of urban ecological space, and is more scientific and reasonable.

[0166] 2. This patent solves the problem of dividing residential areas into multiple grids of the same size in the prior art, which destroys the spatial layout of residential areas, by dividing the residential areas into irregular grids. At present, the grid division of most accessibility calculation processes is based on the same size and shape. The divided grids often divide a residential area into multiple grids, which destroys the morphological pattern of the spatial distribution of the residential area itself. This patent divides the residential area into irregular grids with rectangular grids of different specifications and integrates the field elements to obtain the results of the residential area division, which can greatly ensure the integrity of the spatial layout of the residential area and improve the calculation accuracy of accessibility.

[0167] 3. This patent combines the plane path distance and terrain factors in the calculation of the attenuation coefficient, which solves the problem of insufficient accuracy in the calculation of the attenuation coefficient in the prior art. The actual terrain undulations of the path have a greater impact on the residents' willingness to go to the urban ecological space. That is, the longer the distance and the more complex the path terrain (the greater the elevation difference), the less convenient it is to walk. The lower the residents' willingness to go to the urban ecological space at the end of the path, which means that the accessibility of the urban ecological space at the end of the path is lower. However, most of the current calculations of the attenuation coefficient do not take into account the actual terrain undulations of the path. This patent combines the plane path distance and terrain factors in the calculation of the attenuation coefficient, which can improve the accuracy of the attenuation coefficient calculation.

[0168] The present invention is described with reference to flowcharts and / or block diagrams of methods, devices (systems) and computer program products according to embodiments of the present invention. It should be understood that each process and / or block in the flowchart and / or block diagram, as well as the combination of processes and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 A process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0169] These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer-readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 A process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

[0170] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for implementing the process. Figure 1 A process or multiple processes and / or boxes Figure 1 The steps for the functions specified in one or more boxes.

[0171] The present invention uses specific embodiments to illustrate the principles and implementation methods of the present invention. The description of the above embodiments is only used to help understand the method of the present invention and its core idea. At the same time, for those skilled in the art, according to the idea of ​​the present invention, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as a limitation on the present invention.

[0172] Those skilled in the art will appreciate that the embodiments described herein are intended to help readers understand the principles of the present invention, and should be understood that the protection scope of the present invention is not limited to such specific statements and embodiments. Those skilled in the art can make various other specific variations and combinations that do not deviate from the essence of the present invention based on the technical revelations disclosed by the present invention, and these variations and combinations are still within the protection scope of the present invention.

Claims

1. A method for assessing the accessibility of urban ecological space based on ecosystem services, characterized in that: The following steps are involved: Obtain thematic element data including digital elevation model data, road network data and ecological environment monitoring data, and pre-process to obtain raster data of urban ecological space; Calculate the ecosystem service factor of the urban ecological space based on the raster data of the urban ecological space; Generate irregular residential grid data based on residential spatial distribution characteristics; The shortest true path distance of each grid is calculated based on the selected travel mode and the grid data of the urban ecological space, and the reachable range of each residential irregular grid is determined by combining the residential irregular grid data and the road network data; Based on the shortest true path distance of each grid combined with the terrain factor, the attenuation coefficient of the urban ecological space position within the reach of each residential irregular grid is calculated; According to the ecosystem service factor of the urban ecological space, the attenuation coefficient of the urban ecological space position within the reach of each irregular residential grid, and the raster data of the urban ecological space, the urban ecological space accessibility of each irregular residential grid based on ecosystem services is calculated, specifically: Among them, ACC(d) represents the accessibility of urban ecological space based on ecosystem services of residential irregular grid d, ESSJ(d,m) represents the attenuation coefficient of pixel m in residential irregular grid d, ESfactor(d,m) represents the ecosystem service factor of pixel m in residential irregular grid d, PO(d,n) represents the number of permanent residents of pixel n in residential irregular grid d, dpm represents the number of urban ecological space grid pixels with a set spatial resolution in residential irregular grid d, dpn represents the number of permanent residents pixels with a set spatial resolution in residential irregular grid d, and n represents the permanent residents pixels with a set spatial resolution within the reachable range of residential irregular grid d.

2. The urban ecological space accessibility assessment method based on ecosystem services according to claim 1 is characterized in that: Thematic element data include: The scope of the area to be evaluated, residences, road networks, urban ecological space, permanent population, land use data, nature reserves, normalized difference vegetation index, average daily actual evapotranspiration, digital elevation model data, meteorological data, surface runoff, negative ion concentration, tree height and soil physical and chemical properties data.

3. The urban ecological space accessibility assessment method based on ecosystem services according to claim 1 is characterized in that: Based on the raster data of urban ecological space, the ecosystem service factors of urban ecological space are calculated, including: Calculate the service functions of various ecosystems based on the data of thematic elements; ecosystem service functions include air purification function, water purification function, climate regulation function, negative ion release function, water conservation function, flood regulation function, soil conservation function, and biodiversity maintenance function; Based on the service function quantities of various ecosystems, the ecosystem service factor of the urban ecological space is calculated.

4. The urban ecological space accessibility assessment method based on ecosystem services according to claim 3 is characterized in that: According to the service function of various ecosystems, the ecosystem service factor of urban ecological space is calculated, which is as follows: Among them, ESfactor(p) represents the ecosystem service factor at pixel p, M(p,ec) is the standardized value of the ecosystem service function at pixel p, W(ec) represents the weight of the ecosystem service function, and ecn represents the quantity of ecosystem service function.

5. The urban ecological space accessibility assessment method based on ecosystem services according to claim 1 is characterized in that: According to the residential spatial distribution characteristics, irregular residential grid data is generated, including: Divide the minimum rectangular range that completely surrounds the area to be evaluated into square grid data of the first side length, and superimpose the spatial vector data of the residence, and calculate the residential area ratio within each square grid of the set side length in turn; The residential area ratios in the square grids of the first side length are judged in turn, the square grids with a residential area ratio of 0% are deleted, and the square grids with a residential area ratio less than the first ratio threshold are further divided into square grids of the second side length and the corresponding residential area ratios are calculated; The residential area ratios in the square grids of the second side length are judged in turn, the square grids with a residential area ratio of 0% are deleted, and the square grids with a residential area ratio less than the first ratio threshold are further divided into square grids of the third side length and the corresponding residential area ratios are calculated; The residential area ratios in the square grids of the third side length are judged in turn, and the square grids with a residential area ratio of 0% are deleted to obtain the remaining square grids of the first side length, the second side length, and the third side length; The remaining square grids of the first side length, the second side length and the third side length are merged with adjacent elements, and all non-adjacent irregular grids after fusion are the results of the residential irregular grid division.

6. The urban ecological space accessibility assessment method based on ecosystem services according to claim 1 is characterized in that: The shortest true path distance of each grid is calculated based on the selected travel mode and the grid data of the urban ecological space, and the reachable range of each irregular residential grid is determined by combining the irregular residential grid data and the road network data, including: Select a travel mode, overlay the slope and road network raster data, and calculate the true path distance for each road network raster pixel; Based on the residential irregular grid data, road network data and urban ecological space raster data, taking the geometric center of each residential irregular grid as the starting point, determine the maximum reachable range of each residential irregular grid along the road network under the currently selected travel mode; When the maximum reachable range of the residential irregular grid does not exceed the range of the current residential irregular grid, the reachable range of the current residential irregular grid is set to the urban ecological space grid pixels of all spatial resolutions within the residential irregular grid.

7. The urban ecological space accessibility assessment method based on ecosystem services according to claim 6 is characterized in that: The calculation method for calculating the true path distance of each road network grid pixel is: Among them, RL(k) represents the true path distance at pixel k, AR represents the pixel area, and θ(k) represents the slope at pixel k.

8. The urban ecological space accessibility assessment method based on ecosystem services according to claim 1 is characterized in that: According to the shortest true path distance of each grid combined with the terrain factor, the attenuation coefficient of the urban ecological space position within the reach of each irregular residential grid is calculated, specifically: Among them, ESSJ(d,m) represents the attenuation coefficient of pixel m in the irregular residential grid d, RLw(d,m) represents the shortest true path distance from the geometric center of the irregular residential grid d to pixel m under travel mode w, RLmax(d) represents the maximum value of the shortest true path distance from the geometric center of the irregular residential grid d to any urban ecological space pixel within its reach under travel mode w, and α(d,m) represents the terrain factor.

Images