Road influence domain division method

Abstract

The invention discloses a road influence domain division method, which comprises the following steps of S1, acquiring Landsat remote sensing image data of a research area, performing preprocessing, and extracting ecological indexes; s2, coupling the extracted ecological indexes by using a principal component analysis method, extracting a first principal component, and calculating a remote sensing ecological index; s3, dividing the ecological environment quality into different grades according to the remote sensing ecological index, performing reclassification, and constructing a cost grid map; s4, calculating the Euclidean distance from each pixel in the cost grid map to the nearest road and the cost distance considering the landscape resistance, and analyzing the change rule of the remote sensing ecological index along with the cost distance; s5, extracting abrupt change points by using the change rule of the remote sensing ecological index along with the cost distance, determining a road influence threshold value, and delimiting a road influence domain; according to the method, the road influence threshold value is determined by constructing the cost grid map and combining the Euclidean distance, the cost distance and mutation point detection, so that the road influence domain analysis is realized.

Description

Technical Field

[0001] The present invention relates to the technical field of road influence domain division, and more particularly to a road influence domain division method. Background Art

[0002] The road effect zone (REZ) refers to the area surrounding a road, whose ecology, environment, and socioeconomic status are significantly impacted by road construction and operation. This area typically includes the road itself, its shoulders, slopes, drainage systems, and surrounding ecosystems and communities within a certain range. The specific extent of the REZ depends on factors such as road type, traffic volume, topography, and ecosystem sensitivity. Current research often uses experience or regulations to define a fixed width (e.g., 100 m, 500 m, 1000 m) as the road REZ. This fails to fully account for the spatial heterogeneity of environmental and human factors, such as topography, vegetation cover, and distance from settlements, which vary across space and have inconsistent effects on road ecological impacts. Summary of the Invention

[0003] The purpose of the present invention is to provide a road influence domain division method, which realizes the analysis of road influence domain by constructing a cost grid map and combining Euclidean distance, cost distance and mutation point detection to determine the road influence threshold.

[0004] To achieve the above object, the present invention provides the following technical solutions:

[0005] A road influence domain division method includes the following steps:

[0006] S1. Obtain Landsat remote sensing image data of the study area, perform preprocessing, and extract ecological indicators;

[0007] S2. Using principal component analysis to couple the extracted ecological indicators, extract the first principal component, and use the first principal component to calculate the remote sensing ecological index;

[0008] S3. Divide the ecological environment quality into different levels according to the remote sensing ecological index, reclassify it, and construct a cost grid map;

[0009] S4. Calculate the Euclidean distance from each pixel in the cost grid map to the nearest road and the cost distance considering landscape resistance, and analyze the variation pattern of remote sensing ecological index with cost distance;

[0010] S5. Based on the mutation point detection method, the mutation points are extracted using the change pattern of remote sensing ecological index with cost distance, the road impact threshold is determined, and the road impact domain is delineated.

[0011] Furthermore, in said S1, the preprocessing includes: radiation correction, atmospheric correction, and geometric correction.

[0012] Furthermore, in said S1, the ecological indicators include: greenness, humidity, heat, and dryness.

[0013] Furthermore, the greenness is the normalized vegetation index, and the calculation formula is:

[0014]

[0015] Among them, R red and R nir They represent the reflectance of the red band and near-infrared band of the Landsat image respectively.

[0016] Furthermore, the humidity refers to the surface humidity, which is obtained by extracting the humidity component based on Landsat TM and OLI data. The calculation method is as follows:

[0017] WET(TM)=0.0315R blue +0.2021R green +0.3102R red +0.1594R nir

[0018] -0.6706R nir1 -0.6109R nir2

[0019] WET(OLI)=0.1511R blue +0.1973R green +0.3283R red +0.3407R nir -0.7117R nir1 -0.4599R nir2

[0020] Where: R blue 、R green 、R red 、R nir 、R nir1 、R nir2 They are the reflectances of the blue, green, red, near-infrared, shortwave range 1, and shortwave range 2 bands respectively.

[0021] Furthermore, the heat refers to the thermal environment conditions, specifically:

[0022] Land surface temperature was retrieved from the thermal band of Landsat imagery and processed using a radiative transfer equation algorithm to represent thermal indicators;

[0023] The radiation transfer equation is expressed as follows:

[0024] L λ =[εB(T s )+(1-ε)L ↓ ] / +L ↓

[0025] Among them, L λ is the radiation at the sensor or the top of the atmosphere, ε is the surface emissivity, B(T s ) is the temperature T given by Planck's law s Blackbody radiation, L ↓ and L ↑ are the downwelling and upwelling atmospheric radiations, respectively, and τ is the total atmospheric transmittance between the land surface and the sensor;

[0026] B(T S ) is:

[0027] B(T s )=[Lλ-L ↑ -τ(1-ε)L ↓ ] / τε

[0028] According to Planck's law, we can get T S :

[0029] T s =K2 / ln(K1 / B(T s )+1)

[0030] Among them, K1 and K2 are constants related to the wavelength of the selected band.

[0031] Furthermore, the dryness is specifically the average value of the normalized building index and the bare soil index, expressed as:

[0032]

[0033] SI=[(R mir1 +R red )-(R nir +R blue )] / [(R mir1 +R red )+(R nir +R blue )]

[0034]

[0035] Among them, R red 、R blue 、R green 、R nir 、R mir1are the reflectances of the red band, blue band, green band, near-infrared band and shortwave infrared 1 band, SI is the bare soil index, and IBI is the normalized building index.

[0036] Furthermore, in S2, the extracted ecological indicators are coupled using principal component analysis to extract the first principal component, and the remote sensing ecological index is calculated using the first principal component, specifically:

[0037] Normalizing the ecological indicators, calculating the first principal component based on GEE using feature analysis, and calculating the initial remote sensing ecological indicators using the first principal component; normalizing the initial remote sensing ecological indicators to obtain a final remote sensing ecological index;

[0038] Among them, the expression of the initial remote sensing ecological index RSEI0 is:

[0039] RSEI0=1-PC1[f(NDVI,WET,NDBSI,LST)]

[0040] Among them, PC1 is the first principal component of ecological indicators.

[0041] Furthermore, in S4, the Euclidean distance from each pixel in the cost grid map to the nearest road and the cost distance considering landscape resistance are calculated; wherein the cost distance considering landscape resistance is specifically calculated as follows: the average value of the remote sensing ecological index of different cumulative cost distances is calculated as follows:

[0042]

[0043] Among them, n k is the total number of pixels in the partition, Z k is the division interval of cost distance, RSEI xy is the remote sensing ecological index value within the pixel.

[0044] According to the specific embodiments provided by the present invention, the present invention discloses the following technical effects:

[0045] The present invention uses remote sensing technology and objective methods such as principal component analysis to extract and construct remote sensing ecological indexes, reducing interference from human factors and improving the accuracy and reliability of assessments; through standardization and grading, the ecological environment quality of different regions is made comparable, further improving the accuracy of assessments; constructing a cost grid realizes spatial visualization of ecological environment quality, which is convenient for intuitive understanding and analysis; using Euclidean distance and cost distance calculation, the degree and scope of the impact of roads on the surrounding ecological environment can be intuitively displayed; through the mutation point detection method, the mutation points of the remote sensing ecological index changing with the cost distance are extracted, providing a scientific basis for determining the road impact threshold; the determination of the road impact threshold provides an important reference basis for road planning and ecological protection; combining the cost distance threshold, the road impact domain and the Euclidean distance analysis, the road impact range can be accurately determined, providing decision support for road construction, ecological protection and urban planning; the road impact domain division method provided by the present invention breaks through the limitation of the traditional fixed width of the road impact domain, and its width changes with the changes in the objective geographical space environment, and the result is more scientific. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are merely embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the provided drawings without paying any creative work.

[0047] The road influence domain division method of the present invention will be further described below with reference to the accompanying drawings;

[0048] Figure 1 It is an overall flow chart of the road influence domain division method of the present invention;

[0049] Figure 2 It is a spatial distribution map of the initial remote sensing ecological index in the road influence domain division of the present invention;

[0050] Figure 3 It is a spatial distribution map after raster reclassification in the present invention;

[0051] Figure 4 It is a spatial distribution map for performing Euclidean distance calculation on the total cost grid map in the present invention;

[0052] Figure 5 It is a spatial distribution map for cost distance calculation of the total cost grid map in the present invention;

[0053] Figure 6 It is a curve diagram of the mutation point detection of the remote sensing ecological index and the cost distance in the present invention;

[0054] Figure 7 It is a diagram of the road influence domain division result in the road influence domain division of the present invention. DETAILED DESCRIPTION

[0055] The following embodiments of the present invention are described in further detail with reference to the accompanying drawings and examples. The following examples are used to illustrate the present invention but are not intended to limit the scope of the present invention.

[0056] In order to better understand the purpose, structure and function of the present invention, the present invention is further described in detail below with reference to the accompanying drawings.

[0057] like Figure 1 As shown, the present invention provides a method for dividing a road influence domain, comprising the following steps:

[0058] S1. Obtain Landsat remote sensing image data of the study area, perform preprocessing, and extract ecological indicators;

[0059] S2. Using principal component analysis to couple the extracted ecological indicators, extract the first principal component, and use the first principal component to calculate the remote sensing ecological index;

[0060] S3. Divide the ecological environment quality into different levels according to the remote sensing ecological index, reclassify it, and construct a cost grid map;

[0061] In this example, a key step in calculating the minimum cumulative value based on the minimum cumulative resistance model is the creation of a cost grid, which is used to quantify the resistance of different landform types in the landscape to the road's ecological impact. The purpose of creating a cost grid is to assign a cost value to each grid cell, reflecting the difficulty of the road's impact moving through the area. Several cost allocation strategies can be used:

[0062] 1) Equal allocation. All feature types are assigned the same cost value. This strategy assumes that all landscape types have the same resistance to the road's ecological impact process and is suitable for situations where detailed ecological data is lacking or the landscape heterogeneity in the study area is low.

[0063] 2) Favorable allocation. Land features that are favorable to the road's ecological impact process are assigned lower cost values, while land features that are unfavorable to migration are assigned higher cost values. For example, in this case, areas with poor ecological quality may be assigned lower costs, while areas with good ecological quality may be assigned higher costs.

[0064] 3) Disadvantageous allocation. Lower cost values ​​are assigned to feature types that are unfavorable to the road's ecological impact process, while higher cost values ​​are assigned to feature types that are favorable to migration. This strategy is often used for special research purposes, such as assessing the ability of impacts to migrate in unfavorable environments.

[0065] According to the actual needs of this embodiment, a more favorable allocation system is selected. First, the remote sensing ecological index RSEI is divided into 5 levels at intervals of 0.2, from high to low: good, better, medium, poor and bad. It is further reclassified as: good → 10, better → 9, medium → 3, poor → 2 and poor → 1, thereby constructing a cost grid; Figure 2 and Figure 3 shown.

[0066] S4. Calculate the Euclidean distance from each pixel in the cost grid map to the nearest road and the cost distance considering landscape resistance, and analyze the variation pattern of remote sensing ecological index with cost distance;

[0067] In this embodiment, in ArcGIS, Euclidean distance and cost distance are calculated, which are used to calculate the straight-line distance from each pixel to the nearest road and the weighted distance considering landscape resistance (cost grid), respectively. Figure 4 and Figure 5 shown.

[0068] S5. Based on the mutation point detection method, the mutation points are extracted using the change pattern of remote sensing ecological index with cost distance, the road impact threshold is determined, and the road impact domain is delineated;

[0069] In this embodiment, the road impact domain is delineated based on the selected cost distance 3000 as the road impact threshold. Figure 7 The green part in ). Figure 7 It can be seen that the road influence domain is an irregular strip domain. The influence domain of the northern section of the road is wider, while the influence domain of the southern section of the road is narrower, ranging from 300 to 1000m. Figure 2 and Figure 3 As can be seen, the ecological environment on both sides of the northern section of the road is poor, with a larger impact area, while the ecological environment on both sides of the southern section is better, with a smaller impact area. This shows that the road influence domain demarcation method proposed in this invention is consistent with objective facts, providing new ideas and methods for demarcating road influence areas, and providing a scientific basis for road planning and road ecological protection.

