1. A four-dimensional remote sensing ecological index construction method based on a space geometric principle is characterized by comprising the following steps:

step S1: obtaining remote sensing data of a research area including blue light, green light, red light, near infrared, short wave infrared and thermal infrared bands, preprocessing the remote sensing data, further obtaining a two-dimensional space scatter diagram according to the red light band and the near infrared band, fitting a soil line according to the space scatter diagram and solving a soil line equation;

step S2: calculating a vertical vegetation index PVI;

step S3: calculating a vertical drought index (PDI), setting the limits of two ends of the mixed pixel to be bare soil and vegetation, and subtracting the value of the vegetation coverage PDI from the value of the PDI of the mixed pixel to obtain an improved vertical drought index (MPDI);

step S4: calculating the surface dryness NDSI;

step S5: calculating the surface temperature LST;

step S6: constructing a four-dimensional characteristic space by taking a vertical vegetation index PVI, an improved vertical drought index MPDI, a surface drying index NDSI and a surface temperature index LST as coordinate axes;

step S7: and calculating the remote sensing ecological index SGEI based on the vegetation, the surface drought degree, the surface drying degree and the surface temperature.

2. The method for constructing the four-dimensional remote sensing ecological index based on the spatial geometry principle according to claim 1, wherein the step S1 specifically comprises:

step S11: acquiring research area remote sensing data of blue light, green light, red light, near infrared, intermediate infrared and thermal infrared bands of a research area, preprocessing the remote sensing data, and extracting reflectivity data of the blue light, the green light, the red light, the near infrared, the short wave infrared and the thermal infrared bands;

step S12: constructing a two-dimensional space scatter diagram by taking the reflectivity data of the red light wave band as a horizontal coordinate and the reflectivity data of the near infrared wave band as a vertical coordinate;

step S13: based on the constructed two-dimensional scatter diagram distribution, constructing (X1, Y1) point sets through minimum longitudinal coordinate values corresponding to horizontal coordinates, and performing linear fitting through a least square method to obtain a soil line equation:

ρ_NIR＝M*ρ_Red+I (1)

in the formula, ρ_RedAnd ρ_NIRRespectively representing the reflectance values of the red and near infrared bands (range of value range [0,1]]) (ii) a M and I are the slope and intercept of the soil line equation.

3. The method for constructing the four-dimensional remote sensing ecological index based on the spatial geometry principle according to claim 1, wherein the step S2 specifically comprises: according to the soil line equation determined in step S1, the vertical vegetation index PVI is calculated as follows:

in the formula, ρ_NIRRepresenting the reflectivity, p, of the near infrared band_RedThe reflectivity of a red light wave band is represented; m and I represent the slope and intercept, respectively, of the soil line equation.

4. The method for constructing the four-dimensional remote sensing ecological index based on the spatial geometry principle according to claim 1, wherein the step S3 specifically comprises:

step S31: the vegetation coverage fv is calculated through a normalized vegetation index NDVI, the vegetation coverage fv is calculated through proportion in bare soil areas covered by all vegetation and not covered by plants in one image, the vegetation coverage is calculated by adopting a calculation method of semi-empirical vegetation coverage such as Baret, and the expression is as follows:

in the formula, NDVI_MAXAnd NDVI_MINNDVI under 100% vegetation coverage and NDVI corresponding to bare soil respectively;

step S32: according to the soil line equation determined in step S1, the vertical drought index PDI is calculated as follows:

in the formula, ρ_NIRRepresenting the reflectivity, p, of the near infrared band_RedThe reflectivity of a red light wave band is represented; m and I respectively represent the slope and intercept of a soil line equation;

step S33: using the improved vertical drought index MPDI, according to the soil line equation determined in step S1, the expression of MPDI is:

in the formula, ρ_S,RedRepresenting the reflectivity, rho, of bare earth in the red band_S,NIRThe reflectivity of the bare soil in a near infrared band is represented;

step S34: if the limits of two ends of the mixed pixel are bare soil and vegetation, the reflectivity of the mixed pixel is linearly combined with the bare soil and the vegetation, MPDI is deduced according to the linear combination, and the expression is as follows:

in the formula, ρ_RedAnd ρ_NIRRespectively the reflectivity of red light and near infrared wave band; rho_V,RedAnd ρ_V,NIRVegetation reflectivities of a red light band and a near infrared band respectively; m is the soil baseline slope.

