CN115203624A - Time sequence remote sensing-based comprehensive evaluation method for earth surface environment at any moment - Google Patents
Time sequence remote sensing-based comprehensive evaluation method for earth surface environment at any moment Download PDFInfo
- Publication number
- CN115203624A CN115203624A CN202210810131.0A CN202210810131A CN115203624A CN 115203624 A CN115203624 A CN 115203624A CN 202210810131 A CN202210810131 A CN 202210810131A CN 115203624 A CN115203624 A CN 115203624A
- Authority
- CN
- China
- Prior art keywords
- remote sensing
- pixel
- index
- time sequence
- wave band
- 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.)
- Granted
Links
- 238000011156 evaluation Methods 0.000 title claims abstract description 36
- 238000000034 method Methods 0.000 claims abstract description 46
- 230000008859 change Effects 0.000 claims abstract description 45
- 238000011160 research Methods 0.000 claims abstract description 18
- 238000003384 imaging method Methods 0.000 claims abstract description 7
- 238000007781 pre-processing Methods 0.000 claims abstract description 5
- 238000002310 reflectometry Methods 0.000 claims description 20
- 239000002689 soil Substances 0.000 claims description 17
- 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
- 238000000926 separation method Methods 0.000 claims description 8
- 230000005855 radiation Effects 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
- 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
- 238000012163 sequencing technique Methods 0.000 claims description 3
- 230000001131 transforming effect Effects 0.000 claims description 3
- 238000011895 specific detection Methods 0.000 claims description 2
- 230000000737 periodic effect Effects 0.000 abstract description 6
- 238000001514 detection method Methods 0.000 description 7
- 238000010606 normalization Methods 0.000 description 5
- 230000007613 environmental effect Effects 0.000 description 4
- 238000011161 development Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 238000003915 air pollution Methods 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 238000010276 construction 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
- 238000001228 spectrum 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-shot remote sensing images of the selected area in the research period, and arranging all the 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 representation of the earth 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 method enables the periodic evaluation work of the earth surface environment to be possible, effectively avoids the problem that the periodic evaluation is difficult to carry out due to image loss, and solves the problem that the evaluation result is difficult to compare due to image imaging time difference.
Description
Technical Field
The invention relates to a remote sensing evaluation method of an earth surface environment, in particular to an arbitrary time earth surface environment comprehensive evaluation method based on time sequence remote sensing.
Background
The surface environment is highly relevant to the quality of human life and the sustainable development of the society today. Since half a century, the rapid urbanization process in China has led to a series of environmental problems, such as deforestation, loss of diversity, air and water pollution, land desertification, urban heat island effect, and the like. Since urbanization is an irreversible process in the development process of human civilization, the stress effect of urbanization on the ecological environment is difficult to relieve in a short period of time. Therefore, timely and objectively exploring the change situation of the earth surface environment in the urbanization process is particularly important for formulating a social economic development plan coordinated with the regional environmental capacity, which 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 spectrum information, convenience in data acquisition and the like, and is increasingly becoming an important means for representing surface environments at home and abroad. Early descriptions of the earth's surface environment focused primarily on remote sensing-based computation of a single index, primarily on some particular aspect of the earth's surface environment. Typical examples are Normalized Difference Vegetation Index (NDVI) and surface temperature (LST) often used to describe vegetation coverage and urban heat island conditions. In most cases, the surface environment changes are complex and difficult to represent by a single remote sensing index. Therefore, some scholars propose model schemes for integrating various remote sensing-based indexes, such as Forest Disturbance Index (FDI), global disturbance index (MGDI), ecological Niche Modeling (ENM), and the like. Although these composite indices can more fully reflect the surface environmental conditions, a series of challenges, such as weight selection, still exist. In fact, it is very difficult to determine the magnitude of the effect of each index on the surface environment change and to allocate the weights in proportion, and the weight allocation is difficult to avoid different subjective experience differences, and different weights will directly influence the surface environment evaluation result. In addition, due to the fact that a series of indexes (such as greenness and heat) are involved, the comprehensive index model scheme is very sensitive to the weather in the year, most researches adopt data with good quality and similar seasons in an annual image or perform small-range researches, but available data is often lost in coastal urban areas with high cloud amount, and the accuracy and comparability of the application of the method and the range and the time sequence length of the researched areas are severely limited.
Disclosure of Invention
The invention aims to solve the technical problem of providing a comprehensive evaluation method of the earth surface environment at any moment based on time sequence remote sensing, which has high automation degree and strong robustness, enables the evaluation of the earth surface environment at any moment to be possible, effectively avoids the problem that the periodic evaluation is difficult to be carried out due to image loss, and solves the problem that the evaluation result is difficult to compare due to the difference of the imaging time of the images.
The technical scheme adopted by the invention for solving the technical problems is as follows: a comprehensive evaluation method for any time earth surface environment 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 sequentially 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;
(3) Detecting the change of the initial time sequence model of the pixel, and if the change does not occur, keeping the initial time sequence model; if the change occurs, the time of the change occurrence is taken as a breakpoint of the time sequence model, and the initial time sequence model of the pixel is segmented to obtain a segmented time sequence model so as to reflect the change of the earth surface reflectivity of the corresponding pixel;
(4) On the basis of the time sequence model of the pixel, obtaining a predicted remote sensing image at any moment;
(5) Based on the predicted remote sensing image, taking five remote sensing indexes of vegetation coverage VC, vegetation health VHI, soil humidity Wet, surface civil engineering NDBSI and surface temperature LST as surface environment comprehensive representation of the selected area, and calculating numerical values of the five remote sensing indexes;
(6) Denoising and normalizing the remote sensing index;
(7) And obtaining the comprehensive index of the surface environment.
Further, the specific method for preprocessing the real shot remote sensing image in the step (1) comprises the following steps: and downloading all real-shot optical remote sensing image data which are subjected to geometric correction and radiation correction in a selected area in a research period, or downloading all original real-shot optical remote sensing image data, and then manually carrying out geometric correction and radiation correction, and then removing invalid observed values covered by clouds, cloud shadows or snow in each image by using a quality evaluation waveband 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 relations are constructed from coarse to fine according to the number of effective observed values of each pixel in a time sequence, and if the number of the effective observed values is not less than 12 and is less than 18, the model selects the relation (1); if the number of the effective observed values is more than or equal to 18 and less than 24, selecting a relational expression (2) for the model; if the number of the effective observed values is more than or equal to 24, selecting a relational expression (3) for the model;
wherein ,
representing the 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 of the signal to be detected,
c 0,i and c1,i Is an intercept coefficient and a 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 are determined to be optimal values by a least square method.
Further, the specific detection method in the step (3) is as follows: 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 error resid of the initial time series model predicted value and the subsequent effective observed value of the pixel with the root mean square error RMSE of the initial time series model, i.e.
Wherein: d is the wave band set of the remote sensing image, i is the ith wave band of the remote sensing image, and resid i Is the residual error between the initial time series model predicted value and the subsequent effective observed value of the ith wave band, RMSE i Is the root mean square error of the initial time series model of the picture elements of the ith band,
is a threshold value, and is equal to the inverse Chi-square distribution number of the degree of freedom which is the number of image wave bands;
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 the wave bands is larger than a threshold value, taking the time of the occurrence of the 1 st effective observation value of the pixel in the subsequent 6 effective observation values of each wave band as the change occurrence time, taking the change occurrence time as a breakpoint of the pixel in the time sequence model of each wave band, segmenting the initial time sequence model of the pixel in each wave band, and obtaining a segmented time sequence model according to the relational expression in the step (2) to reflect the surface reflectivity change of the pixel;
if the ratio of the resid to the RMSE of the 6 subsequent effective observation values of the same pixel in all the wave bands is not all larger than the threshold value, the change is considered to be not generated, the 1 st effective observation value of the pixel in the 6 subsequent effective observation values in each wave band is added into the initial time sequence model of the corresponding wave band, and then the initial time sequence model of the pixel is fitted and updated according to the relational expression in the step (2);
the above-mentioned discrimination process is repeated until all the effective observed values are located in the initial time series model or the segmented time series model of the pixel.
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 on a wave band-by-pixel basis, determining the julian day t of the comprehensive evaluation of the earth surface environment according to the requirement, obtaining the predicted earth surface reflectivity of the pixels through the relational expression in the step (2), forming all the pixels into wave bands, forming all the wave bands into images, and obtaining the predicted remote sensing image at any moment.
Further, the five remote sensing indexes in the step (5) are as follows: the relation formula of the vegetation coverage VC, the vegetation health VHI, the soil humidity Wet, the surface civil engineering NDBSI and the surface temperature LST is as follows:
wherein ,ρNIR and ρRed Respectively near infrared band and red light band of the predicted remote sensing image, and NDVI is normalized difference vegetation index, NDVI veg and NDVIsoil The normalized difference vegetation indexes of the vegetation and the soil are respectively the maximum value and the minimum value of the NDVI in the research range;
VHI＝PCA1(NDVI,NRI,NDSVI) (9)
wherein ,ρSWIR1 and ρGreen 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 VHI is a first component PCA1 obtained by transforming main components of NDVI, NRI and NDSVI;
wherein D is the remote sensing image wave band set, i is the ith wave band of the predicted remote sensing image, a i Is a band linear combination coefficient, rho i The image surface reflectivity of the ith wave band is obtained;
wherein ,ρBlue In order to predict the blue light wave band of the remote sensing image, IBI is a built-up area index, and SI is a soil index;
wherein h is Planck constant, h =6.26 × 10 -34 J-sec; c is the speed of light, c =2.998 × 10 8 m/sec; k is Stefan-Boltzmann constant, K = 1.38X 10 -23 J/K; alpha is an intermediate variable, lambda is the central wavelength of a thermal infrared band, and epsilon is the ground surface emergence rate; t is 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 indexes VC, VHI, wet, NDBSI and LST, setting the percentile as a, sequencing pixel values contained in all research periods from small to large, taking a% as the lower limit of an effective value, taking (100-a)% as the upper limit of the effective value, and taking the lower limit of the effective value if the pixel values are smaller than the lower limit of the effective value; if the pixel value is larger than the upper limit of the effective value, the upper limit of the effective value is taken; then, normalizing the five remote sensing indexes through a relational expression (16) so as to shield the influence of the difference of numerical ranges of different remote sensing indexes;
wherein, index is the original remote sensing Index, index lb and Indexub Lower and upper limits, index, respectively, of the remote sensing Index nom Is a normalized remote sensing index.
Further, the specific method of the step (7) is as follows:
carrying out minimum noise separation transformation on each denoised and normalized remote sensing index, and taking the first component as a primary 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 more than or equal to zero, then LSECI 0 The final surface environment comprehensive index LSECI; if the load value representing the vegetation coverage VC is negative, LSECI is carried out 0 Obtaining a final surface environment comprehensive index LSECI by taking the inverse, wherein the relational expression 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, through the construction and change detection of the time sequence model of the pixel, the predicted remote sensing image at any time can be generated, the bottleneck that the usability of a conventional optical image is poor due to cloud and rain coverage is overcome, and the long-time large-range continuous observation of the earth surface in coastal areas with frequent cloud and rain weather becomes possible.
(2) The method provides a comprehensive index LSECI of the earth surface environment, automatically and objectively synthesizes five indexes of vegetation coverage, vegetation health, soil humidity, earth surface civil engineering degree and earth surface temperature by adopting minimum noise separation transformation, solves the uncertainty of the evaluation result of the artificial empowerment of the indexes, can quantify the evaluation result in space and enables the evaluation result to have cross-region comparability.
(3) The method enables the periodic evaluation work of the earth surface environment to be possible, effectively avoids the problem that the periodic evaluation is difficult to carry out due to image loss, and solves the problem that the evaluation result is difficult to compare due to image imaging time difference. The method is high in automation degree and 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 earth surface environment to a fine time scale.
Drawings
Fig. 1 is a process of constructing a pixel timing model and detecting a change according to a first embodiment of the present invention, in which:
(a 1-a 3) capturing land use/cover change (LUCC) by utilizing a piecewise linear harmonic function model, taking the reflectivity of a SWIR1 wave band as an example, green and blue curves correspond to different time sequence model segments before and after the LUCC;
(b 1-b 3) is a Google Earth historical snapshot obtained from the period of the red rectangular mark in (a 1-a 3), and a yellow pin corresponds to a time sequence model and a pixel for change detection;
(c 1-c 2) are predicted Landsat satellite images generated by the model coefficients on 7/month 1 in 2010 and on 7/month 1 in 2018;
and (d 1-d 2) are Globeland30 land use classification graphs, and the red crosses correspond to time sequence models and pixels for change detection.
FIG. 2 is a comparison of normalized NDVI by different methods, wherein:
(a 1-c 1) original NDVI and histogram distribution thereof in 1987 and 2017;
(a 2-c 2) normalized NDVI in 1987 and 2017, calculated using maximum and minimum values (dashed line in a 3) for individual year and their histogram distribution;
(a 3-c 3) are normalized NDVI between 1987 and 2017, generated after normalization by percentage denoise (marked with a dashed line in a 3), and their histogram distribution.
FIG. 3 shows the spatial difference of the LSECI of the environmental quality index between 1986 and 2019 in Hangzhou city, and shows the LSECI of each year corresponding to 7 months and 1 day in each year.
Fig. 4 shows spatial differences of the overall quality index LSECI in 2019 semilunar in hangzhou city, showing that the 10 th and 20 th days of the corresponding months 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 accompanying examples.
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 real-shot optical remote sensing image data which are subjected to geometric correction and radiation correction in a selected area in a research period, or downloading all original real-shot optical remote sensing image data (such as Landsat, sentinel, MODIS and the like), then manually carrying out geometric correction and radiation correction, then utilizing the quality evaluation wave band of the remote sensing image to eliminate invalid observation values covered by clouds, cloud shadows or snow in each image, and then sequentially arranging all preprocessed images according to an 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, which specifically comprises the following steps: the initial time sequence model of each pixel consists of a linear function and a harmonic function, three relations are constructed from coarse to fine according to the number of effective observed values of each pixel in a time sequence, and if the number of the effective observed values is not less than 12 and is less than 18, the model selects the relation (1); if the number of the effective observed values is more than or equal to 18 and less than 24, selecting a relational expression (2) for the model; if the number of the effective observed values is more than or equal to 24, selecting a relational expression (3) for the model;
wherein ,
representing the predicted surface reflectivity of the pixel; n represents the initial timeThe label of the relation formula selected by the sequence model, n =1,2,3; i is the ith wave band of the remote sensing image; t is julian day; w is the frequency of the light emitted by the light source,
c 0,i and c1,i An intercept coefficient and a slope coefficient representing a gradual annual change; a is n,i and bn,i Is the nth harmonic coefficient, representing a periodic annual variation, and c 0,i 、c 1,i 、a n,i and bn,i The four parameters are determined to be optimal values by 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 error 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, that is
Wherein: d is the wave band set of the remote sensing image, i is the ith wave band of the remote sensing image, and resid i Is the residual error between the initial time series model predicted value and the subsequent effective observed value of the ith wave band, RMSE i Is the root mean square error of the initial time series model of the picture elements of the ith band,
is a threshold value, and is equal to the inverse Chi-square distribution number of the degree of freedom which is the number of image wave bands;
judging 6 continuous subsequent effective observed values of the same pixel in the same wave band as a group, if the ratio of the identification to the RMSE of the 6 subsequent effective observed values of the same pixel in all the wave bands is larger than a threshold value, taking the time of the pixel appearing in the 1 st effective observed value in the subsequent 6 effective observed values in each wave band as the change occurrence time, taking the change occurrence time as a breakpoint in a time sequence model of the pixel in each wave band, segmenting the initial time sequence model of the pixel in each wave band, and obtaining a segmented time sequence model according to the relational expression in the step (2) to reflect the change of the surface reflectivity of the pixel;
if the ratio of the resid to the RMSE of the 6 subsequent effective observation values of the same pixel in all the wave bands is not all larger than the threshold value, the change is considered to be not generated, the 1 st effective observation value of the pixel in the 6 subsequent effective observation values in each wave band is added into the initial time sequence model of the corresponding wave band, and then the initial time sequence model of the pixel is fitted and updated according to the relational expression in the step (2);
the judging process is repeated until all effective observed values are positioned in the initial time sequence model or the segmented time sequence model of the pixel;
(4) On the basis of the time series model of the pixel, a predicted remote sensing image at any moment is obtained, and the method specifically comprises the following steps:
positioning an initial time sequence model or a segmented time sequence model of pixels one by one wave band and one by one, determining a julian day t of comprehensive evaluation of the earth surface environment according to the research requirement, wherein t can be any time in the research period, obtaining the predicted earth surface reflectivity of the pixels through the relational expression in the step (2), then forming all the pixels into wave bands, forming all the wave bands into images, and obtaining the predicted remote sensing image at any moment;
(5) Based on the predicted remote sensing image, taking five remote sensing indexes of vegetation coverage VC, vegetation health VHI, soil humidity Wet, surface civil engineering NDBSI and surface temperature LST as the comprehensive representation of the surface environment of the selected area, and calculating the numerical values of the five remote sensing indexes, wherein the specific relational expression is as follows:
wherein ,ρNIR and ρRed Respectively near infrared band and red light band of the predicted remote sensing image, and NDVI is Normalized Difference Vegetation Index (NDVI) veg and NDVIsoil The normalized difference vegetation indexes of the vegetation and the soil are respectively the maximum value and the minimum value of the NDVI in the research range;
VHI＝PCA1(NDVI,NRI,NDSVI) (9)
wherein ,ρSWIR1 and ρGreen Respectively predicting a short wave infrared 1 st wave band and a green 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 main components of NDVI, NRI and NDSVI;
wherein D is a remote sensing image wave band set, i is the ith wave band of the predicted remote sensing image, a i Is a band linear combination coefficient, rho i The image surface reflectivity of the ith wave band is obtained;
wherein ,ρBlue In order to predict the blue light wave band of the remote sensing image, IBI is a built-up area index, and SI is a soil index;
wherein h is Planck constant, h =6.26 × 10 -34 J-sec; c is the speed of light, c =2.998 × 10 8 m/sec; k is Stefan-Boltzmann constant, K = 1.38X 10 -23 J/K; alpha is a middle variable, lambda is the central wavelength of a thermal infrared band, and epsilon is the ground surface emergence rate; t is b Is the brightness temperature;
(6) Denoising the remote sensing indexes VC, VHI, wet, NDBSI and LST by adopting a percentile denoising method, which specifically comprises the following steps: setting a percentile as a, sequencing pixel values contained in all research periods from small to large, taking a% as a lower limit of an effective value, taking (100-a)% as an upper limit of the effective value, and taking the lower limit of the effective value if the pixel values are less than the lower limit of the effective value; if the pixel value is larger than the upper limit of the effective value, the upper limit of the effective value is taken; then, normalizing the five remote sensing indexes through a relational expression (16) so as to shield the influence of the difference of numerical ranges of different remote sensing indexes;
wherein, index is the original remote sensing Index, index lb and Indexub Lower and upper limits, index, respectively, of the remote sensing Index nom The remote sensing index is normalized;
(7) And enabling each denoised and normalized remote sensing index to pass through minimum noise separation transformation (MNF), and taking the first component as a primary surface environment synthetic index LSECI 0 Observing each remote sensing finger in the first componentThe load value of the number is LSECI if the load value representing the vegetation coverage VC is more than or equal to zero 0 The final surface environment comprehensive index LSECI; if the load value representing the vegetation coverage VC is negative, LSECI is carried out 0 Obtaining a final surface environment comprehensive index LSECI by taking the inverse, wherein the relational expression 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 below.
Example one:
taking Hangzhou city as an experimental area and land resource satellite (Landsat) images as a data source, generating an annual (7 months and 1 day per year) earth surface environment comprehensive evaluation result, specifically:
(1) And acquiring an L2C2 Landsat image product published by USGS (United states geological survey) of the United states of America geological survey, wherein each image is systematically corrected in terms of geometry and radiation, and Landsat TM/ETM +/OLI image products in 1985-2020 are collected to obtain 4348 scenes in total. Each Landsat image used 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 value 3,4,5 represents the coverage of the surface with clouds, snow and cloud shadows, respectively. Based on the image, eliminating pixels covered by cloud, snow and cloud shadow in each scene image, and then sequentially arranging the images according to the sequence of imaging time to form an image time sequence;
(2) And for the optical wave band and the thermal infrared wave band, establishing an initial time sequence model of the pixels pixel by pixel:
firstly, 24 effective observation data are collected, an initial time sequence model is established by using a relational expression (3), and values of four parameters of a time sequence model are determined by adopting a least square method, which is shown in an attached figure 1;
(3) The change detection is carried out on the initial time sequence model of the pixel, and other 5 optical bands (green light, red light, near infrared, short wave infrared 1 and short wave infrared 2) are taken as model change detection bands in the example in consideration of the short wavelength of the blue light band and the high potential noise, so that the threshold value V of the change detection inv_x2 =15.1. Sequentially taking the subsequent 6 effective observation values as a group, carrying out change detection by using a relational expression (4), if the change is detected, taking the time of the 1 st effective observation value in the subsequent 6 effective observation values as the change occurrence time, taking the change occurrence time as a breakpoint in the time sequence model of the pixel, and segmenting the initial time sequence model to obtain a segmented time sequence model; otherwise, adding the 1 st of the subsequent 6 effective observed values into the initial time sequence model, fitting and updating the initial time sequence model of the pixel, and repeating the process until all the effective observed values are traversed, as shown in the attached figure 1;
(4) And obtaining a predicted remote sensing image at any moment: in the embodiment, the 7-month and 1-day of each year in 1986 to 2019 is assumed as the comprehensive evaluation time of the earth surface environment, the comprehensive evaluation time is converted into corresponding julian days, the time sequence model corresponding to the julian days is positioned on a wave band-by-pixel basis, the reflectivity of the earth surface of the pixel is obtained through a relational expression (1 to 3), on the basis, all the pixels form wave bands, all the wave bands form images, see the attached figure 1, and the predicted remote sensing image of the 7-month and 1-day of each year in 1986 to 2019 is obtained through the processes, namely 34 scenes in total;
(5) Calculating a remote sensing index: calculating vegetation coverage VC, vegetation health VHI, soil humidity Wet, surface civil engineering NDBSI and surface temperature LST one by one based on 34 scene prediction remote sensing images, wherein the relation of the soil humidity Wet suitable for the 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)
when calculating the surface temperature LST, using Global land30 for land cover products, see FIG. 1, determining the land use type of each pixel, and setting the corresponding surface emergence rate ε: 0.97-forest, 0.91-grassland, 0.92-artificial earth surface, 0.95-bare land, 0.99-water body, 0.97-other vegetation;
(6) Denoising and normalizing the remote sensing index: carrying out percentile denoising normalization on the images of all periods of each remote sensing index by using a relational expression (16), wherein 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% ], compared with the conventional maximum and minimum normalization method, the percentile denoising normalization method can better ensure the relative consistency of the value distribution before and after normalization, see the attached figure 2;
(7) And synthesizing a comprehensive index of the surface environment: for each remote sensing index of each period, performing minimum noise separation transformation by using a relational expression (17) to obtain a first component as a preliminary surface environment comprehensive index LSECI 0 And in the first component, because the load of the vegetation coverage VC is-0.632, the preliminary comprehensive index of the earth surface environment is inverted by using a relational expression (18) to obtain the final comprehensive index of the earth surface environment LSECI. The comprehensive evaluation result of the earth surface environment in 7/month and 1/year in 1986 to 2019 obtained by the above process is shown in figure 3, and the larger the numerical value is, the better the earth surface environment quality is represented.
Example two:
and taking a Hangzhou city as an experimental area and a land resource satellite (Landsat) image as a data source, and generating a ground surface environment comprehensive evaluation result of the semilunar period (10 days and 20 days per month) in 2019.
The remaining steps are the same as in example one, except: in the step (4), the comprehensive evaluation time of the earth surface environment of 10 days and 20 days per month in 2019 is assumed, the comprehensive evaluation time is converted into corresponding julian days, the time sequence model corresponding to the julian days is positioned pixel by waveband, the reflectivity of the earth surface of the pixels is obtained through the relational expression (1-3), then all the pixels form the waveband, all the wavebands form an image, and 24 scenes in total of the predicted remote sensing images of 10 days and 20 days per month in 2019 are obtained; the finally obtained 2019 semilunar surface environment comprehensive evaluation result is shown in fig. 4, which further proves that the comprehensive evaluation method can be used for surface environment comprehensive evaluation at any time.
The scope of the present invention includes, but is not limited to, the above embodiments, and any alternatives, modifications, and improvements that may be easily made by those skilled in the art are within the scope of the present invention as defined by the appended claims.
Claims (8)
1. A comprehensive evaluation method for any time earth surface environment 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 sequentially 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;
(3) Detecting the change of the initial time sequence model of the pixel, and if the change does not occur, keeping the initial time sequence model; if the change occurs, the time of the change occurrence is taken as a breakpoint of the time sequence model, and the initial time sequence model of the pixel is segmented 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 the time sequence model of the pixel;
(5) Based on the predicted remote sensing image, taking five remote sensing indexes of vegetation coverage VC, vegetation health VHI, soil humidity Wet, surface civil engineering NDBSI and surface temperature LST as surface environment comprehensive representation of the selected area, and calculating numerical values of the five remote sensing indexes;
(6) Denoising and normalizing the remote sensing index;
(7) And obtaining the comprehensive index of the surface environment.
2. The method for comprehensively evaluating the earth's surface environment at any moment based on time series remote sensing as claimed in claim 1, wherein the specific method for preprocessing the real-shot remote sensing image in the step (1) is as follows: downloading all real-shot optical remote sensing image data which are subjected to geometric correction and radiation correction in a selected area in a research period, or downloading all original real-shot optical remote sensing image data, then manually carrying out geometric correction and radiation correction, and then utilizing the quality evaluation wave band of the remote sensing image to remove invalid observation values covered by clouds, cloud shadows or snow in each image.
3. The method for comprehensively evaluating the earth's surface environment at any moment based on time series remote sensing as claimed in claim 1, wherein the specific establishment method of the step (2) is as follows: the initial time sequence model of each pixel consists of a linear function and a harmonic function, three relations are constructed from coarse to fine according to the number of effective observed values of each pixel in a time sequence, and if the number of the effective observed values is not less than 12 and is less than 18, the model selects the relation (1); if the number of the effective observed values is more than or equal to 18 and less than 24, selecting a relational expression (2) for the model; if the number of the effective observed values is more than or equal to 24, selecting a relational expression (3) for the model;
wherein ,
representing the 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 the julian day; w is the frequency of the light emitted by the light source,
c 0,i and c1,i As 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 are determined to be optimal values by a least square method.
4. The method for comprehensively evaluating the earth's surface environment at any moment based on time series remote sensing as claimed in claim 3, wherein the specific detection method in step (3) is as follows: 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 error resid of the initial time series model predicted value and the subsequent effective observed value of the pixel with the root mean square error RMSE of the initial time series model, i.e.
Wherein: d is the wave band set of the remote sensing image, i is the ith wave band of the remote sensing image, reset i Is the residual error, RMSE, of the initial time series model predicted value and the subsequent effective observed value of the ith wave band i Is the root mean square error of the initial time series model of the picture elements of the ith band,
is a threshold value, and is equal to the inverse Chi-square distribution number of the degree of freedom which is the number of image wave bands;
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 the wave bands is larger than a threshold value, taking the time of the occurrence of the 1 st effective observation value of the pixel in the subsequent 6 effective observation values of each wave band as the change occurrence time, taking the change occurrence time as a breakpoint of the pixel in the time sequence model of each wave band, segmenting the initial time sequence model of the pixel in each wave band, and obtaining a segmented time sequence model according to the relational expression in the step (2) to reflect the surface reflectivity change of the pixel;
if the ratio of the resid to the RMSE of the 6 subsequent effective observation values of the same pixel in all the wave bands is not all larger than the threshold value, the change is considered to be not generated, the 1 st effective observation value of the pixel in the 6 subsequent effective observation values in each wave band is added into the initial time sequence model of the corresponding wave band, and then the initial time sequence model of the pixel is fitted and updated according to the relational expression in the step (2);
the above-mentioned discrimination process is repeated until all the effective observed values are located in the initial time series model or the segmented time series model of the pixel.
5. The method for comprehensively evaluating the earth's surface environment at any moment based on time series remote sensing as claimed in claim 3, wherein 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 on a wave band-by-wave pixel-by-pixel basis, determining the julian day t of the comprehensive evaluation of the earth surface environment according to the research requirement, obtaining the predicted earth surface reflectivity of the pixels through the relational expression in the step (2), then forming all the pixels into wave bands, forming all the wave bands into images, and obtaining the predicted remote sensing image at any moment.
6. The method for comprehensively evaluating the earth's surface environment at any moment based on time series remote sensing as claimed in claim 1, wherein the five remote sensing indexes in the step (5) are as follows: the relation formula of the vegetation coverage VC, the vegetation health VHI, the soil humidity Wet, the surface civil engineering NDBSI and the surface temperature LST is as follows:
wherein ,ρNIR and ρRed Respectively near infrared band and red light band of the predicted remote sensing image, and NDVI is Normalized Difference Vegetation Index (NDVI) veg and NDVIsoil The normalized difference vegetation indexes of the vegetation and the soil are respectively the maximum value and the minimum value of the NDVI in the research range;
VHI＝PCA1(NDVI,NRI,NDSVI) (9)
wherein ,ρSWIR1 and ρGreen Respectively predicting a short wave infrared 1 st wave band and a green 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 main components of NDVI, NRI and NDSVI;
wherein D is a remote sensing image wave band set, i is the ith wave band of the predicted remote sensing image, a i Is a band linear combination coefficient, rho i The image surface reflectivity of the ith wave band is obtained;
wherein ,ρBlue In order to predict the blue light wave band of the remote sensing image, IBI is a built-up area index, and SI is a soil index;
wherein h is Planck constant, h =6.26 × 10 -34 J-sec; c is the speed of light, c =2.998 × 10 8 m/sec; k is Stefan-Boltzmann constant, K = 1.38X 10 -23 J/K; alpha is a middle variable, lambda is the central wavelength of a thermal infrared band, and epsilon is the ground surface emergence rate; t is a unit of b Is the brightness temperature.
7. The method for comprehensively evaluating the earth's surface environment at any moment based on time series remote sensing as claimed in claim 1, wherein the specific method in step (6) is as follows:
adopting a percentile denoising method for the remote sensing indexes VC, VHI, wet, NDBSI and LST, setting the percentile as a, sequencing pixel values contained in all research periods from small to large, taking a% as the lower limit of an effective value, taking (100-a)% as the upper limit of the effective value, and taking the lower limit of the effective value if the pixel values are smaller than the lower limit of the effective value; if the pixel value is larger than the upper limit of the effective value, the upper limit of the effective value is taken; then, normalizing the five remote sensing indexes through a relational expression (16) so as to shield the influence of the difference of the numerical ranges of different remote sensing indexes;
wherein Index is the original remoteSensory Index, index lb and Indexub Lower and upper limits, index, respectively, of the remote sensing Index nom The index is normalized remote sensing index.
8. The method for comprehensively evaluating the earth's surface environment at any moment based on time series remote sensing as claimed in claim 1, wherein the specific method of the step (7) is as follows:
carrying out minimum noise separation transformation on each denoised and normalized remote sensing index, and taking the first component as a primary 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 more than or equal to zero, then LSECI 0 The final surface environment comprehensive index LSECI; if the load value representing the vegetation coverage VC is negative, LSECI is carried out 0 Obtaining a final surface environment comprehensive index LSECI by taking the inverse, wherein the relational expression 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 true CN115203624A (en) | 2022-10-18 |
| CN115203624B 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) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116466368A (en) * | 2023-06-16 | 2023-07-21 | 成都远望科技有限责任公司 | Dust extinction coefficient profile estimation method based on laser radar and satellite data |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20200311842A1 (en) * | 2019-03-28 | 2020-10-01 | China Waterborne Transport Research Institute | Method for tracking, monitoring and evaluating ecological impact of channel project based on long-term time series satellite remote sensing data |
| US20210027429A1 (en) * | 2019-07-26 | 2021-01-28 | Zhejiang University Of Technology | Noise detection method for time-series vegetation index derived from remote sensing images |
| CN113988626A (en) * | 2021-10-28 | 2022-01-28 | 中国矿业大学(北京) | 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
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20200311842A1 (en) * | 2019-03-28 | 2020-10-01 | China Waterborne Transport Research Institute | Method for tracking, monitoring and evaluating ecological impact of channel project based on long-term time series satellite remote sensing data |
| US20210027429A1 (en) * | 2019-07-26 | 2021-01-28 | Zhejiang University Of Technology | Noise detection method for time-series vegetation index derived from remote sensing images |
| CN113988626A (en) * | 2021-10-28 | 2022-01-28 | 中国矿业大学(北京) | 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 |
Non-Patent Citations (1)
| Title |
|---|
| 李红星 等: "基于遥感生态指数的武汉市生态环境质量评估", 云南大学学报(自然科学版) * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116466368A (en) * | 2023-06-16 | 2023-07-21 | 成都远望科技有限责任公司 | Dust extinction coefficient profile estimation method based on laser radar and satellite data |
| CN116466368B (en) * | 2023-06-16 | 2023-08-22 | 成都远望科技有限责任公司 | Dust extinction coefficient profile estimation method based on laser radar and satellite data |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115203624B (en) | 2023-06-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Zarco‐Tejada et al. | Land cover mapping at BOREAS using red edge spectral parameters from CASI imagery | |
| Song et al. | Predicting temperate conifer forest successional stage distributions with multitemporal Landsat Thematic Mapper imagery | |
| CN112395808A (en) | Biomass remote sensing mapping method combining random forest and collaborative kriging | |
| CN110927120B (en) | Early warning method for vegetation coverage | |
| CN112733596A (en) | Forest resource change monitoring method based on medium and high spatial resolution remote sensing image fusion and application | |
| CN114821349B (en) | Forest biomass estimation method taking harmonic model coefficients and climatic parameters into consideration | |
| CN115372282B (en) | Farmland soil water content monitoring method based on hyperspectral image of unmanned aerial vehicle | |
| CN114241331B (en) | Remote sensing modeling method for ground biomass of reed in wetland by taking UAV as ground and Septinel-2 medium | |
| CN114819737B (en) | Method, system and storage medium for estimating carbon reserves of highway road vegetation | |
| Fitzgerald et al. | Directed sampling using remote sensing with a response surface sampling design for site-specific agriculture | |
| CN115203624B (en) | Comprehensive evaluation method for surface environment at any moment based on time sequence remote sensing | |
| CN112215135B (en) | Mining area mining and treatment effect monitoring method and device | |
| Khudhur et al. | Comparison of the accuracies of different spectral indices for mapping the vegetation covers in Al-Hawija district, Iraq | |
| Singh et al. | Environmental degradation analysis using NOAA/AVHRR data | |
| Das et al. | Mapping vegetation and forest types using Landsat TM in the western ghat region of Maharashtra, India | |
| Fernandes et al. | A multi-scale approach to mapping effective leaf area index in boreal Picea mariana stands using high spatial resolution CASI imagery | |
| Taha | Assessment of urbanization encroachment over Al-Monib island using fuzzy post classification comparison and urbanization metrics | |
| Hakkenberg et al. | Automated continuous fields prediction from Landsat time series: application to fractional impervious cover | |
| Nielsen et al. | The distribution in time and space of savanna fires in Burkina Faso as determined from NOAA AVHRR data | |
| Colaninno et al. | Defining Densities for Urban Residential Texture, through Land use Classification, from Landsat TM Imagery: Case Study of Spanish Mediterranean Coast | |
| CN117196922A (en) | Remote sensing method for continuously evaluating ecological quality of natural protection area | |
| Zhao et al. | Towards an Annually Updated Oil Palm Age Database (OPAD) with Multimodal Remote Sensing Features: Design and Preliminary Result | |
| Gupta et al. | Fraction of vegetation cover and its application in vegetation characterization in the Hazaribagh wildlife sanctuary, Jharkhand, India | |
| Singh et al. | Estimation of changes in land surface temperature using multi-temporal Landsat data of Ghaziabad District, India. | |
| Wang | Monitoring time-series changes in vegetation coverage based on landsat images |
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 |