CN116862798A - Topographic shadow effect correction method for optical satellite remote sensing image - Google Patents
Topographic shadow effect correction method for optical satellite remote sensing image Download PDFInfo
- Publication number
- CN116862798A CN116862798A CN202310831948.0A CN202310831948A CN116862798A CN 116862798 A CN116862798 A CN 116862798A CN 202310831948 A CN202310831948 A CN 202310831948A CN 116862798 A CN116862798 A CN 116862798A
- Authority
- CN
- China
- Prior art keywords
- correction
- factor
- toa
- baf
- reflectivity
- 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
- 238000012937 correction Methods 0.000 title claims abstract description 75
- 238000000034 method Methods 0.000 title claims abstract description 41
- 230000000694 effects Effects 0.000 title claims abstract description 38
- 230000003287 optical effect Effects 0.000 title claims abstract description 15
- 238000002310 reflectometry Methods 0.000 claims abstract description 70
- 230000005855 radiation Effects 0.000 claims abstract description 40
- 230000003595 spectral effect Effects 0.000 claims abstract description 14
- 238000012876 topography Methods 0.000 claims abstract description 14
- 238000005286 illumination Methods 0.000 claims abstract description 12
- 238000004451 qualitative analysis Methods 0.000 claims abstract description 7
- 238000004445 quantitative analysis Methods 0.000 claims abstract description 7
- 238000004458 analytical method Methods 0.000 claims description 13
- 108010014172 Factor V Proteins 0.000 claims description 6
- 238000011179 visual inspection Methods 0.000 claims description 6
- 238000013459 approach Methods 0.000 claims description 4
- 238000011160 research Methods 0.000 claims description 4
- 239000000443 aerosol Substances 0.000 claims description 3
- 238000003702 image correction Methods 0.000 claims description 3
- 238000003384 imaging method Methods 0.000 claims description 3
- 238000001579 optical reflectometry Methods 0.000 claims description 3
- 238000002360 preparation method Methods 0.000 claims description 3
- 230000008569 process Effects 0.000 abstract description 3
- 230000008859 change Effects 0.000 description 4
- 238000001228 spectrum Methods 0.000 description 4
- 230000007547 defect Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 244000208060 Lawsonia inermis Species 0.000 description 1
- 230000002159 abnormal effect Effects 0.000 description 1
- 230000006978 adaptation Effects 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 239000002689 soil Substances 0.000 description 1
- 230000001360 synchronised effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/05—Geographic models
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
- G06T2207/10036—Multispectral image; Hyperspectral image
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30181—Earth observation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30181—Earth observation
- G06T2207/30188—Vegetation; Agriculture
-
- 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)
- Geometry (AREA)
- Software Systems (AREA)
- Remote Sensing (AREA)
- Computer Graphics (AREA)
- Image Processing (AREA)
Abstract
The application discloses a topographic shadow effect correction method of an optical satellite remote sensing image, which aims to solve the problem of surface reflectivity distortion of a shadow area in the remote sensing image caused by the influence of topography. The method is based on the earth surface reflectivity of the pixels and solar radiation received by the earth surface of the illumination area and the shadow area, introduces Shadow Intensity Factors (SIFs), vegetation Index Factors (VIFs) and wave Band Adjustment Factors (BAFs) to construct a terrain shadow effect positive model, and comprehensively evaluates the correction performance of the model through qualitative and quantitative analysis. According to the application, the spectral distortion of the shadow area can be effectively corrected, the original spectral characteristic of the illumination area is not interfered, good decorrelation performance is shown between the surface reflectivity after correction and the illumination condition, the NDVI has high stability after correction, vegetation information is correctly reserved, the phenomena of overcorrection of the shadow area and undercorrection of the shadow area can be restrained, and other data are not required to be introduced in the correction process except DEM.
Description
Technical Field
The application belongs to the technical field of optical satellite remote sensing image preprocessing, and particularly relates to a topographic shadow effect correction method of an optical satellite remote sensing image.
Background
The surface reflectance is a physical quantity that characterizes the solar radiation reflecting capacity and is an inherent property of the ground feature. The method for accurately estimating the surface reflectivity has important significance for quantitative remote sensing theory foundation construction, remote sensing product application and the like. However, in mountainous areas with large relief, the radiation transmission process of solar radiation between the sun-earth surface-sensor is changed, and the effective radiation received by the earth surface pixels is different, so that the same earth surface coverage receives weaker radiation on a hillside and has a lower reflectivity value. This phenomenon is known as the topographic shadow effect. The existence of the shadow effect of the terrain seriously influences the applications of fine classification of the surface coverage, inversion of vegetation parameters, quantitative remote sensing and the like. Therefore, correcting or reducing the effect of terrain shadows is an important step in mountain remote sensing pretreatment.
Currently, correction methods for topographical shading effects can be classified into 3 categories: (1) shadow detection and shadow information compensation method. Such methods are not limited to a certain type of ground object, but rather detect and compensate for information of shadow areas through brightness information (such as DN values). However, the shadow compensation effect is different from one feature to another, and in addition, such methods generally require a large amount of calculation and are rarely applied to mountain regions. (2) recovering information of the shadow region based on the multi-source data. The obvious defects of the method are that the differential assimilation process of the heterologous data is complex, the acquisition difficulty of the multi-source synchronous data is high, and the problems limit the universality of the method. (3) terrain correction modeling method. The method corrects the surface reflectivity by constructing a physical model, a semi-empirical model or an empirical model to eliminate or reduce the effect of the shadow of the terrain. This method is generally simple and easy to implement, especially empirical and semi-empirical models, and therefore, such methods are widely used. However, this approach has 2 drawbacks: 1) In the area with lower solar altitude, the correction effect is poor; 2) The ghost area is liable to cause overcorrection and the drop area is liable to cause undercorrection.
In summary, the topographic shadow effect correction method for the optical satellite remote sensing image has important scientific significance and practical value in improving the correction precision of the earth surface reflectivity of the shadow area of the remote sensing image and recovering the real spectrum information of the remote sensing image in the mountain area.
Disclosure of Invention
The application aims to solve the technical problems and overcome the defects of the prior art, and provides a terrain shading effect correction method for an optical satellite remote sensing image, which adopts the following basic conception:
a topographic shading effect correction method of an optical satellite remote sensing image comprises the following steps:
step S1, data preparation: acquiring multispectral satellite image data (apparent reflectivity image data and earth surface reflectivity image data) of a mountain area and corresponding Digital Elevation Model (DEM) data;
step S2, obtaining parameters in the model: calculating direct solar radiation, sky scattered radiation and topography reflected radiation by using a 6S model through image imaging time, solar zenith angle and azimuth angle, atmospheric mode, aerosol type, visibility, sensor height, spectral parameters and surface reflectivity type parameters, and calculating sky observation factor V by using a DEM d And topography observation factor C t ；
Step S3, shadow intensity factor SIF calculation: calculating an apparent reflectivity shadow index to realize the representation of the shadow intensity in the image, and establishing a shadow intensity factor by utilizing the shadow index in combination with a threshold value;
step S4, calculating a vegetation index factor VIF: NDVI is taken as a vegetation index factor VIF to reflect the condition of the shadow coverage area;
step S5, rough calculation of band adjustment factor BAF: selecting a proper amount of samples with the same quality of the sunny slope and the cloudy slope, taking the BAF value from 0 as a step length, and determining a rough BAF value according to the relation between the average values of the reflectivities of the sunny slope and the cloudy slope samples under different BAF values;
step S6, constructing a terrain shading effect correction model: constructing a topographic shading effect correction model according to the surface reflectivity image data, the direct solar radiation, the sky scattered radiation, the topographic reflected radiation, the sky observation factor, the topographic observation factor, the SIF factor, the VIF factor and the BAF factor to calculate corrected surface reflectivity;
step S7, accurate calculation of a band adjustment factor BAF: on the basis of the step S5, determining an optimal BAF value through the relation between the sunny slope and cloudy slope samples and the cosine cosi of the sun incidence angle under different BAF values, thereby adjusting the accurate BAF value obtained in the step S5;
step S8, evaluating model performance: the performance of the model was assessed comprehensively by qualitative and quantitative analysis.
Further, in the step S2, the sky observation factor V d And topography observation factor C t Calculated according to the following formula:
where S is the slope value of each pel calculated from the DEM data.
Further, the step S3 specifically includes:
step S31: the apparent reflectance shading index is calculated according to the following formula,
in SI TOA Is the apparent reflectance shading index, ρ C_TOA 、ρ G_TOA And ρ NIR_TOA Apparent reflectivities of a Coastal band, a Green band and an NIR band respectively;
step S32: with SI at peak TOA Value as threshold value of earth surface reflectivity correction, SI TOA Greater than the peakThe pixel of the value needs to be corrected, otherwise, the correction is not needed, so that the SIF factor is constructed, and the specific calculation formula of the SIF factor is as follows:
in SI TOA_p Is SI TOA SI at histogram peak TOA Value, SI TOA_max Is SI TOA Is a maximum value of (a).
Further, the VIF factor in step S4 is calculated according to the following formula:
wherein ρ is NIR And ρ R The surface reflectivities of the NIR band and the Red band, respectively.
Further, in the step S5, when the average values of the reflectivities of the sunny slope and the cloudy slope are equal, a rough BAF value is obtained;
further, the corrected surface reflectivity in the step S6 is calculated according to the following formula:
wherein ρ is corr Is the corrected pixel reflectivity ρ orig Is the pixel reflectivity before correction, f (SIF, VIF, BAF) is the self-adjusting direct light reflectivity compensation intensity function calculated by multiplying the values of three factors equal to SIF, VIF, BAF, E d Is the direct solar radiation, E, in said step S2 f Is the solar scattered radiation in said step S2, E is the sum of the direct solar radiation, the sky scattered radiation, the topography reflected radiation in said step S2, ρ a Is the average reflectivity of neighboring pixels of each pixel.
Further, the step S7 specifically includes:
step S71: respectively selecting vegetation sample points with the same quality as a shade slope and a sunny slope in a research area;
step S72: calculating the correction result of each sample point under each BAF value from 0 by taking 0.05 as step length, and calculating the correlation coefficient r of the sample reflectivity of the sunny slope before correction and cosi 1 ；
Step S73: corrected SI TOA_p +2 times SI TOA Standard deviation of values (noted as SI TOA_p +2σ) correlation coefficient r of reflectance and cosi of all samples of the left region 2 ；
Step S74: when r is 1 And r 2 When the difference value of (2) approaches 0, the BAF value is the optimal value;
wherein, in the above steps, cosi is calculated according to the following formula:
wherein Z is the zenith angle of the sun, S is the gradient,is the solar azimuth angle and β is the slope direction.
Further, in the step S8, the qualitative analysis is a visual inspection analysis, and the quantitative analysis includes: and comparing the reflectivity of evergreen broad-leaved forests on the hillside and the sunny hillside, performing decorrelation analysis between the surface reflectivity and illumination conditions, performing NDVI correction stability analysis, and performing penumbra and falling shadow correction performance analysis.
By adopting the technical scheme, compared with the prior art, the application has the following beneficial effects.
According to the application, a shadow effect correction model is built based on solar radiation difference between an illumination area and a shadow area, a shadow intensity factor, a vegetation index factor and a wave band adjustment factor are introduced, so that correction of reflectivity of the shadow area is realized, spectral distortion of the shadow area can be effectively corrected without disturbing original spectral characteristics of the illumination area, good decorrelation performance is shown between the surface reflectivity after correction and illumination conditions, NDVI (non-linear density vi) correction has high stability, vegetation information is correctly reserved, phenomena of overcorrection and insufficient correction of the shadow area can be restrained, and other data are not required to be introduced in the correction process except DEM. The present application provides a simple and effective method to correct the effects of terrain shadows in mountain areas with satisfactory results.
The following describes the embodiments of the present application in further detail with reference to the accompanying drawings.
Drawings
The accompanying drawings, which are included to provide a further understanding of the application and are incorporated in and constitute a part of this specification, illustrate embodiments of the application and together with the description serve to explain the application. It is evident that the drawings in the following description are only examples, from which other drawings can be obtained by a person skilled in the art without the inventive effort. In the drawings:
FIG. 1 is a schematic flow chart of the present application;
FIG. 2 is a visual inspection image (mountain area) before and after correction according to an embodiment of the present application;
FIG. 3 is a visual inspection image (Changjiang area) before and after correction according to an embodiment of the present application;
FIG. 4 is a schematic diagram showing the calculation of SIF factors according to an embodiment of the present application;
FIG. 5 is a graph showing spectra of a shade slope and an adjacent shade slope before and after correction according to an embodiment of the present application;
FIG. 6 is a graph of density scattergrams of shadow index values and cosine values of the angle of incidence of the sun (cosi) and its determination coefficients (R 2 )；
FIG. 7 is a graph showing the NDVI difference distribution before and after correction in accordance with an embodiment of the present application.
It should be noted that these drawings and the written description are not intended to limit the scope of the inventive concept in any way, but to illustrate the inventive concept to those skilled in the art by referring to the specific embodiments.
Detailed Description
For the purpose of making the objects, technical solutions and advantages of the embodiments of the present application more apparent, the technical solutions in the embodiments will be clearly and completely described with reference to the accompanying drawings in the embodiments of the present application, and the following embodiments are used to illustrate the present application, but are not intended to limit the scope of the present application.
Example 1
As shown in fig. 1, the method for correcting the topographic shading effect of the optical satellite remote sensing image according to the embodiment includes the following steps:
step S1, data preparation: multispectral satellite image data (apparent reflectivity image data and surface reflectivity image data) of the mountain region and corresponding digital elevation model DEM data are acquired.
Specifically, landsat8OLI multispectral images of different time periods (3 months, 6 months, 9 months and 12 months) in the wuyishan region of the foggy province and the changjiang region of the henna province were acquired, and since the mountain area was always covered with clouds, an almost cloud-free image nearest to these four time periods was selected as the image data of the present embodiment, as shown in fig. 2 and 3. The image blocks (a) - (d) in fig. 2 are pseudo-color images before correction in the wuyishan region, and it can be seen that shadows in the images are obvious before correction, especially the shadow surfaces of mountains with large elevation drop, and pixels in the shadow regions are black. Furthermore, there is a significant difference in reflectivity between the hillside and the sunny side. The blocks (a) - (d) in fig. 3 are pseudo-color images before correction in the changjiang area, and it can be seen that the area is located in the tropical area with low latitude, and the shadow effect is not serious throughout the year, especially the images of 5 months and 9 months, which helps to evaluate the stability of the method proposed by the present application.
Step S2, obtaining parameters in the model: calculating direct solar radiation, sky scattered radiation and topography reflected radiation by using a 6S model through image imaging time, solar zenith angle and azimuth angle, atmospheric mode, aerosol type, visibility, sensor height, spectral parameters and surface reflectivity type parameters, and calculating sky observation factor V by using a DEM d And topography observation factor C t 。
Sky observation factor V d And topography observation factor C t Calculated according to the following formula:
where S is the slope value of each pel calculated from the DEM data.
Step S3, shadow intensity factor SIF calculation: and calculating an apparent reflectivity shadow index to realize the representation of the shadow intensity in the image, and establishing a shadow intensity factor by utilizing the shadow index in combination with a threshold value.
Step S3 here comprises the steps of:
step S31: the apparent reflectance shading index is calculated according to the following formula,
in SI TOA Is the apparent reflectance shading index, ρ C_TOA 、ρ G_TOA And ρ NIR_TOA Apparent reflectivities of a Coastal band, a Green band and an NIR band respectively;
step S32: with SI at peak TOA Value as threshold value of earth surface reflectivity correction, SI TOA The pixels larger than the peak need correction, otherwise, do not need correction (for the schematic diagram, see fig. 4), so as to construct the SIF factor, and the specific calculation formula of the SIF factor is as follows:
in SI TOA_p Is SI TOA SI at histogram peak TOA Value, SI TOA_max Is SI TOA Is a maximum value of (a).
Step S4, calculating a vegetation index factor VIF: NDVI is used as a vegetation index factor VIF to reflect the shadow coverage situation.
Wherein the VIF factor is calculated according to the following formula:
wherein ρ is NIR And ρ R The surface reflectivities of the NIR band and the Red band, respectively.
Step S5, rough calculation of band adjustment factor BAF: and selecting a proper amount of samples with the same quality of the sunny slope and the cloudy slope, taking the BAF value from 0 as a step length, and determining a rough BAF value through the relation between the average values of the reflectivities of the sunny slope and the cloudy slope samples under different BAF values. When the average reflectivity values of the sunny and cloudy slope samples are equal, a rough BAF value is obtained.
Step S6, constructing a terrain shading effect correction model: and constructing a topographic shading effect correction model by using the surface reflectivity image data, the direct solar radiation, the sky scattered radiation, the topographic reflected radiation, the sky observation factor, the topographic observation factor, the SIF factor, the VIF factor and the BAF factor to calculate corrected surface reflectivity.
The corrected surface reflectivity is calculated according to the following formula:
wherein ρ is corr Is the corrected pixel reflectivity ρ orig Is the pixel reflectivity before correction, f (SIF, VIF, BAF) is the self-adjusting direct light reflectivity compensation intensity function calculated by multiplying the values of three factors equal to SIF, VIF, BAF, E d Is the direct solar radiation, E, in said step S2 f Is the solar scattered radiation in said step S2, E is the sum of the direct solar radiation, the sky scattered radiation, the topography reflected radiation in said step S2, ρ a Is the average reflectivity of neighboring pixels of each pixel.
Step S7, accurate calculation of a band adjustment factor BAF: on the basis of step S5, the optimal BAF value is determined by the relation between the sun slope, the shade slope sample and the sun angle of incidence cosine cosi under different BAF values, so that the exact BAF value obtained in step S5 is adjusted.
The method specifically comprises the following steps:
step S71: respectively selecting vegetation sample points with the same quality as a shade slope and a sunny slope in a research area;
step S72: calculating the correction result of each sample point under each BAF value from 0 by taking 0.05 as step length, and calculating the correlation coefficient r of the sample reflectivity of the sunny slope before correction and cosi 1 ；
Step S73: corrected SI TOA_p +2 times SI TOA Standard deviation of values (noted as SI TOA_p +2σ) correlation coefficient r of reflectance and cosi of all samples of the left region 2 ；
Step S74: when r is 1 And r 2 When the difference value of (2) approaches 0, the BAF value is the optimal value;
wherein, in the above steps, cosi is calculated according to the following formula:
wherein Z is the zenith angle of the sun, S is the gradient,is the solar azimuth angle and β is the slope direction.
Step S8, evaluating model performance: the performance of the model was assessed comprehensively by qualitative and quantitative analysis.
The qualitative analysis is a visual inspection analysis, the quantitative analysis comprising: and comparing the reflectivity of evergreen broad-leaved forests on the hillside and the sunny hillside, performing decorrelation analysis between the surface reflectivity and illumination conditions, performing NDVI correction stability analysis, and performing penumbra and falling shadow correction performance analysis.
Specific visual inspection results are shown in blocks (e) to (h) in fig. 2 and blocks (e) to (h) in fig. 3. The blocks (e) - (h) in fig. 2 are pseudo-color images corrected in the wuyishan region, and it can be obviously seen from the figure that the difference between the corrected sunny slope and the negative slope is obviously reduced, the whole image is smoother, and the phenomena of insufficient correction and excessive correction are not obvious. The spectral characteristics of the sparse grasslands on the mountain tops in the north of the area are greatly influenced by the change of the climates, the spectral characteristics of vegetation are shown in 5-month and 9-month images, the spectral characteristics of bare soil are shown in 3 months and 12 months, and the change of the spectral characteristics is reserved in the corrected images. The blocks (e) - (h) in fig. 3 are respectively pseudo-color images corrected by the changjiang area, and the color tone corrected by 4 periods is very similar because the change of the vegetation weathers of the research area is relatively small and the time interval for image acquisition is very short and the change of the land coverage type is very small. The cultivated land in the southeast of the area has no obvious chromatic aberration before and after correction, which indicates that the correction performance of the TSEC method on each wave band is basically the same.
To further quantitatively evaluate the performance of the method of the present application, a sub-region (white rectangular box in the (a) plot of fig. 3) was selected in the changjiang area for analysis of the corrected positive and negative slope reflectivities. In the sub-areas 8 distinct evergreen broad-leaved forest covered hilly areas are marked (see (a) block in fig. 3), and around these areas four sunny slopes are marked, all of which can be illuminated by the sun. Fig. 5 shows the spectral curves plotted for 635 sunny slope pixels and 987 cloudy slope pixels marked. It can be seen that the negative slope reflectivity in the images of 2 months and 12 months is obviously improved, only the NIR wave band in the images of 2 months is slightly overcorrected, and the correction effect of other wave bands is good. In addition, the 5 month image has a high solar altitude, so that there is almost no shadow area in the image, and therefore, the spectral curves before and after correction almost completely overlap. In general, the spectrum curve of the shade slope and the spectrum curve of the sun slope after the image correction of the 4 time periods are high in matching degree.
FIG. 6 is a plot of the density scatter between the SWIR1 band reflectance and cosi before and after image correction in the Wuyi mountain area. It is apparent that the presence of shadows in the image reduces the reflectance of the band before correction with reduced illumination intensity, especially for 12 months of images, which tends to be most pronounced. After correction, the density points between the reflectivity and cosi show a horizontal distribution trend, which shows that the degree of decorrelation between the reflectivity of the corrected wave band and cosi is high, and R is the same for 4 time periods after correction 2 All lower than 0.06, indicating that the reflectivity of the picture element is hardly affected by the intensity of illumination.
Fig. 7 is a graph of NDVI difference profiles before and after correction for all periods in the mountain area of wuyie. Since NDVI has a certain capability of eliminating shadow effect, for an image in which shadow effect is not serious, the difference between NDVI before and after correction should be small, close to 0. Except for 12 months of images in the Wuyi mountain study area, the other 7 images meet the theory, and the NDVI difference before and after correction is within 0.07. Particularly 5 months and 9 months. For images on the severe shadow effect date (12 months), the difference between NDVI before and after correction is large, even if some areas exceed 0.15, and the number of pixels exceeding 0.05 is also large, which indicates that the pixels in the severe shadow effect area are effectively corrected. Therefore, the method can effectively correct the region with serious shadow effect, hardly interfere with the illumination region and completely retain the original spectral characteristics. In addition, the corrected NDVI has no obvious distortion (namely serious overcorrection or undercorrection) phenomenon, and original vegetation information is reserved.
Table 1 shows the Mean (Mean), standard deviation (Std), coefficient of Variation (CV) and Outliers (Outliers) of NIR bands in images of Wuyi mountain areas before and after correction of the home and fall areas, and compares them with 2 conventional topography correction methods.
TABLE 1 mean, standard deviation, CV and outliers of NIR band reflectivities before and after correction for the ghost and ghost areas
The mean value, standard deviation and variation coefficient of the original shadow area and the shadow falling area after correction by the TSEC method are very close and superior to those of the traditional terrain correction method. The method can well inhibit the overcorrection of the shadow area and improve the undercorrection phenomenon of the shadow area, has high stability for correcting the shadow area, and hardly generates any abnormal value.
In summary, the application can simply and efficiently correct the topographic shadow effect of the satellite remote sensing image.
The foregoing description is only illustrative of the preferred embodiment of the present application, and is not to be construed as limiting the application, but is to be construed as limiting the application to any and all simple modifications, equivalent variations and adaptations of the embodiments described above, which are within the scope of the application, may be made by those skilled in the art without departing from the scope of the application.
Claims (8)
1. A topographic shading effect correction method of an optical satellite remote sensing image is characterized by comprising the following steps:
step S1, data preparation: acquiring multispectral satellite image data (apparent reflectivity image data and earth surface reflectivity image data) of a mountain area and corresponding Digital Elevation Model (DEM) data;
step S2, obtaining parameters in the model: calculating direct solar radiation, sky scattered radiation and topography reflected radiation by using a 6S model through image imaging time, solar zenith angle and azimuth angle, atmospheric mode, aerosol type, visibility, sensor height, spectral parameters and surface reflectivity type parameters, and calculating sky observation factor V by using a DEM d And topography observation factor C t ；
Step S3, shadow intensity factor SIF calculation: calculating an apparent reflectivity shadow index to realize the representation of the shadow intensity in the image, and establishing a shadow intensity factor by utilizing the shadow index in combination with a threshold value;
step S4, calculating a vegetation index factor VIF: NDVI is taken as a vegetation index factor VIF to reflect the condition of the shadow coverage area;
step S5, rough calculation of band adjustment factor BAF: selecting a proper amount of samples with the same quality of the sunny slope and the cloudy slope, taking the BAF value from 0 as a step length, and determining a rough BAF value according to the relation between the average values of the reflectivities of the sunny slope and the cloudy slope samples under different BAF values;
step S6, constructing a terrain shading effect correction model: constructing a topographic shading effect correction model according to the surface reflectivity image data, the direct solar radiation, the sky scattered radiation, the topographic reflected radiation, the sky observation factor, the topographic observation factor, the SIF factor, the VIF factor and the BAF factor to calculate corrected surface reflectivity;
step S7, accurate calculation of a band adjustment factor BAF: on the basis of the step S5, determining an optimal BAF value through the relation between the sunny slope and cloudy slope samples and the cosine cosi of the sun incidence angle under different BAF values, thereby adjusting the accurate BAF value obtained in the step S5;
step S8, evaluating model performance: the performance of the model was assessed comprehensively by qualitative and quantitative analysis.
2. The method for correcting topographic shading effects of an optical satellite remote sensing image according to claim 1, wherein in the step S2, the sky observation factor V d And topography observation factor C t Calculated according to the following formula:
where S is the slope value of each pel calculated from the DEM data.
3. The method for correcting the topographic shading effect of the optical satellite remote sensing image according to claim 1, wherein the step S3 is specifically:
step S31: the apparent reflectance shading index is calculated according to the following formula,
in SI TOA Is the apparent reflectance shading index, ρ C_TOA 、ρ G_TOA And ρ NIR_TOA Apparent reflectivities of a Coastal band, a Green band and an NIR band respectively;
step S32: with SI at peak TOA Value as threshold value of earth surface reflectivity correction, SI TOA The pixels larger than the peak value need correction, otherwise, correction is not needed, so that an SIF factor is constructed, and the specific calculation formula of the SIF factor is as follows:
in SI TOA_p Is SI TOA SI at histogram peak TOA Value, SI TOA_max Is SI TOA Is a maximum value of (a).
4. The method for correcting the topographic shading effect of the remote sensing image of an optical satellite according to claim 1, wherein the VIF factor in the step S4 is calculated according to the following formula:
wherein ρ is NIR And ρ R The surface reflectivities of the NIR band and the Red band, respectively.
5. The method according to claim 1, wherein in the step S5, when the average values of the reflectivities of the sunny slope and the cloudy slope are equal, a rough BAF value is obtained.
6. The method for correcting the topographic shading effect of the remote sensing image of an optical satellite according to claim 1, wherein the corrected surface reflectivity in the step S6 is calculated according to the following formula:
wherein ρ is corr Is the corrected pixel reflectivity ρ orig Is the pixel reflectivity before correction, f (SIF, VIF, BAF) is the self-adjusting direct light reflectivity compensation intensity function calculated by multiplying the values of three factors equal to SIF, VIF, BAF, E d Is the direct solar radiation, E, in said step S2 f Is the solar scattered radiation in said step S2, E is the sum of the direct solar radiation, the sky scattered radiation, the topography reflected radiation in said step S2, ρ a Is the average reflectivity of neighboring pixels of each pixel.
7. The method for correcting the topographic shading effect of the optical satellite remote sensing image according to claim 1, wherein the step S7 is specifically:
step S71, respectively selecting vegetation sample points with the same quality as a shade slope and a sunny slope in a research area;
step S72, calculating the correction result of each sample point under each BAF value from 0 by taking 0.05 as step, and calculating the correlation coefficient r of the positive slope sample reflectivity and cosi before correction 1 ；
Step S73, corrected SI TOA_p +2 times SI TOA Standard deviation of values (noted as SI TOA_p +2σ) correlation coefficient r of reflectance and cosi of all samples of the left region 2 ；
Step S74, when r 1 And r 2 When the difference value of (2) approaches 0, the BAF value is the optimal value;
wherein, in the above steps, cosi is calculated according to the following formula:
wherein Z is the zenith angle of the sun, S is the gradient,is the solar azimuth angle and β is the slope direction.
8. The method according to claim 1, wherein in the step S8, the qualitative analysis is visual inspection analysis; quantitative analysis included: and comparing the decorrelation analysis among the reflectivity of evergreen broad-leaved forests on the hillside and the sunny slope, the surface reflectivity and the illumination condition, and the NDVI correction stability analysis, the principal image and the falling image correction performance analysis.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310831948.0A CN116862798B (en) | 2023-07-07 | 2023-07-07 | Topographic shadow effect correction method for optical satellite remote sensing image |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310831948.0A CN116862798B (en) | 2023-07-07 | 2023-07-07 | Topographic shadow effect correction method for optical satellite remote sensing image |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116862798A true CN116862798A (en) | 2023-10-10 |
| CN116862798B CN116862798B (en) | 2024-01-12 |
Family
ID=88222951
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310831948.0A Active CN116862798B (en) | 2023-07-07 | 2023-07-07 | Topographic shadow effect correction method for optical satellite remote sensing image |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116862798B (en) |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP1653193A1 (en) * | 2004-10-29 | 2006-05-03 | Deutsches Zentrum für Luft- und Raumfahrt e.V. | Method for elimination of shadow effects in remote sensing data over land |
| CN114596234A (en) * | 2022-03-21 | 2022-06-07 | 昆明理工大学 | NDVI terrain shadow effect correction method for complex mountainous region |
| CN114778483A (en) * | 2022-04-25 | 2022-07-22 | 福州大学 | Method for correcting terrain shadow of remote sensing image near-infrared wave band for monitoring mountainous region |
-
2023
- 2023-07-07 CN CN202310831948.0A patent/CN116862798B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP1653193A1 (en) * | 2004-10-29 | 2006-05-03 | Deutsches Zentrum für Luft- und Raumfahrt e.V. | Method for elimination of shadow effects in remote sensing data over land |
| CN114596234A (en) * | 2022-03-21 | 2022-06-07 | 昆明理工大学 | NDVI terrain shadow effect correction method for complex mountainous region |
| CN114778483A (en) * | 2022-04-25 | 2022-07-22 | 福州大学 | Method for correcting terrain shadow of remote sensing image near-infrared wave band for monitoring mountainous region |
Non-Patent Citations (1)
| Title |
|---|
| 吴志杰;徐涵秋;: "卫星影像数据构建山地植被指数与应用分析", 地球信息科学学报, vol. 13, no. 05, pages 656 - 663 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116862798B (en) | 2024-01-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109581372B (en) | Ecological environment remote sensing monitoring method | |
| WO2020015326A1 (en) | Remote sensing image cloud shadow detection method supported by earth surface type data | |
| Li et al. | A physics-based atmospheric and BRDF correction for Landsat data over mountainous terrain | |
| JP4425983B1 (en) | Method and apparatus for evaluating solar radiation | |
| Ma et al. | Determination of regional net radiation and soil heat flux over a heterogeneous landscape of the Tibetan Plateau | |
| CN108986040B (en) | NDVI shadow influence removing method based on remote sensing multispectral image | |
| CN111795936A (en) | Multispectral remote sensing image atmospheric correction system and method based on lookup table and storage medium | |
| WO2004069305A2 (en) | Method for performing automated in-scene based atmospheric compensation for multi- and hyperspectral imaging sensors in the solar reflective spectral region | |
| CN109325973B (en) | Urban river network area water body atmosphere correction method | |
| CN110751727A (en) | Synthetic image construction method based on Landsat long-time sequence | |
| CN104966298A (en) | Night low-cloud heavy-mist monitoring method based on low-light cloud image data | |
| CN114778483A (en) | Method for correcting terrain shadow of remote sensing image near-infrared wave band for monitoring mountainous region | |
| CN113970376A (en) | Satellite infrared load calibration method based on ocean area reanalysis data | |
| CN113379759A (en) | Automatic water body extraction method for optical remote sensing satellite image | |
| Santosa | Evaluation of satellite image correction methods caused by differential terrain illumination | |
| CN116664947A (en) | Blue algae bloom monitoring method and system based on satellite observation data | |
| CN111354054A (en) | Polar region visible light remote sensing self-adaptive mapping method | |
| CN116862798B (en) | Topographic shadow effect correction method for optical satellite remote sensing image | |
| CN108198178B (en) | Method and device for determining atmospheric range radiation value | |
| CN114596234B (en) | NDVI terrain shadow effect correction method for complex mountainous region | |
| CN108399363B (en) | Cloud detection method based on polarization image | |
| Schläpfer et al. | Correction of shadowing in imaging spectroscopy data by quantification of the proportion of diffuse illumination | |
| CN113763272A (en) | Remote sensing inversion method for photosynthetic effective radiation attenuation coefficient of eutrophic lake | |
| CN113218874A (en) | Method and system for obtaining surface target object reflectivity based on remote sensing image | |
| CN103322946B (en) | A kind of method obtaining porosity of maize canopy |
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 |