5. The method for constructing the four-dimensional remote sensing ecological index based on the spatial geometry principle according to claim 1, wherein the step S4 specifically comprises:

the SI represents surface drying caused by bare soil (sand, bare soil and the like), and the calculation method comprises the following steps:

IBI is a building index which represents desiccation caused by construction sites in urban areas and is calculated by the following method:

where rho_SWIR1、ρ_blue、ρ_greenThe surface reflectivities of short wave infrared, blue light and green light wave bands are respectively, and other variables are the same as above;

the calculation method of the surface dryness is as follows:

6. the method for constructing the four-dimensional remote sensing ecological index based on the spatial geometry principle according to claim 1, wherein when the selected data is Landsat5 remote sensing data of 2000 years ago, the step S5 specifically comprises:

(1) image radiometric calibration, namely radiometric calibration is carried out on the Landsat5 thermal infrared band;

(2) calculating vegetation coverage, namely calculating vegetation coverage Fv by using a mixed pixel decomposition method, roughly dividing the land in the image into a water body, vegetation and a building, wherein a specific calculation formula is as follows:

F_v＝[(NDVI-NDVI_Soil)/(NDVI_Veg-NDVI_Soil)] (10)

wherein NDVI is the normalized vegetation index, NDVI_SoilNDVI value, NDVI, for a soil area without vegetation cover_VegThe NDVI value is the NDVI value of the area completely covered by vegetation;

(3) surface emissivity epsilon calculation

The emissivity of the water body pixel is assigned to be 0.995, and the emissivity of the natural surface and town pixels is estimated according to the following formula:

εsurface＝0.9625+0.0614F_v-0.0461F_v ² (11)

εbuilding＝0.9589+0.086F_v-0.0671F_v ² (12)

in the formula, epsilon surface and epsilon building respectively represent the specific radiance of a natural surface pixel and an urban pixel;

(4) spectral radiance L of pixels at the sensor_λComputing

L_λ＝gain·Q_λ+bias (13)

In the formula: l is_λThe spectral radiation value of a pixel of a thermal infrared band at a sensor is W/(m2.ster. mu.m); gain is the gain value of the corresponding wave band; bias is the bias value of the corresponding wave band; q_λCorresponding DN value to the pixel;

(5) calculating the brightness temperature value of the sensor

T_λ＝K2/ln(K1/L_λ+1) (14)

In the formula, T_λIs the value of the light temperature at the sensor, with the unit of K;

LST＝T_λ/[1+(λ·T_λ/ρ)lnε] (15)

in the formula, lambda is the central wavelength or effective wavelength of a thermal infrared band, and the unit is mum; k1 and K2 are scaling parameters, respectively.

7. The method for constructing the ecological index of four-dimensional remote sensing based on the spatial geometry principle according to claim 1, wherein when the selected data is Landsat5 remote sensing data after 2000 years, the step S5 specifically comprises:

(1) image radiometric calibration, namely radiometric calibration is carried out on the Landsat5 thermal infrared band;

(2) calculating vegetation coverage, namely calculating vegetation coverage Fv by using a mixed pixel decomposition method, roughly dividing the land in the image into a water body, vegetation and a building, wherein a specific calculation formula is as follows:

F_v＝[(NDVI-NDVI_Soil)/(NDVI_Veg-NDVI_Soil)] (10)

wherein NDVI is the normalized vegetation index, NDVI_SoilNDVI value for an area without vegetation coverage, NDVI_VegNDVI value for the area that is completely covered by the implant;

(3) surface emissivity epsilon calculation

The emissivity of the water body pixel is assigned to be 0.995, and the emissivity of the natural surface and town pixels is estimated according to the following formula:

εsurface＝0.9625+0.0614F_v-0.0461F_v ² (11)

εbuilding＝0.9589+0.086F_v-0.0671F_v ² (12)

in the formula, epsilon surface and epsilon building respectively represent the specific radiance of a natural surface pixel and an urban pixel;

(4) calculating the radiation brightness value of the black body at the same temperature

Thermal infrared radiation brightness value L received by satellite sensor_λThe method comprises three parts: the atmosphere radiates brightness upwards, and the real radiation brightness of the ground reaches the energy of the satellite sensor after passing through the atmosphere; the atmosphere radiates energy downward that is reflected after reaching the ground. It calculatesThe method comprises the following steps:

L_λ＝[ε·B(T_S)+(1-ε)L↓]·τ+L↑ (13)

here, ε is the ground emissivity, T_SAs true surface temperature, B (T)_S) Is black body at T_Sτ is the transmittance of the atmosphere in the thermal infrared band. The radiant brightness B (T) of the black body with the temperature T in the thermal infrared band_S) Comprises the following steps:

B(T_S)＝[L_λ-L↑-τ·(1-ε)L↓]/τ·ε (14)

(5) inversion of earth surface temperature

At a pickup temperature T_SAfter the radiation brightness of the black body in the thermal infrared band, the real surface temperature T is obtained according to the inverse function of the Planck formula_S：

T_S＝K₂/ln(K₁/B(T_S)+1) (15)

In the formula, B (T)_S) Is the black body radiant brightness.

8. The method for constructing the four-dimensional remote sensing ecological index based on the spatial geometry principle according to claim 1, wherein when the selected data is Landsat8 remote sensing data, the step S5 specifically comprises:

(1) image radiometric calibration, namely radiometric calibration is carried out on the Landsat8 thermal infrared band;

(2) calculating vegetation coverage, namely calculating vegetation coverage Fv by using a mixed pixel decomposition method, roughly dividing the land in the image into a water body, vegetation and a building, wherein a specific calculation formula is as follows:

F_v＝[(NDVI-NDVI_Soil)/(NDVI_Veg-NDVI_Soil)] (10)

wherein NDVI is the normalized vegetation index, NDVI_SoilNDVI value for an area without vegetation coverage, NDVI_VegNDVI value for the area that is completely covered by the implant;

(3) and calculating the earth surface emissivity epsilon, wherein the emissivity of the water body pixel is assigned to be 0.995, and the emissivity of the natural surface and the town pixel is estimated according to the following formula:

εsurface＝0.9625+0.0614F_v-0.0461F_v ² (11)

εbuilding＝0.9589+0.086F_v-0.0671F_v ² (12)

in the formula, epsilon surface and epsilon building respectively represent the specific radiance of a natural surface pixel and an urban pixel;

(4) calculating the radiation brightness value of the black body at the same temperature

Thermal infrared radiation brightness value L received by satellite sensor_λThe method comprises three parts: the atmosphere radiates brightness upwards, and the real radiation brightness of the ground reaches the energy of the satellite sensor after passing through the atmosphere; the atmosphere radiates energy reflected after reaching the ground downwards; the calculation method comprises the following steps:

L_λ＝[ε·B(T_S)+(1-ε)L↓]·τ+L↑ (13)

here, ε is the ground emissivity, T_SAs true surface temperature, B (T)_S) Is black body at T_Sτ is the transmittance of the atmosphere in the thermal infrared band. The radiant brightness B (T) of the black body with the temperature T in the thermal infrared band_S) Comprises the following steps:

B(T_S)＝[L_λ-L↑-τ·(1-ε)L↓]/τ·ε (14)

inputting imaging time and central longitude and latitude in NASA official network (http:// atmcorr. gsfc. NASA. gov /), the parameters required in the above formula are provided;

(5) inversion of earth surface temperature

At a pickup temperature T_SAfter the radiation brightness of the black body in the thermal infrared band, the real surface temperature T is obtained according to the inverse function of the Planck formula_S：

T_S＝K₂/ln(K₁/B(T_S)+1) (15)

In the formula, B (T)_S) Is the black body radiant brightness; for Landsat8 TIRS, K₁＝774.89W/(m²*μm*sr)，K₂＝1321.08K。

9. The method for constructing the four-dimensional remote sensing ecological index based on the spatial geometry principle according to claim 1, wherein the step S7 specifically comprises:

step S71: the values for each axis are normalized to [0,1/2], and then the SGEI is calculated to range from [0,1 ]. Wherein, PVI is normalized in forward direction, LST, MPDI, NDSI are normalized in backward direction

Forward normalization formula:

inverse normalization formula:

in the formula, Normalizedp_iIs a normalized value;

step S72: in the four-dimensional feature space created in step 6, an expression of an index SGEI for representing ecological quality is constructed as follows:

