CN115203624B - Comprehensive evaluation method for surface environment at any moment based on time sequence remote sensing - Google Patents
Comprehensive evaluation method for surface environment at any moment based on time sequence remote sensing Download PDFInfo
- Publication number
- CN115203624B CN115203624B CN202210810131.0A CN202210810131A CN115203624B CN 115203624 B CN115203624 B CN 115203624B CN 202210810131 A CN202210810131 A CN 202210810131A CN 115203624 B CN115203624 B CN 115203624B
- Authority
- CN
- China
- Prior art keywords
- remote sensing
- time sequence
- pixel
- index
- model
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000011156 evaluation Methods 0.000 title claims abstract description 33
- 230000008859 change Effects 0.000 claims abstract description 43
- 238000000034 method Methods 0.000 claims abstract description 42
- 238000011160 research Methods 0.000 claims abstract description 17
- 238000003384 imaging method Methods 0.000 claims abstract description 7
- 238000007781 pre-processing Methods 0.000 claims abstract description 5
- 238000012512 characterization method Methods 0.000 claims abstract description 4
- 238000002310 reflectometry Methods 0.000 claims description 17
- 239000002689 soil Substances 0.000 claims description 17
- 230000014509 gene expression Effects 0.000 claims description 14
- 238000012937 correction Methods 0.000 claims description 12
- 230000036541 health Effects 0.000 claims description 10
- 230000003287 optical effect Effects 0.000 claims description 10
- 230000008569 process Effects 0.000 claims description 10
- 238000010606 normalization Methods 0.000 claims description 8
- 238000000926 separation method Methods 0.000 claims description 8
- 238000009826 distribution Methods 0.000 claims description 7
- IJGRMHOSHXDMSA-UHFFFAOYSA-N Atomic nitrogen Chemical compound N#N IJGRMHOSHXDMSA-UHFFFAOYSA-N 0.000 claims description 6
- 101000812677 Homo sapiens Nucleotide pyrophosphatase Proteins 0.000 claims description 6
- 101000841763 Homo sapiens Ubiquinol-cytochrome-c reductase complex assembly factor 2 Proteins 0.000 claims description 6
- 102100039306 Nucleotide pyrophosphatase Human genes 0.000 claims description 6
- 102100029513 Ubiquinol-cytochrome-c reductase complex assembly factor 2 Human genes 0.000 claims description 6
- 230000005855 radiation Effects 0.000 claims description 6
- 230000009466 transformation Effects 0.000 claims description 5
- 238000013441 quality evaluation Methods 0.000 claims description 4
- 238000012886 linear function Methods 0.000 claims description 3
- 229910052757 nitrogen Inorganic materials 0.000 claims description 3
- 230000001131 transforming effect Effects 0.000 claims description 3
- 238000012163 sequencing technique Methods 0.000 claims description 2
- 230000000737 periodic effect Effects 0.000 abstract description 6
- 238000001514 detection method Methods 0.000 description 6
- 238000011161 development Methods 0.000 description 3
- 230000007613 environmental effect Effects 0.000 description 3
- 241000282414 Homo sapiens Species 0.000 description 2
- 241001272720 Medialuna californiensis Species 0.000 description 2
- 238000004364 calculation method Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000003915 air pollution Methods 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000012217 deletion Methods 0.000 description 1
- 230000037430 deletion Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000005184 irreversible process Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000001303 quality assessment method Methods 0.000 description 1
- 238000011895 specific detection Methods 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000003786 synthesis reaction Methods 0.000 description 1
- 230000002194 synthesizing effect Effects 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
- 238000003911 water pollution Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
- G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
- G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
- G01N2021/1793—Remote sensing
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Data Mining & Analysis (AREA)
- Computational Mathematics (AREA)
- Mathematical Optimization (AREA)
- Immunology (AREA)
- Pathology (AREA)
- Biochemistry (AREA)
- Algebra (AREA)
- Analytical Chemistry (AREA)
- Chemical & Material Sciences (AREA)
- Mathematical Analysis (AREA)
- General Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Pure & Applied Mathematics (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Image Processing (AREA)
- Investigating Or Analysing Materials By Optical Means (AREA)
Abstract
The invention discloses a comprehensive evaluation method of an earth surface environment at any moment based on time sequence remote sensing, which comprises the following specific steps: (1) Preprocessing the real shooting remote sensing image of the selected area in the research period, and arranging all preprocessed images to form an image time sequence; (2) Respectively establishing an initial time sequence model for each pixel in the image; (3) Detecting the change of the initial time sequence model of the pixel; (4) obtaining a predicted remote sensing image at any moment; (5) Taking the five remote sensing indexes as comprehensive characterization of the surface environment of the selected area; (6) denoising and normalizing the remote sensing index; (7) obtaining a comprehensive index of the surface environment; the method has the advantages that the periodic evaluation work of the surface environment is possible, the problem that the periodic evaluation is difficult to develop due to the lack of images is effectively avoided, and the problem that the evaluation results are difficult to compare due to the difference of imaging time of the images is solved.
Description
Technical Field
The invention relates to a remote sensing evaluation method of a ground surface environment, in particular to a comprehensive evaluation method of the ground surface environment at any moment based on time sequence remote sensing.
Background
The surface environment is highly relevant to the quality of life of human beings and the sustainable development of the current society. For half a century, the rapid urban process in China has led to a series of environmental problems such as forest cutting, loss of diversity, air and water pollution, land desertification, urban heat island effect, etc. In view of the fact that urbanization is an irreversible process in the development process of human civilization, the stress of the urbanization on the ecological environment is difficult to relieve in a short period of time. Therefore, the method is particularly important for developing socioeconomic development planning coordinated with regional environment capacity in order to explore the change condition of the surface environment in the urban process in time and objectively, and becomes a core issue of urban ecology research in recent years.
The remote sensing technology has the advantages of wide coverage range, various observation scales, rich spectral information, convenient data acquisition and the like, and is becoming an important means for representing the surface environment at home and abroad increasingly. Early description of the surface environment focused mainly on the calculation of a single index based on remote sensing, with particular emphasis on a particular aspect of the surface environment. Typical examples are Normalized Differential Vegetation Index (NDVI) and surface temperature (LST) are often used to describe vegetation cover and urban heat island conditions. In most cases, the surface environment is very complex to change, and it is difficult to represent by a single remote sensing index. Therefore, some scholars propose model schemes that integrate a variety of remote sensing-based indices, such as Forest Disturbance Index (FDI), global disturbance index (MGDI), and niche modeling (ENM), etc. While these composite indices can more fully reflect the surface environmental conditions, a number of challenges remain, such as weight selection. In practice, it is very difficult to determine the magnitude of the effect of each index on the change of the surface environment and assign weights proportionally, the weight assignment is inevitably different due to subjective experience differences, and the different weights directly affect the surface environment evaluation result. In addition, since the comprehensive index model scheme is very sensitive to the weather of the year due to a series of indexes (such as greenness and heat), most researches adopt to select data with better quality and similar seasons in the annual images or conduct a small-scale study, but available data deletion exists in coastal urban areas with more cloud, so that the application accuracy and comparability of the method and the area range and time sequence length of the study are severely limited.
Disclosure of Invention
The technical problem to be solved by the invention is to provide the comprehensive evaluation method of the earth surface environment at any moment based on time sequence remote sensing, which has high degree of automation and strong robustness, so that the earth surface environment evaluation at any moment is possible, the problem that the periodic evaluation is difficult to develop due to image missing is effectively avoided, and the problem that the evaluation result is difficult to compare due to the difference of imaging time is solved.
The technical scheme adopted for solving the technical problems is as follows: a comprehensive evaluation method of an earth surface environment at any moment based on time sequence remote sensing comprises the following specific steps:
(1) Preprocessing the real-shot remote sensing images of the selected area in the research period, and orderly arranging all preprocessed images according to the imaging time sequence to form an image time sequence;
(2) Respectively establishing an initial time sequence model for each pixel in the image on all wave bands according to the time sequence;
(3) Detecting the change of the initial time sequence model of the pixel, and if the change does not occur, maintaining the initial time sequence model; if the change occurs, taking the change occurrence time as a breakpoint of the time sequence model, segmenting the initial time sequence model of the pixel to obtain a segmented time sequence model so as to reflect the change of the earth surface reflectivity of the corresponding pixel;
(4) Obtaining a predicted remote sensing image at any moment on the basis of a time sequence model of the pixels;
(5) Based on the predicted remote sensing image, five remote sensing indexes of vegetation coverage VC, vegetation health degree VHI, soil humidity Wet, earth surface civil engineering degree NDBSI and earth surface temperature LST are used as the earth surface environment comprehensive characterization of the selected area, and the numerical values of the five remote sensing indexes are calculated;
(6) Denoising and normalizing the remote sensing index;
(7) Obtaining the comprehensive index of the surface environment.
Further, the specific method for preprocessing the real photographing remote sensing image in the step (1) comprises the following steps: and downloading all the real photographing optical remote sensing image data subjected to geometric correction and radiation correction in the selected area in the research period, or manually carrying out geometric correction and radiation correction after downloading all the original real photographing optical remote sensing image data, and then removing invalid observation values covered by cloud, cloud shadow or snow in each scene image by utilizing the quality evaluation wave band of the remote sensing image.
Further, the specific establishment method of the step (2) comprises the following steps: the initial time sequence model of each pixel consists of a linear function and a harmonic function, three relational expressions are built from thick to thin according to the number of effective observation values of each pixel in a time sequence, and if the number of the effective observation values is less than or equal to 12 and less than 18, the model adopts the relational expression (1); if the number of the effective observed values is less than or equal to 18 and less than 24, the model adopts a relation (2); if the number of the effective observation values is more than or equal to 24, selecting a relational expression (3) for the model;
wherein ,
representing predicted surface reflectivity of the pixel; n represents the index of the relation selected by the initial time series model, n=1, 2,3; i is the ith wave band of the remote sensing image; t is julian day; w is the frequency, < >>
c 0,i and c1,i For intercept coefficient and slope coefficient, a n,i and bn,i Is the nth harmonic coefficient, and c 0,i 、c 1,i 、a n,i and bn,i The four parameters all determine the optimal value through a least square method.
Further, the specific detection method in the step (3) comprises the following steps: in order to accurately detect the change of the initial time series model of the pixel, the occurrence time of the change is determined by comparing the residual resid of the predicted value and the subsequent effective observed value of the initial time series model of the pixel with the root mean square error RMSE of the initial time series model, namely
Wherein: d is a band set of the remote sensing image, i is an ith band of the remote sensing image, and rest i RMSE for the residual of the initial time series model predictor and the subsequent effective observer for the ith band i The root mean square error for the initial time series model of the picture elements of the ith band,
the degree of freedom is the inverse chi-square distribution number of the image wave band number;
judging 6 continuous subsequent effective observation values of the same pixel in the same wave band as a group, if the ratio of the resid to the RMSE of the 6 subsequent effective observation values of the same pixel in all wave bands is greater than a threshold value, taking the occurrence time of the 1 st effective observation value in the subsequent 6 effective observation values of the pixel in each wave band as change occurrence time, taking the change occurrence time as a breakpoint in a time sequence model of the pixel in each wave band, segmenting an initial time sequence model of the pixel in each wave band, and obtaining a segmented time sequence model according to the relation in the step (2) so as to reflect the surface reflectivity change of the pixel;
if the ratio of the resid and the RMSE of the 6 subsequent effective observations of the same pixel in all the wave bands is not all greater than the threshold value, considering that the change does not occur, adding the 1 st effective observation of the pixel in the 6 subsequent effective observations in each wave band into an initial time sequence model of the corresponding wave band, and then fitting and updating the initial time sequence model of the pixel according to the relation in the step (2);
the above-described discrimination process is repeated until all of the valid observations are located in the initial time series model or the segmented time series model of the picture element.
Further, the specific method of the step (4) is as follows: and (3) positioning an initial time sequence model or a segmented time sequence model of the pixels band by band, determining the julian day t of the comprehensive evaluation of the surface environment according to the requirement, obtaining the predicted surface reflectivity of the pixels by the relational expression of the step (2), forming all the pixels into bands, forming all the bands into images, and obtaining the predicted remote sensing images at any moment.
Further, the five remote sensing indexes in the step (5): the relation among vegetation coverage VC, vegetation health VHI, soil humidity Wet, earth surface civil engineering degree NDBSI and earth surface temperature LST is as follows:
wherein ,ρNIR and ρRed Respectively predicting near infrared band and red light band of the remote sensing image, wherein NDVI is normalized difference vegetation index, and NDVI veg and NDVIsoil Normalized difference vegetation indexes of vegetation and soil are respectively the maximum and minimum values of NDVI in a research range;
VHI＝PCA1(NDVI,NRI,NDSVI) (9)
wherein ,ρSWIR1 and ρGreen Respectively predicting the 1 st wave band and the green light wave band of the short wave infrared of the remote sensing image, wherein NRI is the reflectivity index of nitrogen element,NDSVI is normalized differential vegetation decay index, vegetation health VHI is first component PCA1 by transforming NDVI, NRI, NDSVI principal component;
wherein D is a remote sensing image wave band set, i is a predicted remote sensing image ith wave band, a i Is the linear combination coefficient of wave bands, ρ i The reflection rate of the image ground surface of the ith wave band;
wherein ,ρBlue For predicting blue light wave bands of the remote sensing image, IBI is a built-up area index, and SI is a soil index;
where h is the planck constant, h=6.26×10 -34 J-sec; c is the speed of light, c= 2.998 ×10 8 m/sec; k is the stent-boltzmann constant, k=1.38x10 -23 J/K; alpha is an intermediate variable, lambda is the central wavelength of a thermal infrared band, and epsilon is the surface emissivity; t (T) b Is the brightness temperature.
Further, the specific method of the step (6) is as follows:
adopting a percentile denoising method for the remote sensing index VC, VHI, wet, NDBSI, LST, setting a percentile as a, sequencing pixel values contained in all research periods from small to large, taking a% as an effective value lower limit, taking (100-a)% as an effective value upper limit, and taking the effective value lower limit if the pixel value is smaller than the effective value lower limit; if the pixel value is greater than the effective value upper limit, taking the effective value upper limit; then normalizing the five remote sensing indexes by a relational expression (16) to shield the influence of the difference of the numerical ranges of different remote sensing indexes;
wherein Index is the original remote sensing Index, index lb and Indexub Respectively a lower limit and an upper limit of a remote sensing Index, index nom Is the normalized remote sensing index.
Further, the specific method of the step (7) is as follows:
each remote sensing index after denoising normalization is subjected to minimum noise separation transformation, and a first component of each remote sensing index is taken as a preliminary surface environment comprehensive index LSECI 0 Observing the load value of each remote sensing index in the first component, and if the load value representing the vegetation coverage VC is greater than or equal to zero, LSECI 0 The final surface environment comprehensive index LSECI; if the load value representing vegetation coverage VC is negative, LSECI will be used 0 Obtaining a final surface environment comprehensive index LSECI by inverting, wherein the relation is as follows:
LSECI 0 ＝MNF1(VC,VHI,Wet,NDBSI,LST) (17)
where MNF1 is the first component of the minimum noise separation transform and q (VC) is the load value of vegetation coverage VC in the first component.
Compared with the prior art, the invention has the advantages that:
(1) According to the method, the predicted remote sensing image at any moment can be generated by constructing and detecting the time sequence model of the pixels, the bottleneck that the usability is poor due to the fact that a conventional optical image is covered by cloud and rain is overcome, and long-time and large-range continuous observation of coastal areas with frequent cloud and rain weather on the earth surface is possible.
(2) The comprehensive index LSECI of the earth surface environment is provided, five indexes of vegetation coverage, vegetation health, soil humidity, earth surface civil engineering and earth surface temperature are automatically and objectively integrated by adopting minimum noise separation transformation, the uncertainty of the evaluation result of artificial weighting of the indexes is solved, the evaluation result can be quantized in space, and the evaluation result has cross-region comparability.
(3) The method enables the periodic evaluation of the surface environment to be possible, effectively avoids the problem that the periodic evaluation is difficult to develop due to the lack of images, and simultaneously solves the problem that the evaluation results are difficult to compare due to the difference of imaging time of the images. The method is high in automation degree and strong in robustness, is suitable for combining remote sensing data of various platforms and different space and time resolutions, and is expected to expand comprehensive evaluation of the surface environment to a fine time scale.
Drawings
FIG. 1 is a block diagram illustrating a pixel timing model construction and change detection process according to an embodiment of the present invention, wherein:
(a 1-a 3) capturing land utilization/coverage variation (LUCC) by using a piecewise linear harmonic function model, taking the reflectivity of SWIR1 wave band as an example, the green and blue curves correspond to different time sequence model segments before and after the LUCC;
(b 1-b 3) are *** earth history snapshots obtained from the red rectangular marked period in (a 1-a 3), the yellow pin corresponding to the time sequence model and the pixels of the change detection;
(c 1-c 2) are predicted Landsat satellite images generated from model coefficients for 2010, 7, 1 and 2018, 7, 1;
(d 1-d 2) is a Globeland30 land utilization classification map, the red cross corresponding to the time series model and the pixels of the change detection.
Fig. 2 is a comparison of normalized NDVI by different methods, wherein:
(a 1-c 1) are the raw NDVI and histogram distributions thereof in 1987 and 2017;
(a 2-c 2) normalized NDVI in 1987 and 2017, calculated using the maximum and minimum values of a single year (marked by dotted lines in a 3) and their histogram distributions;
(a 3-c 3) are normalized NDVI from 1987 to 2017, generated after normalization with percent denoising (marked with dashed lines in a 3), and their histogram distribution.
FIG. 3 is a spatial variation of the environmental integrated quality index LSECI from 1986 to 2019 in Hangzhou, showing the annual LSECI for 7 months and 1 day of the corresponding year.
Fig. 4 is a spatial variation of the integrated quality index LSECI of half-month 2019 in hangzhou, showing that the 10 th and 20 th days of the corresponding month correspond to the first to fourth columns in spring, summer, autumn and winter, respectively.
Detailed Description
The invention is described in further detail below with reference to the embodiments of the drawings.
As shown in the figure, the comprehensive evaluation method of the earth surface environment at any moment based on time sequence remote sensing comprises the following specific steps:
(1) Downloading all the real photographing optical remote sensing image data subjected to geometric correction and radiation correction in the selected area in the research period, or downloading all the original real photographing optical remote sensing image data (such as Landsat, sentinel, MODIS and the like), then manually performing geometric correction and radiation correction, removing invalid observed values covered by cloud, cloud shadow or snow in each image by utilizing the quality evaluation wave band of the remote sensing image, and orderly arranging all the preprocessed images according to the imaging time sequence to form an image time sequence;
(2) Respectively establishing an initial time sequence model for each pixel in the image on all wave bands according to the time sequence, wherein the initial time sequence model comprises the following concrete steps: the initial time sequence model of each pixel consists of a linear function and a harmonic function, three relational expressions are built from thick to thin according to the number of effective observation values of each pixel in a time sequence, and if the number of the effective observation values is less than or equal to 12 and less than 18, the model adopts the relational expression (1); if the number of the effective observed values is less than or equal to 18 and less than 24, the model adopts a relation (2); if the number of the effective observation values is more than or equal to 24, selecting a relational expression (3) for the model;
wherein ,
representing predicted surface reflectivity of the pixel; n represents the index of the relation selected by the initial time series model, n=1, 2,3; i is the ith wave band of the remote sensing image; t is julian day; w is the frequency, < >>
c 0,i and c1,i Representing gradual annual changes for intercept coefficients and slope coefficients; a, a n,i and bn,i Is the nth harmonic coefficient representing periodic intra-annual changes, and c 0,i 、c 1,i 、a n,i and bn,i The four parameters all determine the optimal value through a least square method;
(3) And detecting the change of the initial time sequence model of the pixel:
in order to accurately detect the change of the initial time series model of the pixel, the occurrence time of the change is determined by comparing the residual resid of the predicted value and the subsequent effective observed value of the initial time series model of the pixel with the root mean square error RMSE of the initial time series model, namely
Wherein: d is a band set of the remote sensing image, i is an ith band of the remote sensing image, and rest i RMSE for the residual of the initial time series model predictor and the subsequent effective observer for the ith band i The root mean square error for the initial time series model of the picture elements of the ith band,
the degree of freedom is the inverse chi-square distribution number of the image wave band number;
judging 6 continuous subsequent effective observation values of the same pixel in the same wave band as a group, if the ratio of the resid to the RMSE of the 6 subsequent effective observation values of the same pixel in all wave bands is greater than a threshold value, taking the occurrence time of the 1 st effective observation value in the subsequent 6 effective observation values of the pixel in each wave band as change occurrence time, taking the change occurrence time as a breakpoint in a time sequence model of the pixel in each wave band, segmenting an initial time sequence model of the pixel in each wave band, and obtaining a segmented time sequence model according to the relation in the step (2) so as to reflect the surface reflectivity change of the pixel;
if the ratio of the resid and the RMSE of the 6 subsequent effective observations of the same pixel in all the wave bands is not all greater than the threshold value, considering that the change does not occur, adding the 1 st effective observation of the pixel in the 6 subsequent effective observations in each wave band into an initial time sequence model of the corresponding wave band, and then fitting and updating the initial time sequence model of the pixel according to the relation in the step (2);
the above-mentioned distinguishing process is repeated until all the effective observed values are positioned in the original time sequence model or the sectional time sequence model of the pixel;
(4) On the basis of a time sequence model of the pixels, a predicted remote sensing image at any moment is obtained, specifically:
the initial time sequence model or the sectional time sequence model of the pixels is positioned band by band and pixel by pixel, the julian day t of the comprehensive evaluation of the earth surface environment is determined according to the research requirement, t can be any time in the research period, the predicted earth surface reflectivity of the pixels is obtained through the relation of the step (2), then all the pixels are formed into bands, all the bands are formed into images, and the predicted remote sensing images at any moment are obtained;
(5) Based on the predicted remote sensing image, five remote sensing indexes of vegetation coverage VC, vegetation health degree VHI, soil humidity Wet, earth surface civil engineering degree NDBSI and earth surface temperature LST are used as the earth surface environment comprehensive characterization of the selected area, and the numerical values of the five remote sensing indexes are calculated, wherein the specific relation is as follows:
wherein ,ρNIR and ρRed Respectively predicting near infrared band and red light band of the remote sensing image, wherein NDVI is normalized difference vegetation index, and NDVI veg and NDVIsoil Normalized difference vegetation indexes of vegetation and soil are respectively the maximum and minimum values of NDVI in a research range;
VHI＝PCA1(NDVI,NRI,NDSVI) (9)
wherein ,ρSWIR1 and ρGreen Respectively predicting short wave infrared 1 st wave band and green light wave band of the remote sensing image, wherein NRI is reflectivity index of nitrogen element, NDSVI is normalized difference vegetation decay index, and vegetation health degree VHI is obtained by transforming NDVI, NRI, NDSVI main componentA first component PCA1;
wherein D is a remote sensing image wave band set, i is a predicted remote sensing image ith wave band, a i Is the linear combination coefficient of wave bands, ρ i The reflection rate of the image ground surface of the ith wave band;
wherein ,ρBlue For predicting blue light wave bands of the remote sensing image, IBI is a built-up area index, and SI is a soil index;
where h is the planck constant, h=6.26×10 -34 J-sec; c is the speed of light, c= 2.998 ×10 8 m/sec; k is the stent-boltzmann constant, k=1.38x10 -23 J/K; alpha is an intermediate variable, lambda is the central wavelength of a thermal infrared band, and epsilon is the surface emissivity; t (T) b Is the brightness temperature;
(6) The remote sensing index VC, VHI, wet, NDBSI, LST is denoised by adopting a percentile denoising method, specifically: setting a percentile as a, sorting pixel values contained in all research periods from small to large, taking a% as an effective value lower limit, taking (100-a)% as an effective value upper limit, and taking the effective value lower limit if the pixel value is smaller than the effective value lower limit; if the pixel value is greater than the effective value upper limit, taking the effective value upper limit; then normalizing the five remote sensing indexes by a relational expression (16) to shield the influence of the difference of the numerical ranges of different remote sensing indexes;
wherein Index is the original remote sensing Index, index lb and Indexub Respectively a lower limit and an upper limit of a remote sensing Index, index nom The normalized remote sensing index;
(7) The denoised and normalized remote sensing indexes are subjected to minimum noise separation transformation (MNF) and the first component is taken as a preliminary surface environment integrated index LSECI 0 Observing the load value of each remote sensing index in the first component, and if the load value representing the vegetation coverage VC is greater than or equal to zero, LSECI 0 The final surface environment comprehensive index LSECI; if the load value representing vegetation coverage VC is negative, LSECI will be used 0 Obtaining a final surface environment comprehensive index LSECI by inverting, wherein the relation is as follows:
LSECI 0 ＝MNF1(VC,VHI,Wet,NDBSI,LST) (17)
where MNF1 is the first component of the minimum noise separation transform and q (VC) is the load value of vegetation coverage VC in the first component.
The results of the method are verified as follows.
Example one:
the Hangzhou city is taken as an experimental area, land resource satellite (Landsat) images are taken as data sources, and a comprehensive evaluation result of the annual (7 months and 1 day each year) surface environment is generated, which is specifically as follows:
(1) The L2C2 Landsat image products published by USGS of the United states geological survey were acquired, each image had been systematically geometrically and radioactively corrected, and the Landsat TM/ETM+/OLI image products of 1985-2020 were collected for 4348 scenes in total. Each Landsat image uses 8 bands including 6 optical bands (blue, green, red, near infrared, short wave infrared 1 and short wave infrared 2), 1 thermal infrared band and 1 quality assessment band. In the quality evaluation band, the values 3,4,5 represent the coverage of the earth surface by cloud, snow, cloud shadows, respectively. Taking the method as a basis, removing pixels covered by cloud, snow and cloud shadows in each image, and sequentially arranging the images according to imaging time sequence to form an image time sequence;
(2) For the optical band and the thermal infrared band, an initial time sequence model of the pixels is built pixel by pixel:
firstly, collecting 24 effective observation data, establishing an initial time sequence model by using a relation (3), and determining the values of four parameters of a time sequence model by using a least square method, wherein the values are shown in the figure 1;
(3) The initial time sequence model of the pixel is subjected to change detection, and in consideration of the short wavelength of the blue light wave band and high potential noise, other 5 optical wave bands (green light, red light, near infrared, short wave infrared 1 and short wave infrared 2) are taken as model change detection wave bands in the embodiment, so that the threshold V of change detection inv_x2 =15.1. Sequentially taking the subsequent 6 effective observation values as a group, detecting the change by using a relation (4), taking the time of occurrence of the 1 st effective observation value in the subsequent 6 effective observation values as change occurrence time if the change is detected, taking the change occurrence time as a breakpoint in a time sequence model of the pixel, and segmenting an initial time sequence model to obtain a segmented time sequence model; otherwise, adding the 1 st of the following 6 effective observations into the initial time sequence model, fitting and updating the initial time sequence model of the pixel, as shown in figure 1, and repeating the process until all the effective observations are traversed;
(4) Obtaining a predicted remote sensing image at any moment: the example assumes that 7 month 1 day of each year in 1986-2019 is ground surface environment comprehensive evaluation time, converts the time into corresponding julian days, locates a time sequence model corresponding to the julian days above each band and each pixel, obtains the reflectivity of the pixel ground surface through the relation (1-3), and on the basis, forms all pixels into bands, forms all bands into images, see figure 1, and obtains the 7 month 1 day predictive remote sensing images of each year in 1986-2019 through the above processes;
(5) Remote sensing index calculation: calculating vegetation coverage VC, vegetation health degree VHI, soil humidity Wet, earth surface civil engineering degree NDBSI and earth surface temperature LST on the basis of the 34-scene predicted remote sensing image scene by scene, wherein the relation of the soil humidity Wet suitable for Landsat images is as follows:
WET＝0.1484ρ Blue +0.3068ρ Green +0.2437ρ Red +0.1886ρ NIR -0.7184ρ SWIR1 -0.5352ρ SWIR2 (19)
in calculating the surface temperature LST, the land utilization type of each pixel is determined by using the Globeland30 global land cover product, see fig. 1, and the corresponding surface emission rate epsilon is set: 0.97-forest, 0.91-grassland, 0.92-artificial ground surface, 0.95-bare land, 0.99-water body and 0.97-other vegetation;
(6) Remote sensing index denoising and normalization: the relation (16) is utilized to carry out percentile denoising normalization on the images of all periods of each remote sensing index, the percentile a is set to be 0.5, namely the effective interval of the pixel values of the remote sensing indexes is [0.5%,99.5% ], and compared with the conventional maximum and minimum normalization method, the percentile denoising normalization method can better ensure the relative consistency of the distribution of the values before and after normalization, and is shown in the figure 2;
(7) Synthesizing the comprehensive index of the surface environment: for each remote sensing index of each period, carrying out minimum noise separation transformation by using a relation (17) to obtain a first component serving as a preliminary surface environment comprehensive index LSECI 0 The contribution degree of the first component to the whole information of the remote sensing indexes is 95.3%, which shows that the first component can comprehensively represent the earth surface condition represented by each remote sensing index, and in the first component, the load of vegetation coverage VC is minus 0.632, so that the relation (18) is utilized to invert the preliminary earth surface environment comprehensive index to obtain the final earthTable environment synthesis index LSECI. The comprehensive evaluation result of the surface environment, which is obtained by adopting the above process and is shown in the figure 3, is 7 months and 1 days each year, and the larger the numerical value is, the better the surface environment quality is represented.
Example two:
and (3) generating a comprehensive surface environment assessment result of the half-moon of 2019 (10 days and 20 days of each month) by taking Hangzhou city as an experimental area and taking a land resource satellite (Landsat) image as a data source.
The remaining steps are the same as in example one, except that: in the step (4), the example assumes that 10 days and 20 days of each month in 2019 are ground surface environment comprehensive evaluation time, converts the comprehensive evaluation time into corresponding julian days, positions a time sequence model corresponding to the julian days above by each band and each pixel, obtains the reflectivity of the ground surface of the pixel through the relation (1-3), then forms all the pixels into bands, forms all the bands into images, and obtains 24 scenes of predicted remote sensing images of 10 days and 20 days of each month in 2019; the final comprehensive evaluation result of the surface environment in the half-moon of 2019 is shown in fig. 4, which further proves that the invention can perform comprehensive evaluation of the surface environment at any time.
The protection scope of the present invention includes, but is not limited to, the above embodiments, the protection scope of which is subject to the claims, and any substitutions, modifications, and improvements made by those skilled in the art are within the protection scope of the present invention.
Claims (6)
1. A comprehensive evaluation method of an earth surface environment at any moment based on time sequence remote sensing is characterized by comprising the following specific steps:
(1) Preprocessing the real-shot remote sensing images of the selected area in the research period, and orderly arranging all preprocessed images according to the imaging time sequence to form an image time sequence;
(2) Respectively establishing an initial time sequence model for each pixel in the image on all wave bands according to the time sequence, wherein the initial time sequence model comprises the following concrete steps: the initial time sequence model of each pixel consists of a linear function and a harmonic function, three relational expressions are built from thick to thin according to the number of effective observation values of each pixel in a time sequence, and if the number of the effective observation values is less than or equal to 12 and less than 18, the model adopts the relational expression (1); if the number of the effective observed values is less than or equal to 18 and less than 24, the model adopts a relation (2); if the number of the effective observation values is more than or equal to 24, selecting a relational expression (3) for the model;
wherein ,
representing predicted surface reflectivity of the pixel; n represents the index of the relation selected by the initial time series model, n=1, 2,3; i is the ith wave band of the remote sensing image; t is julian day; w is the frequency, < >>
c 0,i and c1,i For intercept coefficient and slope coefficient, a n,i and bn,i Is the nth harmonic coefficient, and c 0,i 、c 1,i 、a n,i and bn,i The four parameters all determine the optimal value through a least square method;
(3) Detecting the change of the initial time sequence model of the pixel, and if the change does not occur, maintaining the initial time sequence model; if the change occurs, taking the change occurrence time as a breakpoint of the time sequence model, segmenting the initial time sequence model of the pixel to obtain a segmented time sequence model so as to reflect the change of the earth surface reflectivity of the corresponding pixel; the method comprises the following steps:
in order to accurately detect the change of the initial time series model of the pixel, the occurrence time of the change is determined by comparing the residual resid of the predicted value and the subsequent effective observed value of the initial time series model of the pixel with the root mean square error RMSE of the initial time series model, namely
Wherein: d is a band set of the remote sensing image, i is an ith band of the remote sensing image, and rest i RMSE for the residual of the initial time series model predictor and the subsequent effective observer for the ith band i The root mean square error for the initial time series model of the picture elements of the ith band,
the degree of freedom is the inverse chi-square distribution number of the image wave band number;
judging 6 continuous subsequent effective observation values of the same pixel in the same wave band as a group, if the ratio of the resid to the RMSE of the 6 subsequent effective observation values of the same pixel in all wave bands is greater than a threshold value, taking the occurrence time of the 1 st effective observation value in the subsequent 6 effective observation values of the pixel in each wave band as change occurrence time, taking the change occurrence time as a breakpoint in a time sequence model of the pixel in each wave band, segmenting an initial time sequence model of the pixel in each wave band, and obtaining a segmented time sequence model according to the relation in the step (2) so as to reflect the surface reflectivity change of the pixel;
if the ratio of the resid and the RMSE of the 6 subsequent effective observations of the same pixel in all the wave bands is not all greater than the threshold value, considering that the change does not occur, adding the 1 st effective observation of the pixel in the 6 subsequent effective observations in each wave band into an initial time sequence model of the corresponding wave band, and then fitting and updating the initial time sequence model of the pixel according to the relation in the step (2);
the above-mentioned distinguishing process is repeated until all the effective observed values are positioned in the original time sequence model or the sectional time sequence model of the pixel;
(4) Obtaining a predicted remote sensing image at any moment on the basis of an initial time sequence model or a segmented time sequence model of the pixels;
(5) Based on the predicted remote sensing image, five remote sensing indexes of vegetation coverage VC, vegetation health degree VHI, soil humidity Wet, earth surface civil engineering degree NDBSI and earth surface temperature LST are used as the earth surface environment comprehensive characterization of the selected area, and the numerical values of the five remote sensing indexes are calculated;
(6) Denoising and normalizing the remote sensing index;
(7) Obtaining the comprehensive index of the surface environment.
2. The method for comprehensively evaluating the surface environment at any moment based on time sequence remote sensing as claimed in claim 1, wherein the specific method for preprocessing the real photographing remote sensing image in the step (1) is as follows: and downloading all the real photographing optical remote sensing image data subjected to geometric correction and radiation correction in the selected area in the research period, or manually carrying out geometric correction and radiation correction after downloading all the original real photographing optical remote sensing image data, and then removing invalid observation values covered by cloud, cloud shadow or snow in each scene image by utilizing the quality evaluation wave band of the remote sensing image.
3. The method for comprehensively evaluating the surface environment at any moment based on time sequence remote sensing as claimed in claim 1, wherein the specific method in the step (4) is as follows: and (3) positioning an initial time sequence model or a segmented time sequence model of the pixels band by band, determining the julian day t of the comprehensive evaluation of the surface environment according to the research requirement, obtaining the predicted surface reflectivity of the pixels by the relational expression of the step (2), forming all the pixels into bands, forming all the bands into images, and obtaining the predicted remote sensing images at any moment.
4. The method for comprehensively evaluating the surface environment at any moment based on time sequence remote sensing as claimed in claim 1, wherein the five remote sensing indexes in the step (5) are as follows: the relation among vegetation coverage VC, vegetation health VHI, soil humidity Wet, earth surface civil engineering degree NDBSI and earth surface temperature LST is as follows:
wherein ,ρNIR and ρRed Respectively predicting near infrared band and red light band of the remote sensing image, wherein NDVI is normalized difference vegetation index, and NDVI veg and NDVIsoil Normalized difference vegetation indexes of vegetation and soil are respectively the maximum and minimum values of NDVI in a research range;
VHI＝PCA1(NDVI,NRI,NDSVI) (9)
wherein ,ρSWIR1 and ρRGreen Respectively predicting a short wave infrared 1 st wave band and a green light wave band of a remote sensing image, wherein NRI is a reflectivity index of nitrogen element, NDSVI is a normalized difference vegetation decay index, and vegetation health degree VHI is a first component PCA1 obtained by transforming NDVI, NRI, NDSVI main components;
wherein D is a remote sensing image wave band set, i is a predicted remote sensing image ith wave band, a i Is the linear combination coefficient of wave bands, ρ i The reflection rate of the image ground surface of the ith wave band;
wherein ,ρBlue For predicting blue light wave bands of the remote sensing image, IBI is a built-up area index, and SI is a soil index;
where h is the planck constant, h=6.26×10 -34 J-sec; c is the speed of light, c= 2.998 ×10 8 m/sec; k is the stent-boltzmann constant, k=1.38x10 -23 J/K; alpha is an intermediate variable, lambda is the central wavelength of a thermal infrared band, and epsilon is the surface emissivity; t (T) b Is the brightness temperature.
5. The method for comprehensively evaluating the surface environment at any moment based on time sequence remote sensing as claimed in claim 1, wherein the specific method in the step (6) is as follows:
adopting a percentile denoising method for the remote sensing index VC, VHI, wet, NDBSI, LST, setting a percentile as a, sequencing pixel values contained in all research periods from small to large, taking a% as an effective value lower limit, taking (100-a)% as an effective value upper limit, and taking the effective value lower limit if the pixel value is smaller than the effective value lower limit; if the pixel value is greater than the effective value upper limit, taking the effective value upper limit; then normalizing the five remote sensing indexes by a relational expression (16) to shield the influence of the difference of the numerical ranges of different remote sensing indexes;
wherein Index is the original remote sensing Index, index lb and Indexub Respectively a lower limit and an upper limit of a remote sensing Index, index nom Is the normalized remote sensing index.
6. The method for comprehensively evaluating the surface environment at any moment based on time sequence remote sensing as claimed in claim 1, wherein the specific method of the step (7) is as follows:
each remote sensing index after denoising normalization is subjected to minimum noise separation transformation, and a first component of each remote sensing index is taken as a preliminary surface environment comprehensive index LSECI 0 Observing the load value of each remote sensing index in the first component, and if the load value representing the vegetation coverage VC is greater than or equal to zero, LSECI 0 The final surface environment comprehensive index LSECI; if the load value representing vegetation coverage VC is negative, LSECI will be used 0 Obtaining a final surface environment comprehensive index LSECI by inverting, wherein the relation is as follows:
LSECI 0 ＝MNF1(VC,VHI,Wet,NDBSI,LST) (17)
where MNF1 is the first component of the minimum noise separation transform and q (VC) is the load value of vegetation coverage VC in the first component.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210810131.0A CN115203624B (en) | 2022-07-11 | 2022-07-11 | Comprehensive evaluation method for surface environment at any moment based on time sequence remote sensing |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210810131.0A CN115203624B (en) | 2022-07-11 | 2022-07-11 | Comprehensive evaluation method for surface environment at any moment based on time sequence remote sensing |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115203624A CN115203624A (en) | 2022-10-18 |
| CN115203624B true CN115203624B (en) | 2023-06-16 |
Family
ID=83579312
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210810131.0A Active CN115203624B (en) | 2022-07-11 | 2022-07-11 | Comprehensive evaluation method for surface environment at any moment based on time sequence remote sensing |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115203624B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116466368B (en) * | 2023-06-16 | 2023-08-22 | 成都远望科技有限责任公司 | Dust extinction coefficient profile estimation method based on laser radar and satellite data |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110009213B (en) * | 2019-03-28 | 2021-10-22 | 交通运输部水运科学研究所 | Method for tracking, monitoring and evaluating ecological influence of channel engineering based on long-time sequence satellite remote sensing data |
| CN110472525B (en) * | 2019-07-26 | 2021-05-07 | 浙江工业大学 | Noise detection method for time series remote sensing vegetation index |
| CN113988626B (en) * | 2021-10-28 | 2022-04-12 | 中国矿业大学(北京) | Mining area ecological environment remote sensing comprehensive evaluation index realization method |
| CN114241334A (en) * | 2021-12-20 | 2022-03-25 | 南京大学 | Multi-factor integrated arid region ecological environment remote sensing monitoring method |
-
2022
- 2022-07-11 CN CN202210810131.0A patent/CN115203624B/en active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN115203624A (en) | 2022-10-18 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Berberoglu et al. | Assessing different remote sensing techniques to detect land use/cover changes in the eastern Mediterranean | |
| CN109308688B (en) | Visible light and near-infrared band thick cloud and shadow removing method | |
| CN111982822B (en) | Long-time sequence high-precision vegetation index improvement algorithm | |
| CN111582575B (en) | Method for identifying urban thermal environment formation development leading factors under multiple space-time scales | |
| CN115330159A (en) | Desert ecological restoration effect evaluation method and device, terminal equipment and storage medium | |
| CN115203624B (en) | Comprehensive evaluation method for surface environment at any moment based on time sequence remote sensing | |
| CN114819737B (en) | Method, system and storage medium for estimating carbon reserves of highway road vegetation | |
| CN115690596A (en) | MODIS data-based regional ecological environment quality comprehensive evaluation method | |
| CN112329790A (en) | Rapid extraction method for urban impervious surface information | |
| Solaimani et al. | Land use/cover change detection based on remote sensing data (A case study; Neka Basin) | |
| CN113887493A (en) | Black and odorous water body remote sensing image identification method based on ID3 algorithm | |
| Morakinyo et al. | Mapping of land cover and estimation of their emissivity values for gas flaring sites in the Niger Delta | |
| Lili et al. | Land Use/Cover Change from 2001 to 2010 and its Socioeconomic Determinants in Guangdong Province, A Rapid Urbanization Area of China. | |
| CN112215135A (en) | Mining area mining and treatment effect monitoring method and device | |
| CN107576399A (en) | Towards bright the temperature Forecasting Methodology and system of MODIS forest fire detections | |
| CN110287915B (en) | City impervious bed extraction method based on Landsat remote sensing image | |
| Xie et al. | Estimation of forest stand parameters using SPOT-5 satellite images and topographic information | |
| Nielsen et al. | The distribution in time and space of savanna fires in Burkina Faso as determined from NOAA AVHRR data | |
| Li et al. | New automated method for extracting river information using optimized spectral threshold water index | |
| Watson | Analysis of urban heat island climates along the I-85/I-40 corridor in central North Carolina | |
| Alexandridis et al. | LAI measurement with hemispherical photographs at variable conditions for assessment of remotely sensed estimations | |
| Kekula et al. | AN EMPIRICAL STUDY OF RELATIONSHIPS BETWEEN URBAN LIGHTING INDICATORS AND NIGHT-TIME LIGHT RADIANCE. | |
| Walder et al. | Land use change impact on urban land surface temperatures: A GIS-supported satellite-based case study | |
| CN117196922A (en) | Remote sensing method for continuously evaluating ecological quality of natural protection area | |
| Horton et al. | Using snow depth observation to provide insight into the quality of regional-scale snowpack simulations for avalanche forecasting |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |