A method and system for estimating light energy resource distribution based on cloud frequency

Abstract

A method and system for estimating light energy resource distribution based on cloud frequency, the method comprising the following steps: obtaining a long-time cloud mask image, and calculating a cloud frequency image based on the long-time cloud mask image; calculating a sunshine percentage image of the area to be estimated according to the cloud frequency image of the area to be estimated; obtaining historical sunshine data of multiple meteorological stations in the area to be estimated, and fitting the relationship between the clear sky index and the sunshine percentage of each meteorological station based on the historical sunshine data using different sunshine radiation models; determining an optimal sunshine radiation model corresponding to the area to be estimated based on the fitting accuracy of different models and the stability of empirical coefficients; calculating a terrain shielding factor image of the area to be estimated; and calculating light energy resource distribution data of the area to be estimated based on the optimal sunshine radiation model, the sunshine percentage image of the area to be estimated, and the terrain shielding factor image.

Description

Technical Field

[0001] The present invention relates to the technical field of light energy resource estimation, and in particular to a method and system for estimating light energy resource distribution based on cloud frequency. Background Art

[0002] Solar energy resources are the most abundant of all renewable and sustainable energy resources in the world, and can reduce the ozone layer depletion and climate deterioration caused by traditional energy generation. Solar energy resources play a very important role in new energy strategies. In the past, the evaluation of solar energy resource distribution mainly relied on radiation measurement equipment at meteorological stations. However, due to the high cost and complex measurement technology of solar radiation measurement instruments (such as direct beam solar radiometers and scattered solar radiometers), only a few ground stations have the conditions to measure solar radiation.

[0003] Since the 1920s, many solar radiation estimation models based on conventional meteorological parameters have been proposed, which can be generally divided into theoretical models, empirical models and machine learning models. The theoretical model has a strong physical basis, but due to the complex atmospheric conditions, the model structure is also complex. It is not easy to obtain the atmospheric parameters that affect solar radiation in actual operation. Machine learning models can better solve complex nonlinear problems. However, the model complexity is high and the model training time is long, which is not conducive to the estimation of solar radiation in a large area. Therefore, the application of theoretical models and machine learning models in solar radiation estimation is limited. The empirical model is one of the most commonly used models, and its purpose is to establish the interaction between meteorological elements and radiation through linear or nonlinear regression equations. Among them, the empirical model for estimating light energy resources based on sunshine percentage is simple in form, and the input data is easy to obtain. It has the advantages of high efficiency and short calculation time in estimating light energy resources.

[0004] In the past, when estimating light energy resources based on sunshine percentage, the sunshine hours meteorological data of ground observation stations were usually used. The sunshine hours station observation data is highly credible, but the meteorological stations are unevenly distributed and have poor spatial continuity. With the development of satellite remote sensing technology in recent years, satellites can provide objective data with large-scale periodicity. At present, light energy resource products generated by satellite remote sensing have taken into account the weakening effect of cloud frequency on solar radiation. However, the method of directly estimating solar radiation based on cloud frequency of satellite remote sensing has low accuracy. In addition, the resolution and accuracy of cloud products used in current solar radiation products are insufficient, resulting in the inability to reflect the true distribution of light energy resources in local areas in many areas with highly variable clouds. Summary of the invention

[0005] In view of the above analysis, the embodiments of the present invention aim to provide a method and system for estimating light energy resource distribution based on cloud frequency, so as to solve the problems of low resolution and precision and inaccurate light energy resource distribution estimation of existing cloud products.

[0006] On the one hand, an embodiment of the present invention provides a method for estimating light energy resource distribution based on cloud frequency, comprising the following steps:

[0007] Acquire a long-time cloud mask image, and calculate a cloud frequency image based on the long-time cloud mask image; calculate a sunshine percentage image of the area to be estimated based on the cloud frequency image of the area to be estimated;

[0008] Obtain historical sunshine data of multiple meteorological stations in the area to be estimated, and based on the historical sunshine data, use different sunshine radiation models to fit the relationship between the clear sky index and the sunshine percentage of each meteorological station; determine the optimal sunshine radiation model corresponding to the area to be estimated based on the fitting accuracy of different models and the stability of empirical coefficients;

[0009] Calculate the terrain occlusion factor image of the area to be estimated;

[0010] Based on the optimal sunshine radiation model, the sunshine percentage image of the area to be estimated and the terrain shielding factor image, the light energy resource distribution data of the area to be estimated is calculated.

[0011] Based on the further improvement of the above method, the terrain occlusion factor image of the area to be estimated is calculated, including:

[0012] Calculate the occlusion state of each pixel in the area to be estimated at each moment;