the SGEI is a four-dimensional remote sensing ecological index based on the greenness, the soil dryness, the earth surface temperature and the earth surface dryness, the numerical range is [0,1], the closer the SGEI is to 1, the better the ecological condition is, and on the contrary, the worse the ecological condition is.

Background

Regional ecological evaluation indexes have been gradually emphasized by people in recent years, and remote sensing has been widely applied as an important means for quantitatively evaluating regional ecological environments. Remote sensing obtains data repeatedly in a certain period, simultaneously can synchronously observe in a large area, truly reflects regional ecological environment change, and rapidly obtains target object attribute characteristics, thereby carrying out ecological evaluation.

The most common evaluation of ecological conditions by remote sensing data is to evaluate one aspect of the ecology from a single perspective by using remote sensing data, but the ecological conditions cannot be evaluated as a whole. Some researchers integrate a plurality of indexes through a principal component transformation or factor analysis method to construct a comprehensive remote sensing index reflecting the ecological condition, such as RSEI (Xuequu, 2013), EEM (Zhang,2016) and the like, but the physical meanings of the indexes are unclear; VMTEI (zhongziyan et al, 2020) utilizes the principles of space geometry, but only takes into account the 3 indices vegetation cover, surface moisture and surface temperature. Based on multi-band remote sensing data, the invention couples vertical vegetation index PVI representing greenness, improved vertical drought index MPDI (PVI and MPDI are vertical in space) representing dryness, surface dryness NDSI and surface temperature LST, introduces a space geometry principle to construct space geometry ecological index SGEI, improves the physical significance of the ecological index SGEI, and can comprehensively reflect the space distribution of the ecological environment condition of a region.

Disclosure of Invention

In view of the above, the present invention provides a method for constructing a four-dimensional remote sensing ecological index based on a spatial geometry principle, so as to overcome the defects in the prior art.

In order to achieve the purpose, the invention adopts the following technical scheme:

a four-dimensional remote sensing ecological index construction method based on a space geometric principle comprises the following steps:

step S1: obtaining remote sensing data of a research area including blue light, green light, red light, near infrared, short wave infrared and thermal infrared bands, preprocessing the remote sensing data, further obtaining a two-dimensional space scatter diagram according to the red light band and the near infrared band, fitting a soil line according to the scatter diagram and solving a soil line equation;

step S2: calculating a vertical vegetation index PVI;

step S3: calculating a vertical drought index (PDI), setting the limits of two ends of the mixed pixel to be bare soil and vegetation, and subtracting the value of the vegetation coverage PDI from the value of the PDI of the mixed pixel to obtain an improved vertical drought index (MPDI);

step S4: calculating the surface dryness NDSI;

step S5: calculating the surface temperature LST;

step S6: constructing a four-dimensional characteristic space by taking a vertical vegetation index PVI, an improved vertical drought index MPDI, a surface drying index NDSI and a surface temperature index LST as coordinate axes;

step S7: and calculating the remote sensing ecological index SGEI based on the vegetation, the surface drought degree, the surface drying degree and the surface temperature.

Further, the step S1 is specifically:

step S11: acquiring research area remote sensing data of blue light, green light, red light, near infrared, short wave infrared and thermal infrared bands of a research area, carrying out pretreatment such as geometric correction, radiometric calibration, atmospheric correction and the like on the remote sensing data, and extracting reflectivity data of the blue light, the green light, the red light, the near infrared, the short wave infrared and the thermal infrared bands;

step S12: constructing a two-dimensional space scatter diagram by taking the reflectivity data of the red light wave band as a horizontal coordinate and the reflectivity data of the near infrared wave band as a vertical coordinate;

step S13: based on the distribution of the constructed space two-dimensional scatter diagram, constructing (X1, Y1) point sets through minimum longitudinal coordinate values corresponding to horizontal coordinates, and fitting the soil line through a least square method to obtain an equation of the soil line:

ρ_NIR＝M*ρ_Red+I (1)

in the formula, ρ_RedAnd ρ_NIRRespectively representing the reflectance values of the red and near infrared bands (range of value range [0,1]]) (ii) a M and I refer to the slope and intercept of the soil line equation。

Further, the step S2 is specifically: according to the soil line equation determined in step S1, the vertical vegetation index PVI is calculated as follows:

in the formula, ρ_NIRRepresenting the reflectivity, p, of the near infrared band_RedThe reflectivity of a red light wave band is represented; m and I represent the slope and intercept, respectively, of the soil line equation.

Further, the step S3 is specifically:

step S31: vegetation coverage index f_vThe vegetation coverage is calculated by a normalized vegetation index NDVI, the vegetation coverage is calculated by proportion of bare soil areas covered by all vegetation in the image and bare soil areas not covered by plants, a semi-empirical vegetation coverage calculation method such as Baret is selected as the vegetation coverage, and the expression is as follows:

in the formula, NDVI_MAXAnd NDVI_MINNDVI under 100% vegetation coverage and NDVI corresponding to bare soil respectively;

step S32: according to the soil line equation determined in step S1, the vertical drought index PDI is calculated as follows:

in the formula, ρ_NIRRepresenting the reflectivity, p, of the near infrared band_RedThe reflectivity of a red light wave band is represented; m and I respectively represent the slope and intercept of a soil line equation;

step S33: using the improved vertical drought index MPDI, according to the soil line equation determined in step S1, the expression of MPDI is:

in the formula, ρ_S,RedRepresenting the reflectivity, rho, of bare earth in the red band_S,NIRThe reflectivity of the bare soil in a near infrared band is represented;

step S34: if the limits of two ends of the mixed pixel are bare soil and vegetation, the reflectivity of the mixed pixel is linearly combined with the bare soil and the vegetation, MPDI is deduced according to the linear combination, and the expression is as follows:

in the formula, ρ_RedAnd ρ_NIRRespectively the reflectivity of red light and near infrared wave bands which are corrected by atmosphere in the remote sensing image; rho_V,RedAnd ρ_V,NIRVegetation reflectivities of a red light band and a near infrared band respectively; m is the soil baseline slope.

Further, the step S4 is specifically:

the SI represents surface drying caused by bare soil (sand, bare soil and the like), and the calculation method comprises the following steps:

IBI is a building index which represents desiccation caused by construction sites in urban areas and is calculated by the following method:

where rho_SWIR1、ρ_blue、ρ_greenThe surface reflectivities of the short wave infrared, blue and green wave bands are respectively, and other variables are the same as above.

The calculation method of the surface dryness is as follows:

further, when the selected data is remote sensing data Landsat5 before 2000, the step S5 specifically includes:

(1) image radiometric calibration, namely radiometric calibration is carried out on the Landsat5 thermal infrared band;

(2) calculating vegetation coverage, namely calculating vegetation coverage Fv by using a mixed pixel decomposition method, roughly dividing the land in the image into a water body, vegetation and a building, wherein a specific calculation formula is as follows:

F_v＝[(NDVI-NDVI_Soil)/(NDVI_Veg-NDVI_Soil)] (10)

wherein NDVI is the normalized vegetation index, NDVI_SoilNDVI value, NDVI, for a soil area without vegetation cover_VegThe NDVI value is the NDVI value of the area completely covered by vegetation;

(3) surface emissivity epsilon calculation

The emissivity of the water body pixel is assigned to be 0.995, and the emissivity of the natural surface and town pixels is estimated according to the following formula:

εsurface＝0.9625+0.0614F_v-0.0461F_v ² (11)

εbuilding＝0.9589+0.086F_v-0.0671F_v ² (12)

in the formula, epsilon surface and epsilon building respectively represent the specific radiance of a natural surface pixel and an urban pixel;

(4) spectral radiance L of pixels at the sensor_λComputing

L_λ＝gain·Q_λ+bias (13)

In the formula: l is_λThe spectral radiation value of a pixel of a thermal infrared band at a sensor is W/(m2.ster. mu.m); gain is the gain value of the corresponding wave band; bias is the bias value of the corresponding wave band; q_λCorresponding DN value to the pixel;

(5) calculating the brightness temperature value of the sensor

T_λ＝K2/ln(K1/L_λ+1) (14)

In the formula, T_λIs the value of the light temperature at the sensor, with the unit of K;

LST＝T_λ/[1+(λ·T_λ/ρ)lnε] (15)

in the formula, lambda is the central wavelength or effective wavelength of a thermal infrared band, and the unit is mum; k1 and K2 are scaling parameters, respectively.

Further, when the selected data is remote sensing data Landsat5 after 2000 years, the step S5 specifically includes:

(1) image radiometric calibration, namely radiometric calibration is carried out on the Landsat5 thermal infrared band;

(2) calculating vegetation coverage, calculating vegetation coverage F by using mixed pixel decomposition method_vThe land types in the image are roughly divided into water bodies, vegetation and buildings, and the specific calculation formula is as follows:

F_v＝[(NDVI-NDVI_Soil)/(NDVI_Veg-NDVI_Soil)] (10)

wherein NDVI is the normalized vegetation index, NDVI_SoilNDVI value for an area without vegetation coverage, NDVI_VegNDVI value for the area that is completely covered by the implant;

(3) surface emissivity epsilon calculation

The emissivity of the water body pixel is assigned to be 0.995, and the emissivity of the natural surface and town pixels is estimated according to the following formula:

εsurface＝0.9625+0.0614F_v-0.0461F_v ² (11)

εbuilding＝0.9589+0.086F_v-0.0671F_v ² (12)

in the formula, epsilon surface and epsilon building respectively represent the specific radiance of a natural surface pixel and an urban pixel;

(4) calculating the radiation brightness value of the black body at the same temperature

Thermal infrared radiation brightness value L received by satellite sensor_λThe method comprises three parts: the brightness of the radiation from the atmosphere upwards, the true brightness of the radiation from the ground passing through the atmosphereEnergy reaching the satellite sensor after the layer; the atmosphere radiates energy downward that is reflected after reaching the ground. The calculation method comprises the following steps:

L_λ＝[ε·B(T_S)+(1-ε)L↓]·τ+L↑ (13)

here, ε is the ground emissivity, T_SAs true surface temperature, B (T)_S) Is black body at T_Sτ is the transmittance of the atmosphere in the thermal infrared band. The radiant brightness B (T) of the black body with the temperature T in the thermal infrared band_S) Comprises the following steps:

B(T_S)＝[L_λ-L↑-τ·(1-ε)L↓]/τ·ε (14)

(5) inversion of earth surface temperature

At a pickup temperature T_SAfter the radiation brightness of the black body in the thermal infrared band, the real surface temperature T is obtained according to the inverse function of the Planck formula_S：

T_S＝K2/ln(K1/B(T_S)+1) (15)

In the formula, B (T)_S) Is the black body radiant brightness.

Further, when the selected data is Landsat8 remote sensing data, the step S5 specifically includes:

(1) image radiometric calibration, namely radiometric calibration is carried out on the Landsat8 thermal infrared band;

(2) calculating vegetation coverage, calculating vegetation coverage F by using mixed pixel decomposition method_vThe land types in the image are roughly divided into water bodies, vegetation and buildings, and the specific calculation formula is as follows:

F_v＝[(NDVI-NDVISoil)/(NDVIVeg-NDVISoil)] (10)

in the formula, NDVI is a normalized vegetation index, NDVISUI is an NDVI value of an area without vegetation coverage, and NDVIVeg is an NDVI value of an area which is completely planted and covered;

(3) and calculating the earth surface emissivity epsilon, wherein the emissivity of the water body pixel is assigned to be 0.995, and the emissivity of the natural surface and the town pixel is estimated according to the following formula:

εsurface＝0.9625+0.0614F_v-0.0461F_v ² (11)

εbuilding＝0.9589+0.086F_v-0.0671F_v ² (12)

in the formula, epsilon surface and epsilon building respectively represent the specific radiance of a natural surface pixel and an urban pixel;

(4) calculating the radiation brightness value of the black body at the same temperature

Thermal infrared radiation brightness value L received by satellite sensor_λThe method comprises three parts: the atmosphere radiates brightness upwards, and the real radiation brightness of the ground reaches the energy of the satellite sensor after passing through the atmosphere; the atmosphere radiates energy reflected after reaching the ground downwards; the calculation method comprises the following steps:

L_λ＝[ε·B(T_S)+(1-ε)L↓]·τ+L↑ (13)

here, ε is the ground emissivity, T_SAs true surface temperature, B (T)_S) Is black body at T_Sτ is the transmittance of the atmosphere in the thermal infrared band. The radiant brightness B (T) of the black body with the temperature T in the thermal infrared band_S) Comprises the following steps:

B(T_S)＝[L_λ-L↑-τ·(1-ε)L↓]/τ·ε (14)