[0070] In S1, Landsat remote sensing image data is acquired and preprocessed to ensure data quality, and ecological indicators are extracted; the preprocessing includes radiation correction, atmospheric correction, and geometric correction.

[0071] In the above-mentioned S1, Landsat remote sensing image data is obtained, pre-processed, and ecological indicators are extracted; wherein the ecological indicators include: greenness, humidity, heat, and dryness.

[0072] Greenness (NDVI): Normalized Difference Vegetation Index, reflecting vegetation coverage. NDVI values ​​range from -1 to 1. Generally speaking, the better the vegetation growth, the higher the NDVI. Therefore, NDVI can reflect the advantages and disadvantages of the ecological environment. The calculation formula of NDVI is as follows

[0073]

[0074] R red and R nir They represent the reflectance of the red band and near-infrared band of the Landsat image respectively.

[0075] Wetness: The moisture component (WET) extracted through the Tasseled Cap Transform reflects surface humidity. The calculation method for extracting WET based on Landsat TM and OLI data is as follows:

[0076] WET(TM)=0.0315R blue +0.2021R green +0.3102R red +0.1594R nir

[0077] -0.6706R nir1 -0.6109R nir2

[0078] WET(OLI)=0.1511R blue +0.1973R green +0.3283R red +0.3407R nir -0.7117R nir1 -0.4599R nir2

[0079] Where: R blue 、R green 、R red 、R nir 、R nir1 、R nir2 They are the reflectance of the blue, green, red, near-infrared, shortwave range 1, and shortwave range 2 bands, respectively. Before extraction, the water body index MNDWI is used to mask the water body so that WET reflects the actual land surface moisture conditions.

[0080] Thermal (LST): Land surface temperature, reflecting the thermal environment. The LST is retrieved from the thermal bands of Landsat satellite images (Landsat TM band 6, Landsat OLI band 10) and processed using the Radiative Transfer Equation (RTE) algorithm to represent the thermal indicator.

[0081] The formula of the radiative transfer equation (RTE) algorithm is as follows:

[0082] L λ =[εB(T s )+(1-ε)L ↓ ] / +L ↓

[0083] Among them, L λ is the radiation at the sensor or the top of atmosphere (TOA) radiation, ε is the surface emissivity, B(T s ) is the temperature T given by Planck's law s (T s =LST) blackbody radiation, L ↓ and L ↑ are the downwelling and upwelling atmospheric radiation, respectively, and τ is the total atmospheric transmittance between the land surface and the sensor. B(T S ) is expressed as:

[0084] B(T s )=[L λ -L ↑ -τ(1-ε)L ↓ ] / τε

[0085] According to Planck's law, we can get T S :

[0086] T s =K2 / ln(K1 / B(T s )+1)

[0087] Where K1=607.76Wm -2 μm -2 sr -1 and K2=1260.56Km -2 μm -2 sr -1 If obtained from TM, K1=774.89Wm -2 μm -2 sr -1 and K2=1321.08Km -2 μm -2 sr -1 For Thermal Infrared Sensor (TIRS) Band 10. Atmospheric parameter L ↓ , L↑, and τ were obtained from NASA's website (http: / / atmcorr.gsfc.nasa.gov / ).

[0088] Dryness (NDBSI): A combination of the Normalized Building Index and the Soil Index, reflecting the condition of buildings and bare soil. NDBSI dryness is usually expressed as the Bare Soil Index (SI). The index-based Building Land Index (IBI) precisely reflects the condition of building land;

[0089] The dryness is expressed as the average value of the normalized building index and the bare soil index, and the expression is:

[0090]

[0091] SI=[(R mir1 +R red )-(R nir +R blue )] / [(R mir1 +R red )+(R nir +R blue )]

[0092]

[0093] Among them, R red 、R blue 、R green 、R nir 、R mir1 They are the reflectances of the red band, blue band, green band, near-infrared band and short-wave infrared 1 band respectively.

[0094] In S2, the extracted ecological indicators are coupled using principal component analysis to extract the first principal component, and the remote sensing ecological index is calculated using the first principal component, specifically:

[0095] Normalizing the ecological indicators, calculating the first principal component using eigenanalysis based on GEE (GEE website: https: / / developers.google.com / earth-engine / array-s_eigen_analysis), and using the first principal component to calculate the initial remote sensing ecological indicators; normalizing the initial remote sensing ecological indicators to obtain the final remote sensing ecological index;

[0096] Among them, the expression of the initial remote sensing ecological index RSEI0 is:

[0097] RSEI0=1-PC1[f(NDVI,WET,NDBSI,LST)]

[0098] Among them, PC1 is the first principal component of the four indicators (NDVI, WET, NDSI and LST).

[0099] In S4, the Euclidean distance from each pixel in the cost grid map to the nearest road and the cost distance considering landscape resistance are calculated; wherein the cost distance considering landscape resistance is specifically: the average value of the remote sensing ecological index of different cumulative cost distances is calculated, and the calculation formula is:

[0100]

[0101] Among them, n k is the total number of pixels in the partition, Z k The cost distance interval (such as: 0-500m, 500-1km...), RSEI xy is the RSEI value within the pixel. In the calculation process, invalid values ​​within the partition (such as points with no value or outliers outside the theoretical range) need to be excluded.

[0102] In this embodiment, mutation point detection (binary detection method) is used to extract mutation points ( Figure 6 ).Depend on Figure 6 As can be seen, as the cost distance increases, the remote sensing ecological index first increases rapidly, then stabilizes after reaching the first mutation point (3000), continuing to stabilize until 10500. After a period of low values, it stabilizes again. Based on this, a cost distance of 3000 is selected as the road impact threshold in this case.

[0103] In summary, in order to improve the accuracy of road influence domain division and fully consider the different restrictive effects of different geographical spatial environments, a new method for road influence domain division based on the minimum cumulative resistance model and remote sensing ecological index is proposed ( Figure 1 ), providing strong support for road planning, management, and ecological protection. This method fully considers the objective distribution patterns of the geographic spatial environment. For example, areas with steep slopes and high forest coverage have good ecological environments, and road construction has a greater resistance to them; conversely, the resistance is smaller. Based on this, reclassification is performed based on remote sensing ecological indices to construct a resistance surface (where the ecological environment is good, the resistance value is large; where the ecological environment is poor, the resistance value is large). Based on the minimum cumulative resistance model, it expands outward from the road. When the ecological environment is poor in all areas it passes through, the cumulative resistance value is relatively small, and the range of influence reaching a certain fixed resistance value is large. Conversely, when the ecological environment is good in all areas it passes through, the cumulative resistance value is relatively large, and the range of influence reaching a certain fixed resistance value is small. Based on this, further combined with signal analysis to detect mutation points and extract cumulative resistance thresholds, it is possible to objectively divide the road influence domain. This road influence domain division method breaks through the limitation of the traditional fixed width of the road influence domain. Its width changes with changes in the objective geographic spatial environment, resulting in a more scientific result.

[0104] The above description of the disclosed embodiments is intended to enable one skilled in the art to implement or use the present invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention is not limited to the embodiments shown herein, but is intended to conform to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for dividing a road influence domain, characterized in that: The following steps are involved: S1. Obtain Landsat remote sensing image data of the study area, perform preprocessing, and extract ecological indicators; S2. Using principal component analysis to couple the extracted ecological indicators, extract the first principal component, and use the first principal component to calculate the remote sensing ecological index; S3. Divide the ecological environment quality into different levels according to the remote sensing ecological index, reclassify it, and construct a cost grid map; S4. Calculate the Euclidean distance from each pixel in the cost grid map to the nearest road and the cost distance considering landscape resistance, and analyze the variation pattern of remote sensing ecological index with cost distance; S5. Based on the mutation point detection method, the mutation points are extracted using the change pattern of remote sensing ecological index with cost distance, the road impact threshold is determined, and the road impact domain is delineated.

2. The road influence domain division method according to claim 1, characterized in that: In S1, the preprocessing includes: radiation correction, atmospheric correction, and geometric correction.

3. The road influence domain division method according to claim 1, characterized in that: In S1, the ecological indicators include: greenness, humidity, heat, and dryness.

4. The road influence domain division method according to claim 3, characterized in that: The greenness is the normalized vegetation index, and the calculation formula is: Among them, R red and R nir They represent the reflectance of the red band and near-infrared band of the Landsat image respectively.

5. The road influence domain division method according to claim 3, characterized in that: The humidity refers to the surface humidity, which is obtained by extracting the humidity component based on Landsat TM and OLI data. The calculation method is as follows: WET(TM)=0.0315R blue +0.2021R green +0.3102R red +0.1594R nir -0.6706R nir1 -0.6109R nir2 WET(OLI) <h2 style=";text-align:left;direction:ltr">=0.1511R<h2 style=";text-align:left;direction:ltr"> blue <h2 style=";text-align:left;direction:ltr"> +0.1973R<h2 style=";text-align:left;direction:ltr"> green <h2 style=";text-align:left;direction:ltr"> +0.3283R<h2 style=";text-align:left;direction:ltr"> red <h2 style=";text-align:left;direction:ltr"> +0.3407R<h2 style=";text-align:left;direction:ltr"> nir <h2 style=";text-align:left;direction:ltr"> -0.7117R<h2 style=";text-align:left;direction:ltr"> nir1 <h2 style=";text-align:left;direction:ltr"> -0.4599R<h2 style=";text-align:left;direction:ltr"> nir2 Where: R blue 、R green 、R red 、R nir 、R nir1 、R nir2 They are the reflectances of the blue, green, red, near-infrared, shortwave range 1, and shortwave range 2 bands, respectively.

6. The road influence domain division method according to claim 3, characterized in that: The heat refers to the thermal environment conditions, specifically: Land surface temperature was retrieved from the thermal band of Landsat imagery and processed using a radiative transfer equation algorithm to represent thermal indicators; The radiation transfer equation is expressed as follows: L λ =[εB(T s )+(1-ε)L ↓ ] / +L ↓ , Among them, L λ is the radiation at the sensor or the top of the atmosphere, ε is the surface emissivity, B(T s ) is the temperature T given by Planck's law s Blackbody radiation, L ↓ and L ↑ are the downwelling and upwelling atmospheric radiations, respectively, and τ is the total atmospheric transmittance between the land surface and the sensor; B(T S ) is: B(T s )=[L λ -L ↑ -τ(1-ε)L ↓ ] / te According to Planck's law, we can get T S : T s =K2 / ln(K1 / B(T s )+1) Among them, K1 and K2 are constants related to the wavelength of the selected band.

7. The road influence domain division method according to claim 3, characterized in that: The dryness is specifically the average value of the normalized building index and the bare soil index, expressed as: SI=[(R mir1 +R red )-(R nir +R blue )] / [(R mir1 +R red )+(R nir +R blue )] Among them, R red 、R blue 、R green 、R nir 、R mir1 are the reflectances of the red band, blue band, green band, near-infrared band and shortwave infrared 1 band, SI is the bare soil index, and IBI is the normalized building index.

8. The road influence domain division method according to claim 1, characterized in that: In S2, the extracted ecological indicators are coupled using principal component analysis to extract the first principal component, and the remote sensing ecological index is calculated using the first principal component, specifically: Normalizing the ecological indicators, calculating the first principal component based on GEE using feature analysis, and calculating the initial remote sensing ecological indicators using the first principal component; normalizing the initial remote sensing ecological indicators to obtain a final remote sensing ecological index; Among them, the expression of the initial remote sensing ecological index RSEI0 is: RSEI0=1-PC1[f(NDVI,WET,NDBSI,LST)] Among them, PC1 is the first principal component of ecological indicators.

9. The road influence domain division method according to claim 1, characterized in that: In S4, the Euclidean distance from each pixel in the cost grid map to the nearest road and the cost distance considering landscape resistance are calculated; wherein the cost distance considering landscape resistance is specifically: the average value of the remote sensing ecological index of different cumulative cost distances is calculated, and the calculation formula is: Among them, n k is the total number of pixels in the partition, Z k is the cost distance division interval, RSEI xy is the remote sensing ecological index value within the pixel.

Images