[0013] Hillshade t =(cos(Z s )*cos(Slo))+(sin(Z s )*sin(Slo)*cos(A s - Asp);

[0014]

[0015] Among them, X ij,t Indicates the occlusion state of the position (i, j) at time t, Hillshade t Represents the illumination at the position (i, j) at time t, Z s represents the solar zenith angle at the position (i, j) at time t, Slo represents the slope at the position (i, j), and A s It refers to the solar azimuth at the position (i, j) at time t, and Asp is the slope direction at the position (i, j);

[0016] For each pixel in the estimated area, the occlusion states at multiple moments are averaged to obtain the mean occlusion state of the pixel. The mean occlusion states of all pixels constitute the terrain occlusion factor image.

[0017] Further, based on the optimal sunshine radiation model, the sunshine percentage image of the area to be estimated and the terrain shielding factor image, light energy resource distribution data of the area to be estimated is calculated, including:

[0018] Based on the optimal sunshine radiation model and the sunshine percentage image of the area to be estimated, the initial sunshine radiation of each pixel in the area to be estimated is calculated to obtain initial light energy distribution data;

[0019] According to the formula Calculate the light energy resource distribution data of the area to be estimated; represents the occlusion factor image, ⊙ represents bitwise multiplication, and H represents the initial light energy distribution data.

[0020] Furthermore, after acquiring the long-time cloud mask image and before calculating the cloud frequency image based on the long-time cloud mask image, the method further includes:

[0021] The SWIR / NDSI threshold method is used to identify the ice and snow areas in each cloud mask image, and the ice and snow areas in each cloud mask image are removed;

[0022] The Band2 / 6 ratio method is used to identify the sand and rock areas in each cloud mask image, and the sand and rock areas in each cloud mask image are removed.

[0023] Furthermore, after calculating the cloud frequency image based on the long-time cloud mask image and before calculating the sunshine percentage image of the area to be estimated based on the cloud frequency image of the area to be estimated, the method further includes:

[0024] A stationary noise removal algorithm is used to remove the orbital artifact areas in the cloud frequency image;

[0025] Based on the albedo data product and the water mask product, the following formula is used to identify abnormal areas:

[0026] Area prob =LAV or WS or WA

[0027] LAV=alb min ≥0.08 and alb cv ≥30 and slope≤0 and water=0,

[0028] WS=alb min ≥0.025 and alb cv≥30 and water=1,

[0029] WA=alb min ≥0.075 and water=1,

[0030] Among them, Area prob represents an abnormal area, LAV represents a land area with high albedo and high variability; WS represents a water area with high variability; WA represents a water area with high albedo, slope represents elevation data; alb min Indicates the minimum albedo, alb cv is the coefficient of variation, alb sd Indicates the standard deviation of albedo data, alb mn represents the average value of albedo data;

[0031] The values ​​in the abnormal area corresponding to the cloud frequency image are removed, and the abnormal area is filled using linear interpolation.

[0032] Furthermore, the historical sunshine data includes solar radiation and sunshine hours;

[0033] Based on the historical sunshine data, different sunshine radiation models are used to fit the relationship between the clear sky index and sunshine percentage of each meteorological station, including:

[0034] Based on the solar radiation, sunshine hours, astronomical radiation and available sunshine hours, the relationship between the clear sky index and sunshine percentage of each season at each meteorological station is fitted by using the least squares method based on different sunshine radiation models.

[0035] Furthermore, based on the fitting accuracy of different models and the stability of empirical coefficients, the optimal solar radiation model corresponding to the estimated area is determined, including:

[0036] The coefficient of determination, mean absolute error and root mean square error were used to test the errors of different models.

[0037] Calculate the intra-seasonal standard deviation and inter-seasonal standard deviation of the model coefficients corresponding to each season of different models;

[0038] The model with small error and standard deviation is selected as the optimal solar radiation model.

[0039] Furthermore, a sunshine percentage image of the area to be estimated is calculated based on the cloud frequency image of the area to be estimated, including calculating the sunshine percentage image of the area to be estimated according to the formula S / S0=1-C, wherein C represents the cloud frequency and S / S0 represents the sunshine percentage.

[0040] On the other hand, an embodiment of the present invention provides a light energy resource distribution estimation system based on cloud frequency, the estimation system comprising:

[0041] A sunshine percentage image acquisition module is used to acquire a long-time cloud mask image, calculate a cloud frequency image based on the long-time cloud mask image, and calculate a sunshine percentage image of the area to be estimated based on the cloud frequency image of the area to be estimated;

[0042] The optimal sunshine radiation model determination module is used to obtain the historical sunshine data of multiple meteorological stations in the area to be estimated, and based on the historical sunshine data, use different sunshine radiation models to fit the relationship between the clear sky index and the sunshine percentage of each meteorological station; determine the optimal sunshine radiation model corresponding to the area to be estimated based on the fitting accuracy of different models and the stability of empirical coefficients;

[0043] A terrain occlusion factor calculation module calculates a terrain occlusion factor image of the area to be estimated;

[0044] The light energy resource distribution estimation module is used to calculate the light energy resource distribution data of the area to be estimated based on the optimal sunshine radiation model, the sunshine percentage image of the area to be estimated and the terrain shielding factor image.

[0045] Furthermore, the terrain occlusion factor image of the area to be estimated is calculated, including:

[0046] Calculate the occlusion state of each pixel in the area to be estimated at each moment;

[0047] Hillshade t =(cos(Z s )*cos(Slo))+(sin(Z s )*sin(Slo)*cos(A s - Asp);

[0048]

[0049] Among them, X ij,t Indicates the occlusion state of the position (i, j) at time t, Hillshade t Represents the illumination at the position (i, j) at time t, Z s represents the solar zenith angle at the position (i, j) at time t, Slo represents the slope at the position (i, j), and A s It refers to the solar azimuth at the position (i, j) at time t, and Asp is the slope direction at the position (i, j);

[0050] For each pixel in the estimated area, the occlusion states at multiple moments are averaged to obtain the mean occlusion state of the pixel. The mean occlusion states of all pixels constitute the terrain occlusion factor image.

[0051] The present invention calculates cloud frequency by using a high temporal and spatial resolution and high precision cloud mask generated by satellite remote sensing, and calculates the sunshine percentage of the area to be estimated based on the cloud frequency, thereby eliminating the problem of poor spatial continuity and inaccurate data caused by uneven distribution of meteorological stations. By calculating the sunshine percentage according to the cloud frequency, the light energy distribution data of the area to be estimated is calculated based on the sunshine percentage and the optimal sunshine radiation model, and the calculation result is accurate. By considering terrain shielding, the terrain shielding factor is added to the light energy resource calculation model, thereby further improving the accuracy of light energy resource distribution calculation.

[0052] In the present invention, the above-mentioned technical solutions can also be combined with each other to achieve more preferred combination solutions. Other features and advantages of the present invention will be described in the subsequent description, and some advantages can become obvious from the description, or can be understood by practicing the present invention. The purpose and other advantages of the present invention can be realized and obtained through the contents particularly pointed out in the description and the drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] The drawings are only for the purpose of illustrating particular embodiments and are not to be considered limiting of the present invention. Like reference symbols denote like components throughout the drawings.

[0054] Figure 1 This is a flow chart of a method for estimating light energy resource distribution based on cloud frequency according to an embodiment of the present invention;

[0055] Figure 2 A block diagram of a light energy resource distribution estimation system based on cloud frequency according to an embodiment of the present invention;

[0056] Figure 3 A schematic diagram of global land cloud frequencies over ten days provided by an embodiment of the present invention;

[0057] Figure 4 A schematic diagram of the annual average cloud frequency provided by an embodiment of the present invention;

[0058] Figure 5 A schematic diagram of average cloud frequencies of different cloud products in the Amazon region in 2010 provided by an embodiment of the present invention;

[0059] Figure 6 A monthly distribution map of solar energy resources in China provided by an embodiment of the present invention;

[0060] Figure 7 The annual light energy resource distribution map of China provided by the embodiment of the present invention;

[0061] Figure 8 An optical resource distribution map of Jinchuan County, Sichuan Province provided in an embodiment of the present invention;

[0062] Fig. 9 Schematic diagram of error checking results of different seasons and years for different models according to an embodiment of the present invention;

[0063] Fig.10 It is a schematic diagram of the intra-seasonal standard deviation and inter-seasonal standard deviation results of the Sen model and the Bahel model according to an embodiment of the present invention. DETAILED DESCRIPTION

[0064] The preferred embodiments of the present invention are described in detail below in conjunction with the accompanying drawings, wherein the accompanying drawings constitute a part of this application and are used together with the embodiments of the present invention to illustrate the principles of the present invention, but are not used to limit the scope of the present invention.

[0065] An embodiment of the present invention discloses a method for estimating light energy resource distribution based on cloud frequency, comprising the following steps:

[0066] S1. Acquire a long-time cloud mask image, and calculate a cloud frequency image based on the long-time cloud mask image; calculate a sunshine percentage image of the area to be estimated according to the cloud frequency image of the area to be estimated;

[0067] S2. Obtain historical sunshine data of multiple meteorological stations in the area to be estimated, and based on the historical sunshine data, use different sunshine radiation models to fit the relationship between the clear sky index and the sunshine percentage of each meteorological station; determine the optimal sunshine radiation model corresponding to the area to be estimated based on the fitting accuracy of different models and the stability of empirical coefficients;

[0068] S3, calculating the terrain occlusion factor image of the area to be estimated;

[0069] S4. Calculate the light energy resource distribution data of the area to be estimated based on the optimal sunshine radiation model, the sunshine percentage image of the area to be estimated and the terrain shielding factor image.

[0070] During implementation, the long-term time series cloud mask image is a long-term time series cloud mask image, for example, 16 years of time series cloud mask images from 2001 to 2016 are extracted.

[0071] By using high temporal and spatial resolution and high precision cloud masks generated by satellite remote sensing to calculate cloud frequency, the sunshine percentage of the area to be estimated is calculated based on the cloud frequency, thereby eliminating the problem of poor spatial continuity and inaccurate data caused by uneven distribution of meteorological stations. By calculating the sunshine percentage according to the cloud frequency, the light energy distribution data of the area to be estimated is calculated based on the sunshine percentage and the optimal sunshine radiation model, and the calculation results are accurate. By considering terrain shading, the terrain shading factor is added to the light energy resource calculation model, thereby further improving the accuracy of light energy resource distribution calculation.

[0072] Specifically, long-term cloud mask images are extracted from high-resolution cloud products. When implemented, daily cloud markers are extracted from daily data of MOD09GA and MYD09GA and converted into separate cloud masks. MOD09 can generate cloud masks with a high spatiotemporal resolution of 1 km.

[0073] The cloud detection algorithm of the MOD09 cloud mask itself relies on two reflectance tests (MIRA and 1.38 mic reflectivity) and a thermal test. MIRA can detect low or high reflectance clouds, while 1.38 mic can detect high clouds. In addition, a thermal test is used to identify pixels with abnormal high infrared reflectance (e.g., fires, sun glints, and high albedo surfaces). After these two reflectance and thermal tests, the MOD09 cloud mask still has many uncertainties, such as bright rocks and sand being misclassified as clouds. The identification of ice and snow in the MOD09 cloud algorithm is achieved by considering the ground temperature. However, at high latitudes or high altitudes, the surface temperature cannot effectively distinguish ice and snow from clouds. In the MOD09 cloud mask, ice and snow and high reflectivity areas are the main components of cloud misdetection. In the past, the MOD09 cloud mask was not considered in the application of the cloud mask itself, resulting in obvious abnormal cloud frequency in the generated long-term cloud climate data. Therefore, before performing the cloud frequency calculation, it is necessary to add a spectral test to minimize the misclassification of ice and snow and bright areas as clouds.

[0074] Therefore, in step S1, after acquiring the long-time cloud mask image and before calculating the cloud frequency image based on the long-time cloud mask image, the method further includes:

[0075] S11. Using the SWIR / NDSI threshold method to identify the ice and snow area in each cloud mask image, and removing the ice and snow area in each cloud mask image;

[0076] Ice and snow have high reflectivity in the visible light band and strong absorption in the shortwave infrared band, while clouds have high reflectivity in both the visible light band and the shortwave infrared band. Therefore, based on this characteristic, the SWIR / NDSI threshold method is used to identify ice and snow in the cloud mask. The specific identification formula is:

[0077] Snow / ice Test=Band 7>0.025 and NDSI<0.8

[0078] NDSI=(Band 4-Band 6) / (Band 4+Band 6)

[0079] Band 7 is the shortwave infrared band, and NDSI is the normalized snow index. In principle, Band 7 should be less than 0.03. In order to prevent thin clouds from being removed, the embodiment of the present invention sets Band 7 to 0.025. Generally speaking, the NDSI value of clouds is around 0, but the value of ice clouds will be larger, but will not exceed 0.8. Therefore, during implementation, areas with Band 7>0.025 and NDSI<0.8 are identified as ice and snow areas. Band 1-7 refers to the reflectivity of bands 1-7.

[0080] MOD09GA comes with reflectivity data. The ice and snow areas are identified according to the above recognition formula, and the cloud marks of the corresponding areas are removed in the cloud mask image, that is, the ice and snow areas in each cloud mask are removed.

[0081] S12. Using the Band2 / 6 ratio method, the sand and stone areas in each cloud mask image are identified, and the sand and stone areas in each cloud mask image are removed.

[0082] In addition to the confusion between ice and clouds, bright rocks and sand often appear in cloud masks, especially in desert areas. The reflectivity of sand and rocks in Band 6 is often higher than that in Band 2, resulting in the ratio of Band 2 to Band 6 of sand and stone usually being less than 1, while the reflectivity of clouds is the opposite. Based on this feature, the ratio of Band 2 to Band 6 can be used to eliminate sand and rocks with high reflectivity. The formula is as follows:

[0083] B2 / B6 Test=Band 2 / Band 6>0.85

[0084] The threshold is adjusted to 0.85 here because sand and gravel will be exposed under thin clouds, which will reduce the ratio. According to the above recognition formula, the sand and gravel areas with Band 2 / Band 6>0.85 are identified, and the cloud marks of the corresponding areas are removed in the cloud mask image, that is, the sand and gravel areas in each cloud mask are removed.

[0085] By removing the ice, snow and sand areas from the cloud mask image, the accuracy of the cloud mask is improved, laying the foundation for the subsequent generation of high temporal and spatial resolution and high-precision long-term cloud frequency data.

[0086] Based on the improved long-time cloud mask image, the cloud frequency image is calculated.

[0087] During implementation, since the MODIS sensors are carried on Terra and Aqua satellites, the cloud mask data obtained include MOD09 data and MYD09 data. The MOD09 cloud mask and MYD09 cloud mask are calculated separately during calculation. The average cloud frequency of each decade from 2001 to 2016 is calculated, and then the 16-year data of the same decade are averaged to obtain 72 16-year average ten-day cloud frequencies. The MOD09 cloud frequency and MYD09 cloud frequency of the same period are averaged to obtain 36 long-term ten-day cloud frequency images, thereby obtaining a long-term ten-day cloud frequency product with high temporal and spatial resolution and high accuracy. The ten-day cloud frequency is the number of days with clouds in 10 days / 10.

[0088] The generated cloud frequency image still has some abnormal areas, such as MODIS orbit artifacts and abnormal cloud frequencies in high albedo areas. Therefore, before calculating the sunshine percentage image of the area to be estimated based on the cloud frequency image of the area to be estimated, it also includes:

[0089] S13. Use a stationary noise removal algorithm to remove track artifact areas in the cloud frequency image.

[0090] The MODIS orbit will cause a systematic bias in the daily global coverage product near the equator, which will cause longitudinal orbit artifacts in the long-term cloud frequency. In order to remove these stripe noises, the present invention uses a stationary noise removal algorithm (VSNR) to process the produced long-term cloud frequency data. VSNR can retain image details and remove stripes and most of the white noise.

[0091] Since the stripe noise shapes and proportions of the MOD09 cloud frequency image and the MYD09 cloud frequency image are different, step S13 is performed before the MOD09 cloud frequency and the MYD09 cloud frequency of the same period are averaged to obtain a 36-period long-term ten-day cloud frequency image.

[0092] In the implementation, images with obvious streak noise were screened out from 72 cloud frequency images, and the problem areas in these images were cropped out as the input data of VSNR. In the specific processing of long-term cloud frequency data, VSNR removed the known streak noise by setting the shape and ratio of the noise. For the MOD09 cloud frequency image, the size ratio of the noise was set to y = 200 and x = 5, and the angle was set to θ = 15. For the MYD09 cloud frequency image, the angle of the noise was set to θ = -15 to minimize the orbital artifacts in the eastern Australian and central Pacific regions. In the northwest of Africa and the northern Indian Ocean region, the noise was set to y = 40, x = 5 and θ = 15 (MOD09 cloud frequency image) / θ = -15 (MYD09 cloud frequency image). The output image after removing the streak noise was re-mosaiced into the long-term cloud frequency image.

[0093] S14. Based on the albedo data product and the water mask product, the following formula is used to identify abnormal areas:

[0094] Area prob =LAV or WS or WA

[0095] LAV=alb min ≥0.08 and alb cv ≥30 and slope≤0 and water=0,

[0096] WS=alb min ≥0.025 and alb cv ≥30 and water=1,

[0097] WA=alb min ≥0.075 and water=1,

[0098] Among them, Area prob represents an abnormal area, LAV represents a land area with high albedo and high variability; WS represents a water area with high variability; WA represents a water area with high albedo, slope represents elevation data; alb min Indicates the minimum albedo, alb cv is the coefficient of variation, alb sd Indicates the standard deviation of albedo data, alb mn Represents the average value of albedo data.

[0099] During implementation, the 1km MODIS MCD43B3 product, MODIS MOD44W water mask product and SRTM DEM can be used to identify abnormal areas. Since the cloud mask data is obtained based on 16 years of data from 2000 to 2016, the albedo data of the same period is extracted. Therefore, alb in the above formula is min The minimum albedo of the pixel within 16 years, alb sd is the standard deviation of the albedo data of the pixel point for 16 years, alb mn It represents the average albedo data of the pixel point in 16 years. Slope represents the elevation data of the pixel point, which can be calculated based on SRTM DEM. For water, 0 and 1 represent non-water surface and water surface respectively.

[0100] These problem areas are generally divided into three categories: (1) land with high albedo and high variability (LAV), (2) water with high variability (WS), and (3) water with high albedo (WA).

[0101] S13, removing the values ​​in the abnormal area corresponding to the cloud frequency image, and filling the abnormal area with a linear interpolation method, so as to obtain an accurate cloud frequency image.

[0102] Specifically, in step S1, the sunshine percentage image of the area to be estimated is calculated based on the cloud frequency image of the area to be estimated, including calculating the sunshine percentage image of the area to be estimated according to the formula S / S0=1-C, where C represents the cloud frequency and S / S0 represents the sunshine percentage. The sunshine percentage is the percentage of the actual sunshine hours S to the available sunshine hours S0 in a period of time.

[0103] In the absence of visual errors and cloud estimation errors, the daytime cloud cover is inversely proportional to the sunshine percentage. Therefore, in the absence of observation errors, the relationship between the sunshine percentage (S / S0) and the cloud frequency (C) is linear, that is, C = 1-S / S0. Therefore, the sunshine percentage of each pixel in the area to be estimated can be calculated according to the formula S / S0 = 1-C, so as to generate the sunshine percentage image of the area to be estimated according to the cloud frequency image, that is, to obtain 36 periods of sunshine percentage images.

[0104] The present invention uses the ten-day cloud frequency to generate the sunshine percentage. The daily observation data of cloud frequency, whether it is ground observation or satellite remote sensing data, has a problem of observation time, that is, the time term is discrete within a single day. The present invention uses the observation twice a day for calculation. Although the step size is not as small as the ground observation data, it uses the observation time range of every 10 days in 16 years. The time range is large enough to make up for the scarcity on the single-day time scale with the surplus on the long time series scale.

[0105] The empirical model for estimating light energy resources based on sunshine percentage is to construct a linear or nonlinear relationship between the ratio of total radiation H to astronomical radiation H0 (clear sky index) and sunshine percentage.

[0106] In step S2, the historical sunshine data includes solar radiation and sunshine hours.

[0107] Based on the historical sunshine data, different sunshine radiation models are used to fit the relationship between the clear sky index and sunshine percentage of each meteorological station, including:

[0108] S21. Based on the solar radiation, sunshine hours, astronomical radiation and available sunshine hours, the relationship between the clear sky index and sunshine percentage of each season at each meteorological station is fitted using the least squares method based on different sunshine radiation models.

[0109] Specifically, the astronomical radiation and sunshine hours of the weather station can be calculated based on basic astronomical and geographical data. The specific calculation process is referred to the prior art and will not be repeated here. The solar radiation and sunshine hours are measured by each weather station.

[0110] During implementation, the solar radiation models include AP model, Iqbal model, Bahel model and Sen model. The least squares fitting method is used to fit the relationship between the clear sky index and the sunshine percentage in each season of each meteorological station based on different solar radiation models.

[0111] For example, the formula of Sen model is expressed as Taking season as time period, the least square method is used to fit the solar radiation, sunshine hours, astronomical radiation and available sunshine hours of each meteorological station in the season to solve the optimal coefficients a, b, c of the model. The optimal coefficients of different meteorological stations in different seasons may be different. Among them, H represents solar radiation, H0 represents astronomical radiation, S represents actual sunshine hours, and S0 represents available sunshine hours.

[0112] In step S2, the optimal solar radiation model corresponding to the area to be estimated is determined based on the fitting accuracy of different models and the stability of the empirical coefficients, including:

[0113] S22, use the coefficient of determination, mean absolute error and root mean square error to test the error of different models;

[0114] S23, calculating the intra-seasonal standard deviation and inter-seasonal standard deviation of the model coefficients corresponding to each season of different models;

[0115] S24. Select a model with small error and standard deviation as the optimal solar radiation model.

[0116] For example, the error checking results of each model in different seasons and years are as follows: Fig. 9 As shown. Taking MAE as an example, the light gray one is the Sen model, which has the smallest spring, summer, autumn, winter and annual MAE values. The smaller the MAE value, the smaller the error. The same is true for RMSE. The larger the coefficient of determination R2, the smaller the error, and the Sen model is the largest.

[0117] For example, Fig.10 The intra-seasonal standard deviation and inter-seasonal standard deviation (i.e., annual standard deviation) of the Sen model and the Bahel model are shown in Figure 2. The standard deviation of the coefficients of the Sen model in all seasons is much smaller than that of the Bahel model, so the coefficients are more stable. The Sen model can be selected as the optimal model based on the error and standard deviation.

[0118] When performing least squares fitting on meteorological data, empirical coefficients are solved according to seasons rather than ten-day periods. This is because the empirical coefficients solved based on 16 years and 10 days of data cannot guarantee their stability and the accuracy of the final results. The empirical coefficients of the model divided by season are suitable for the prediction of light resource distribution in each ten-day period of the season. The empirical coefficients are more stable and the estimated light resource distribution is more accurate.

[0119] After determining the optimal model for the area to be estimated, for other pixels in the area to be estimated except the location of the meteorological station, the model coefficient of the pixel is obtained by performing inverse distance weighted interpolation calculation based on the coefficients of all meteorological stations.

[0120] After the model coefficient of each pixel is determined, the solar radiation of the pixel can be calculated according to the sunshine percentage of the pixel calculated in step S1 and the model coefficient of the pixel.

[0121] For example, the optimal model is the Sen model. H0, a, b, c and All are known, so H can be calculated. It should be noted that since the model coefficients of different seasons are different, the 36-period sunshine percentage images obtained above are brought into the model of the corresponding season to calculate the daily radiation.

[0122] When sunlight passes through clouds and reaches the ground, it is affected by the local terrain of the ground, and the solar radiation reaching the ground is redistributed on the local surface. In mountainous areas with complex terrain, different slopes often receive different amounts of radiation, which is caused by slope gradient, slope aspect and ground shielding. Traditional solar radiation estimation models do not take into account the shielding of solar radiation by terrain.

[0123] Therefore, the present invention needs to consider the shielding factor when calculating the final solar radiation, so as to make the solar radiation result more accurate. Specifically, calculating the terrain shielding factor image of the area to be estimated includes:

[0124] Calculate the occlusion state of each pixel in the area to be estimated at each moment;

[0125] Hillshade t =(cos(Z s )*cos(Slo))+(sin(Z s )*sin(Slo)*cos(A s - Asp);

[0126]

[0127] Among them, X ij,tIndicates the occlusion state of the position (i, j) at time t, Hillshade t Represents the illumination at the position (i, j) at time t, Z s represents the solar zenith angle at the position (i, j) at time t, Slo represents the slope at the position (i, j), and A s It refers to the solar azimuth at the position (i, j) at time t, and Asp is the slope direction at the position (i, j);

[0128] For each pixel in the estimated area, the occlusion states at multiple moments are averaged to obtain the mean occlusion state of the pixel. The mean occlusion states of all pixels constitute the terrain occlusion factor image.

[0129] During implementation, for example, every 10 minutes can be regarded as a moment, and the occlusion state map of all moments within the time range of 36 ten days can be calculated. For each pixel point, the sum of the occlusion state values ​​of all moments is calculated, and divided by the number of moments included in each ten days (for example, 6×24×10=1440 moments), and the average occlusion state of each ten days is obtained. The average occlusion state of all pixels forms the occlusion factor image.

[0130] Specifically, in step S4, based on the optimal sunshine radiation model, the sunshine percentage image of the area to be estimated and the terrain shielding factor image, the light energy resource distribution data of the area to be estimated is calculated, including:

[0131] Based on the optimal sunshine radiation model and the sunshine percentage image of the area to be estimated, the initial sunshine radiation of each pixel in the area to be estimated is calculated to obtain initial light energy distribution data;

[0132] According to the formula Calculate the light energy resource distribution data of the area to be estimated; represents the occlusion factor image, ⊙ represents bitwise multiplication, and H represents the initial light energy distribution data.

[0133] The effects of the embodiments of the present invention are described in detail below.

[0134] Figure 3 is an example of a global ten-day land cloud frequency plot showing non-global land areas excluding the polar regions. From left to right and top to bottom, the 2nd, 5th, 8th, 11th, 14th, 17th, 20th, 23rd, 26th, 29th, 32nd and 35th ten-day periods of the year are shown, respectively, for a 16-year average cloud frequency. Figure 4 The annual average cloud frequency of the present invention is shown. The statistical results of cloud distribution on land except the polar regions show that the annual average cloud cover is about 50.4%, which meets the definition of the global land annual average cloud frequency of 50% to 60%.

[0135] Other cloud coverage products are introduced as comparative data, and the Amazon River region is taken as an example to demonstrate the accuracy of the cloud frequency data of the present invention. Figure 5 The Amazon River Basin Figure 5 a), MOD / MYD35( Figure 5 b) and PATMOS-X( Figure 5 c) Cloud distribution. The present invention and MOD / MYD35 show that the cloud frequency in the river area is lower than that in the surrounding area. However, PATMOS-X shows a completely opposite result to GLHCC and MOD / MYD35, with a higher cloud frequency over the river. By studying the true color images of the Amazon basin, it was determined that the cloud frequency of PATMOS-X is inconsistent with the facts. Cumulus clouds are everywhere over the rainforests in the Amazon basin, but there are no clouds over the rivers. This phenomenon reveals the transpiration of tropical rainforests, which is also common in tropical areas such as the Congo River Basin. In addition, the cloud frequency of the present invention is also different from that of MOD35. See the enlarged image in the lower right corner of the two for details, which shows the area in the red box. The river channel in the MOD35 cloud frequency shows a high cloud frequency, while the cloud frequency of the present invention is normal. This is because the MOD35 cloud product identifies the exposed river channel as a cloud.

[0136] Taking China as an example, the light energy resource distribution based on cloud frequency is estimated, and the estimation results are displayed in monthly and annual forms, such as Figure 6 and Figure 7 According to statistics, the annual average daily total solar radiation in China is 5213MJ / m 2 , in line with the annual total solar radiation of 5000MJ / m 2 The above definition.

[0137] In order to prove that the light energy resource estimation of the present invention will bring more accurate estimation results, the solar radiation data of the World Bank and the irradiance data of the photovoltaic geographic information system were selected for visual comparison. Figure 8 a is the DEM of Jinchuan County with a resolution of 90m. Figure 8 b is the light resource distribution of the present invention (after normalization), Figure 8 c is the World Bank's GHI (normalized), Figure 8d is the irradiance of PVGIS (normalized). It can be found from the figure that due to the difference in resolution, the World Bank's GHI is more detailed than the light energy resources and irradiance of PVGIS of the present invention. The distribution of light energy resources in the present invention is closer to the irradiance of PVGIS, with relatively abundant light resources in the river valley area and relatively scarce light resource distribution on the hillside. However, the GHI of the region released by the World Bank is contrary to the light resources and irradiance of PVGIS of the present invention, and there is an obvious positive correlation between light energy resources and altitude. After investigation, it was found that the distribution of light resources of the present invention is more in line with the actual situation than the irradiance of PVGIS. Because the river valley in this area is a typical dry and hot river valley, it is characterized by dryness, heat, and long sunshine time. The high resolution of GHI comes from SRTM DEM data, rather than satellite data or other data sources used to calculate the cloud index. Although most of the current solar radiation products take into account the influence of cloud cover, the cloud cover data used often have low spatial resolution or low accuracy. Therefore, although the GHI data have a higher resolution and can reflect more terrain details, the distribution of light energy resources is not accurate in areas with clouds with high spatial variability.

[0138] The present invention uses the SWIR / NDSI threshold method and the Band 2 / 6 index method to improve the MOD09 cloud mask, reducing the confusion between ice and snow, bright areas and clouds in the cloud mask, so that the long-term cloud climate data can more accurately reflect the cloud frequency distribution. Based on the high temporal and spatial resolution and high-precision cloud frequency data, the sunshine percentage is integrated, the optimal model is used, and the distribution of light energy resources is estimated in combination with the terrain shielding factor, which greatly improves the accuracy of light energy resource data, especially in areas with high spatial variability clouds.

[0139] A specific embodiment of the present invention discloses a light energy resource distribution estimation system based on cloud frequency, such as Figure 2 As shown, it includes the following modules:

[0140] A sunshine percentage image acquisition module is used to acquire a long-time cloud mask image, calculate a cloud frequency image based on the long-time cloud mask image, and calculate a sunshine percentage image of the area to be estimated based on the cloud frequency image of the area to be estimated;

[0141] The optimal sunshine radiation model determination module is used to obtain the historical sunshine data of multiple meteorological stations in the area to be estimated, and based on the historical sunshine data, use different sunshine radiation models to fit the relationship between the clear sky index and the sunshine percentage of each meteorological station; determine the optimal sunshine radiation model corresponding to the area to be estimated based on the fitting accuracy of different models and the stability of empirical coefficients;

[0142] A terrain occlusion factor calculation module calculates a terrain occlusion factor image of the area to be estimated;

[0143] The light energy resource distribution estimation module is used to calculate the light energy resource distribution data of the area to be estimated based on the optimal sunshine radiation model, the sunshine percentage image of the area to be estimated and the terrain shielding factor image.

[0144] Preferably, calculating the terrain occlusion factor image of the area to be estimated includes:

[0145] Calculate the occlusion state of each pixel in the area to be estimated at each moment;

[0146] Hillshade t =(cos(Z s )*cos(Slo))+(sin(Z s )*sin(Slo)*cos(A s - Asp);

[0147]

[0148] Among them, X ij,t Indicates the occlusion state of the position (i, j) at time t, Hillshade t Represents the illumination at the position (i, j) at time t, Z s represents the solar zenith angle at the position (i, j) at time t, Slo represents the slope at the position (i, j), and A s It refers to the solar azimuth at the position (i, j) at time t, and Asp is the slope direction at the position (i, j);

[0149] For each pixel in the estimated area, the occlusion states at multiple moments are averaged to obtain the mean occlusion state of the pixel. The mean occlusion states of all pixels constitute the terrain occlusion factor image.

[0150] The above method embodiments and system embodiments are based on the same principle, and their related parts can be used for reference, and can achieve the same technical effect. The specific implementation process refers to the above embodiment, which will not be repeated here.

[0151] Those skilled in the art will appreciate that all or part of the processes of the above-mentioned embodiments can be implemented by instructing related hardware through a computer program, and the program can be stored in a computer-readable storage medium, wherein the computer-readable storage medium is a disk, an optical disk, a read-only storage memory, or a random access memory, etc.

[0152] The above description is only a preferred specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by any technician familiar with the technical field within the technical scope disclosed by the present invention should be covered within the protection scope of the present invention.

Claims

1. A method for estimating light energy resource distribution based on cloud frequency, characterized in that: The following steps are involved: Acquire a long-time cloud mask image, and calculate a cloud frequency image based on the long-time cloud mask image; calculate a sunshine percentage image of the area to be estimated based on the cloud frequency image of the area to be estimated; Obtain historical sunshine data of multiple meteorological stations in the area to be estimated, and based on the historical sunshine data, use different sunshine radiation models to fit the relationship between the clear sky index and sunshine percentage of each meteorological station; Determine the optimal solar radiation model corresponding to the area to be estimated based on the fitting accuracy of different models and the stability of empirical coefficients; Calculate the terrain occlusion factor image of the area to be estimated; Calculating light energy resource distribution data of the area to be estimated based on the optimal sunshine radiation model, the sunshine percentage image of the area to be estimated and the terrain shielding factor image; Calculate the terrain occlusion factor image of the area to be estimated, including: Calculate the occlusion state of each pixel in the area to be estimated at each moment; Hillshade t =(cos(Z s )*cos(Slo))+(sin(Z s )*sin(Slo)*cos(A s -Asp)); Among them, X ij,t Indicates the occlusion state of the position (i, j) at time t, Hillshade t Represents the illumination at the position (i, j) at time t, Z s represents the solar zenith angle at the position (i, j) at time t, Slo represents the slope at the position (i, j), and A s It refers to the solar azimuth at the position (i, j) at time t, and Asp is the slope direction at the position (i, j); For each pixel in the estimated area, the occlusion states at multiple moments are averaged to obtain the average occlusion state of the pixel. The average occlusion states of all pixels constitute a terrain occlusion factor image. Calculating light energy resource distribution data of the area to be estimated based on the optimal sunshine radiation model, the sunshine percentage image of the area to be estimated, and the terrain shielding factor image, including: Based on the optimal sunshine radiation model and the sunshine percentage image of the area to be estimated, the initial sunshine radiation of each pixel in the area to be estimated is calculated to obtain initial light energy distribution data; According to the formula Calculate the light energy resource distribution data of the area to be estimated; represents the occlusion factor image, ⊙ represents bitwise multiplication, and H represents the initial light energy distribution data; The historical sunshine data includes solar radiation and sunshine hours; Based on the historical sunshine data, different sunshine radiation models are used to fit the relationship between the clear sky index and sunshine percentage of each meteorological station, including: Based on the solar radiation, sunshine hours, astronomical radiation and available sunshine hours, the relationship between the clear sky index and sunshine percentage of each season at each meteorological station is fitted using the least square method based on different sunshine radiation models; Based on the fitting accuracy of different models and the stability of empirical coefficients, the optimal solar radiation model corresponding to the estimated area is determined, including: The coefficient of determination, mean absolute error and root mean square error were used to test the errors of different models. Calculate the intra-seasonal standard deviation and inter-seasonal standard deviation of the model coefficients corresponding to each season of different models; Select the model with small error and standard deviation as the optimal solar radiation model; Calculating the sunshine percentage image of the area to be estimated based on the cloud frequency image of the area to be estimated includes calculating the sunshine percentage image of the area to be estimated according to the formula S / S0=1-C, where C represents the cloud frequency and S / S0 represents the sunshine percentage.

2. The method for estimating light energy resource distribution based on cloud frequency according to claim 1, characterized in that: After acquiring the long-time cloud mask image and before calculating the cloud frequency image based on the long-time cloud mask image, the method further includes: The SWIR / NDSI threshold method is used to identify the ice and snow areas in each cloud mask image, and the ice and snow areas in each cloud mask image are removed; The Band2 / 6 ratio method is used to identify the sand and rock areas in each cloud mask image, and the sand and rock areas in each cloud mask image are removed.

3. The method for estimating light energy resource distribution based on cloud frequency according to claim 1, characterized in that: After calculating the cloud frequency image based on the long time series cloud mask image and before calculating the sunshine percentage image of the area to be estimated based on the cloud frequency image of the area to be estimated, the method further includes: A stationary noise removal algorithm is used to remove the orbital artifact areas in the cloud frequency image; Based on the albedo data product and the water mask product, the following formula is used to identify abnormal areas: Area prob =LAV or WS or WA LAV=alb min ≥0.08 and alb cv ≥30 and slope≤0 and water=0, WS=alb min ≥0.025 and alb cv ≥30 and water=1, WA=alb min ≥0.075 and water=1, Among them, Area prob represents an abnormal area, LAV represents a land area with high albedo and high variability; WS represents a water area with high variability; WA represents a water area with high albedo, slope represents elevation data; alb min Indicates the minimum albedo, alb cv is the coefficient of variation, alb sd Indicates the standard deviation of albedo data, alb mn represents the average value of albedo data; The values ​​in the abnormal area corresponding to the cloud frequency image are removed, and the abnormal area is filled using linear interpolation.

4. A light energy resource distribution estimation system based on cloud frequency, characterized in that: The estimation system comprises: A sunshine percentage image acquisition module is used to acquire a long-time cloud mask image, calculate a cloud frequency image based on the long-time cloud mask image, and calculate a sunshine percentage image of the area to be estimated based on the cloud frequency image of the area to be estimated; The optimal sunshine radiation model determination module is used to obtain the historical sunshine data of multiple meteorological stations in the area to be estimated, and adopt different sunshine radiation models based on the historical sunshine data to simulate the The relationship between the clear sky index and sunshine percentage of each meteorological station is combined; the optimal sunshine radiation model corresponding to the estimated area is determined based on the fitting accuracy of different models and the stability of empirical coefficients; A terrain occlusion factor calculation module calculates a terrain occlusion factor image of the area to be estimated; A light energy resource distribution estimation module, used to calculate light energy resource distribution data of the area to be estimated based on the optimal sunshine radiation model, the sunshine percentage image of the area to be estimated and the terrain shielding factor image; Calculate the terrain occlusion factor image of the area to be estimated, including: Calculate the occlusion state of each pixel in the area to be estimated at each moment; Hillshade t =(cos(Z s )*cos(Slo))+(sin(Z s )*sin(Slo)*cos(A s -Asp)); Among them, X ij,t Indicates the occlusion state of the position (i, j) at time t, Hillshade t Represents the illumination at the position (i, j) at time t, Z s represents the solar zenith angle at the position (i, j) at time t, Slo represents the slope at the position (i, j), and A s It refers to the solar azimuth at the position (i, j) at time t, and Asp is the slope direction at the position (i, j); For each pixel in the estimated area, the occlusion states at multiple moments are averaged to obtain the average occlusion state of the pixel. The average occlusion states of all pixels constitute a terrain occlusion factor image. Calculating light energy resource distribution data of the area to be estimated based on the optimal sunshine radiation model, the sunshine percentage image of the area to be estimated, and the terrain shielding factor image, including: Based on the optimal sunshine radiation model and the sunshine percentage image of the area to be estimated, the initial sunshine radiation of each pixel in the area to be estimated is calculated to obtain initial light energy distribution data; According to the formula Calculate the light energy resource distribution data of the area to be estimated; represents the occlusion factor image, ⊙ represents bitwise multiplication, and H represents the initial light energy distribution data; The historical sunshine data includes solar radiation and sunshine hours; Based on the historical sunshine data, different sunshine radiation models are used to fit the relationship between the clear sky index and sunshine percentage of each meteorological station, including: Based on the solar radiation, sunshine hours, astronomical radiation and available sunshine hours, the relationship between the clear sky index and sunshine percentage of each season at each meteorological station is fitted using the least square method based on different sunshine radiation models; Based on the fitting accuracy of different models and the stability of empirical coefficients, the optimal solar radiation model corresponding to the estimated area is determined, including: The coefficient of determination, mean absolute error and root mean square error were used to test the errors of different models. Calculate the intra-seasonal standard deviation and inter-seasonal standard deviation of the model coefficients corresponding to each season of different models; Select the model with small error and standard deviation as the optimal solar radiation model; Calculating the sunshine percentage image of the area to be estimated based on the cloud frequency image of the area to be estimated includes calculating the sunshine percentage image of the area to be estimated according to the formula S / S0=1-C, where C represents the cloud frequency and S / S0 represents the sunshine percentage.

Images