inputting imaging time and central longitude and latitude in NASA official network (http:// atmcorr. gsfc. NASA. gov /), the parameters required in the above formula are provided;

(5) inversion of earth surface temperature

At a pickup temperature T_SAfter the radiation brightness of the black body in the thermal infrared band, the real surface temperature T is obtained according to the inverse function of the Planck formula_S：

T_S＝K₂/ln(K₁/B(T_S)+1) (15)

Wherein B (TS) is the black body radiation brightness; for Landsat8 TIRS, K₁＝774.89W/(m²*μm*sr)，K₂＝1321.08K。

10. The method for constructing the four-dimensional remote sensing ecological index based on the spatial geometry principle according to claim 1, wherein the step S7 specifically comprises:

step S71: the values for each axis are normalized to [0,1/2], and then the SGEI is calculated to range from [0,1 ]. Wherein, PVI is normalized in forward direction, LST, MPDI, NDSI are normalized in backward direction

Forward normalization formula:

inverse normalization formula:

in the formula, Normalizedp_iIs a normalized value;

step S72: in the three-dimensional feature space created in step 6, an expression of an index SGEI for representing ecological quality is constructed as follows:

the SGEI is a four-dimensional remote sensing ecological index based on the greenness, the soil dryness, the earth surface temperature and the earth surface dryness, the numerical range is [0,1], the closer the SGEI is to 1, the better the ecological condition is, and on the contrary, the worse the ecological condition is.

Compared with the prior art, the invention has the following beneficial effects:

the invention not only considers a plurality of indexes closely related to ecology, but also improves the physical significance of SGEI by introducing the space geometry principle. In addition, the SGEI can be used as a quantitative index to describe the spatial ecological quality, and can also be used for visualizing the ecological quality space and carrying out space-time analysis.

Drawings

FIG. 1 is a flow chart of a construction method of the present invention;

FIG. 2 is a diagram of a coordinate system during construction of a spatial geometric ecological index SGEI according to an embodiment of the present invention;

FIG. 3 is a plot of fitted soil line equations in one embodiment of the present invention;

FIG. 4 is a spatial distribution diagram of PVI, MPDI, NDSI, and LST in accordance with an embodiment of the present invention;

FIG. 5 is a spatial distribution diagram of the SGEI in an embodiment of the present invention.

Detailed Description

The invention is further explained below with reference to the drawings and the embodiments.

Referring to fig. 1, the invention provides a method for constructing a four-dimensional remote sensing ecological index (SGEI) based on a space geometric principle, comprising the following steps:

step S1: obtaining remote sensing data of a research area Landsat including blue light, green light, red light, near infrared, short wave infrared and thermal infrared bands, preprocessing the remote sensing data, fitting a soil line according to the red light band and a near infrared two-dimensional space scatter diagram, and solving a soil line equation:

the specific steps of step S1 are as follows:

step S11: and acquiring research area remote sensing data of blue light, green light, red light, near infrared, short wave infrared and thermal infrared bands of the research area Landsat, carrying out pretreatment such as geometric correction, radiometric calibration, atmospheric correction and the like on the remote sensing data, and extracting reflectivity data of the blue light, green light, red light, near infrared, short wave infrared and thermal infrared bands.

Step S12: and constructing a two-dimensional space scatter diagram according to the data of the reflectivity of the red light wave band and the reflectivity of the near infrared wave band.

Step S13: according to the distribution of the constructed spatial two-dimensional scatter diagram, fitting a soil line equation:

ρ_NIR＝1.1373*ρ_Red-0.0825 (1)

where rho_NIRAnd ρ_RedRespectively representing the reflectance values of the near infrared band and the red light band (both 0 to 1)

Step S2: and selecting the vertical distance from any point to the soil line according to the slope and the intercept of the soil line to represent the covering condition of the plants in the area, wherein the distances from the reflectivities of different types of the plants to the soil line are different, and calculating the vertical vegetation index PVI.

The specific steps of step S2 are as follows:

according to the soil line equation determined in step S13, the vertical vegetation index PVI is calculated as follows:

in the formula, ρ_NIRRepresenting the reflectivity, p, of the near infrared band_RedThe reflectivity of a red light wave band is represented; m and I represent the slope and intercept, respectively, of the soil line equation. According to equation (1), in this example, M is 1.1373 and I is 0.0825, the specific expression of PVI is:

step S3: and subtracting the value of vegetation coverage PDI from the value of PDI of the mixed pixel to obtain MPDI.

The specific steps of step S3 are as follows:

calculating the vegetation coverage, wherein the expression is as follows:

in the formula, NDVI_MAXAnd NDVI_MINNDVI under 100% vegetation coverage and NDVI corresponding to bare soil, respectively.

According to the soil line equation determined in step 13, the MPDI expression is:

in the formula, ρ_S,RedRepresenting the reflectivity, rho, of bare earth in the red band_S,NIRThe reflectivity of the bare soil in the near infrared band is shown.

The improvement part mainly considers the influence of implantation coverage, and the MPDI expression is as follows:

in the formula, ρ_RedAnd ρ_NIRRespectively the reflectivity rho of red light and near infrared wave band corrected by atmosphere in the remote sensing image_V,RedAnd ρ_V,NIRThe vegetation reflectivities of a red light wave band and a near infrared wave band respectively. Obtaining the numerical expression:

will rho_V,RedAnd ρ_V,NIR0.05 and 0.5 were used, respectively. M is the soil baseline slope, and if the soil baseline slope is calculated for years, the average slope of the soil baseline can be considered.

Step S4: and calculating the surface drying degree NDSI, and synthesizing by using An impervious surface building Index IBI (An Index-base brick-up Index) and a bare soil Index SI (soil Index).

Step S41: the SI represents surface drying caused by bare soil (sand, bare soil and the like), and the calculation method comprises the following steps:

step S42: IBI represents the building index and the desiccation caused by the construction land of an urban area, and the calculation method comprises the following steps:

where rho_SWIR1、ρ_blue、ρ_greenThe surface reflectivities of the short wave infrared, blue and green wave bands are respectively, and other variables are the same as above.

Step S43: the method for calculating the surface dryness comprises the following steps:

step S5: the surface temperature LST is calculated.

In this example, Landsat8 TIRS was used to calculate LST.

Step S6: and constructing a four-dimensional feature space by using a vertical vegetation index PVI, an improved vertical drought index MPDI, a surface drying degree index NDSI and a surface temperature index LST.

Step S7: and calculating the remote sensing ecological index SGEI based on the vegetation, the surface drought degree, the surface drying degree and the surface temperature.

In order to make the index more reasonable and facilitate comparison, the value of each axis is normalized to [0,1/2], and then the SGEI is calculated to make the value range of the SGEI be [0,1 ]. The PVI is normalized in the forward direction, and the LST, MPDI, NDSI are normalized in the backward direction.

Forward normalization formula:

inverse normalization formula:

in the formula, Normalizedp_iIs a normalized value.

In the four-dimensional feature space created in step 7, an expression of an index SGEI for representing ecological quality is constructed as follows:

the SGEI is a four-dimensional remote sensing ecological index based on the greenness, the soil dryness, the earth surface temperature and the earth surface dryness, the numerical range is [0,1], the closer the SGEI is to 1, the better the ecological condition is, and on the contrary, the worse the ecological condition is.

In the embodiment, the spatial geometry remote sensing ecological index SGEI is calculated by utilizing Landsat8 data of 2013 month 10 in inland area of Fuzhou city, Fujian province. FIG. 4 is a diagram of the spatial distribution of PVI, MPDI, NDSI and LST, and FIG. 5 is a diagram of the spatial distribution of SGEI in an embodiment of the present invention.

As can be seen from the figure, the vegetation coverage of the inland areas of Fuzhou is high except for the main urban area; the earth surface dryness is higher in main urban areas of two banks of the Minjiang and higher in southeast areas of the Fuzhou; the surface dryness and the vegetation index show a high negative correlation, and the surface temperature and the surface dryness show a relatively obvious positive correlation; the distribution of the space geometric ecological index SGEI accords with the space distribution of a vertical vegetation index, an improved vertical drought index, a ground surface drying degree and a ground surface temperature, and overall, the higher the greenness, the lower the ground surface drought degree, the lower the ground surface drying degree and the lower the ground surface temperature, the better the ecological condition represented by the area is, and otherwise, the worse the ecological condition is.

The above description is only a preferred embodiment of the present invention, and all equivalent changes and modifications made in accordance with the claims of the present invention should be covered by the present invention.