A method for irradiance optimization and correction based on terrain and mountain shadow discrimination and complex terrain.

Abstract

This invention relates to the field of solar radiation resource assessment technology. Its purpose is to provide a method for optimizing and correcting irradiance based on terrain and mountain shadow discrimination and complex terrain, addressing the problem of large errors in existing calculations of total horizontal radiation data for complex terrain due to simplistic calculations. The technical solution can be summarized as follows: In complex terrain areas, firstly, Hillshade grids and illumination conditions are identified. Then, based on data availability, direct normal radiation and diffuse radiation are obtained. For direct normal radiation, terrain visibility, albedo parameters, and the Minnaert non-Lambertian correction coefficient k are introduced to calculate the direct radiation correction value. For diffuse radiation, a sky visibility factor is introduced for correction. Finally, the optimized and corrected total horizontal radiation is obtained by combining shadow masks and horizontal geometric correction coefficients. This invention improves the accuracy of horizontal total radiation correction and is applicable to solar radiation resource assessment.

Description

Technical Field

[0001] This invention relates to the field of solar radiation resource assessment technology, and in particular to a method for optimizing and correcting irradiance based on terrain and mountain shadow discrimination and complex terrain. Background Technology

[0002] Total horizontal radiation is a key input in solar energy resource assessment and surface energy balance calculation. In flat areas, total horizontal radiation can represent the level of regional solar energy resources well; however, in undulating terrain areas, such as ridges and valleys, the aspect and slope of the slope will significantly change the path of sunlight to the surface, resulting in strong spatial differences in the light reception status of different grids at the same time.

[0003] The impact of complex terrain on radiation is not entirely consistent across all radiation values. Direct radiation is closely related to the sun's azimuth, solar altitude, slope normal, and terrain shading, and is strongly affected by mountain shadows. Scattered radiation primarily depends on the size of the visible sky hemisphere at the target location; reflections from surrounding slopes and ground features also contribute some radiation energy to the target grid. Therefore, directly applying a single factor to correct the total horizontal radiation can easily mix different physical mechanisms, leading to distorted correction results.

[0004] Common practices in existing engineering applications include directly subtracting shaded areas from binary shading results, or uniformly multiplying the total radiation by a topographic coefficient (i.e., single-factor correction). While these methods are computationally simple and can quickly remove shaded areas, they cannot simultaneously handle direct shading, diffuse shading, neighboring reflections, and non-Lambertian surface differences. In mountain photovoltaic site selection, canyon topographic radiation assessment, and the preparation of meteorological-driven data for complex high-altitude terrain, a correction method is needed that can handle radiation components and topographic mechanisms separately to correct radiation values ​​to approximate the true photometric values. Existing technologies face several problems when applied in this context:

[0005] 1. The total radiation is corrected as a whole without distinguishing between direct and diffuse radiation. The direct component is controlled by both solar geometry and slope geometry, while the diffuse component is controlled by the degree of sky openness. If the same correction factor is used for both, it will weaken the physical interpretability of the model.

[0006] 2. There is a lack of connection between shadow discrimination and radiometric correction. Simply outputting the shadow map cannot complete the radiometric correction, while using only the continuous illumination factor may retain outliers in the backlit shadow area. It is necessary to incorporate Hillshade illumination relationship and shadow mask into the calculation. Hillshade refers to the terrain illumination raster calculated based on solar azimuth angle, solar altitude angle, slope and aspect, which is used to reflect the local slope illumination state of the target raster at a given time. Shadow mask refers to the illumination area discrimination formed by assigning 1 to areas with a Hillshade or normalized slope incident cosine greater than zero, and assigning 0 or invalid values ​​to other areas.

[0007] 3. When only GHI (Global Horizontal Irradiance) data is available, existing methods often cannot establish a differential correction process for direct and scattered radiation, resulting in limited usable meteorological data.

[0008] 4. The surrounding terrain reflection, surface albedo, and non-Lambertian characteristics are not considered uniformly, which can easily lead to systematic biases under different underlying surface conditions such as snow cover, bare soil, or vegetation cover.

[0009] 5. Some slope radiation correction methods may change the physical meaning of the data on the horizontal plane, and additional conversion is still required when using them for meteorological models, photovoltaic resource assessment or spatial interpolation.

[0010] In summary, existing methods for calculating total horizontal radiation data for complex terrain suffer from the problem of simplistic calculations leading to significant errors in the results. Summary of the Invention

[0011] The purpose of this invention is to provide an irradiance optimization and correction method based on terrain and mountain shadow discrimination and complex terrain, so as to solve the problem that the existing calculation of total horizontal radiation data for complex terrain is simple and results in large errors.

[0012] The technical solution adopted by this invention to solve the above-mentioned technical problems provides a method for irradiance optimization and correction based on terrain and mountain shadow discrimination and complex terrain, including:

[0013] Acquire the Digital Elevation Model (DEM), the latitude and longitude of the target area, the time information required for calculation, the Hillshade output path, the SVF (Sky View Factor) raster, the Vf (Vicinity Openness factor) raster, the albedo parameter, the Minnaert non-Lambertian correction factor k, and the radiative data.

[0014] Generate and use Hillshade raster as the output spatial reference, and perform coordinate system, resolution, width and height uniformity processing on SVF raster, Vf raster and radiometric data;

[0015] Calculate the corresponding solar altitude angle, solar zenith angle, and solar azimuth angle based on the required time information and the latitude and longitude of the target area.

[0016] Based on the digital elevation model, the east-west and north-south elevation gradients are calculated using third-order neighborhood difference to obtain slope and aspect data. Combined with the solar zenith angle and solar azimuth angle, the slope incident light relationship of each grid is calculated. The Hillshade grid is filled according to the Hillshade output path and a shadow mask is generated.

[0017] For any target grid corresponding to the target region in the Hillshade grid, the horizon occlusion angle of the target grid in multiple directions is calculated, and the cosine square average value is calculated based on the horizon occlusion angle in each direction and used as the SVF to fill the SVF grid.

[0018] For any target grid corresponding to the target region in the Hillshade grid, the Vf grid is filled by multi-directional horizon angle calculation under finite search distance and finite number of directions.

[0019] The direct normal radiation, scattered radiation, and total horizontal radiation before correction were obtained from the radiation data.

[0020] The direct correction value of the direct component is calculated based on the direct normal radiation, the total horizontal radiation before correction, the albedo parameter, SVF, Vf, and the Minnaert non-Lambertian correction factor k.

[0021] The scattering correction value of the scattering component is calculated based on the scattered radiation and SVF.

[0022] The direct radiation correction value and the scattering correction value are added together, and combined with the shadow mask and the horizontal plane geometric correction coefficient, to calculate the optimized and corrected total horizontal plane radiation.

[0023] In some embodiments, to provide feasible radiation data, the radiation data is the total horizontal radiation before correction;

[0024] Or the radiation data are the total horizontal radiation and direct normal radiation before correction;

[0025] Or the radiation data may be direct normal radiation and scattered radiation.

[0026] In some embodiments, to provide a feasible method for obtaining direct normal radiation, scattered radiation, and total horizontal radiation before correction from radiation data, the method of obtaining direct normal radiation, scattered radiation, and total horizontal radiation before correction from radiation data includes:

[0027] When the radiation data is the total horizontal radiation before correction, the Erbs model is used to perform empirical decomposition on the total horizontal radiation before correction to obtain the estimated values ​​of direct normal radiation and scattered radiation, which are respectively used as direct normal radiation and scattered radiation.

[0028] When the radiation data is the total horizontal radiation and direct normal radiation before correction, or when the radiation data is direct normal radiation and scattered radiation, the scattered radiation is estimated based on the total horizontal radiation and direct normal radiation, or the total horizontal radiation before correction is calculated based on the direct normal radiation and scattered radiation. The calculation formula includes:

[0029] ;

[0030] Where GHIbase refers to the total horizontal radiation before correction, DNI refers to the direct normal radiation, h refers to the solar altitude angle, and DHI refers to the diffuse radiation.

[0031] In some embodiments, to provide a feasible method for calculating solar altitude angle, solar zenith angle, and solar azimuth angle, the method for calculating the solar altitude angle, solar zenith angle, and solar azimuth angle based on the required time information and the latitude and longitude of the target area includes:

[0032] The corresponding solar declination and solar hour angle are calculated based on the required time information and the latitude and longitude of the target area.

[0033] The solar altitude angle is calculated based on the latitude of the target area, the corresponding solar declination, and the solar hour angle. The calculation formula includes:

[0034] ;

[0035] Where h refers to the solar altitude angle, φ refers to the latitude of the target area, δ is the corresponding solar declination, and ω is the corresponding solar hour angle;

[0036] The solar zenith angle is calculated from the solar altitude angle, and the calculation formula includes:

[0037] ;

[0038] Where θz refers to the solar zenith angle;

[0039] The solar altitude angle is calculated using the arctangent function based on the solar hour angle, the latitude of the target region, and the solar declination.

[0040] In some embodiments, to provide a feasible method for calculating the slope incident lighting relationship of each grid cell and filling the Hillshade grid according to the Hillshade output path, the calculation formula for calculating the slope incident lighting relationship of each grid cell and filling the Hillshade grid according to the Hillshade output path includes:

[0041] ;

[0042] Where H refers to the relative light intensity of the target grid under the current sun position, θz refers to the solar zenith angle, α refers to the slope aspect, β refers to the slope, and A refers to the solar azimuth angle.

[0043] In some embodiments, to provide a feasible method for calculating the horizon occlusion angles of any target raster corresponding to a target region in a Hillshade raster in multiple directions, and calculating the cosine squared average value as the SVF based on the horizon occlusion angles in each direction, the method of calculating the horizon occlusion angles of any target raster corresponding to a target region in a Hillshade raster and calculating the cosine squared average value as the SVF based on the horizon occlusion angles in each direction includes:

[0044] For any target grid p corresponding to the target region in the Hillshade grid, the elevation of adjacent grids is read step by step along the search distance d in the i-th azimuth direction, and the candidate elevation angle γi,d is calculated. The calculation formula includes:

[0045] ;

[0046] Where Zi,d refers to the elevation value of the target grid p when it searches to a distance d along the i-th azimuth direction, Zp refers to the elevation of the target grid p, and r refers to the grid resolution;

[0047] Take the maximum value among all candidate elevation angles in the i-th direction as the maximum horizon occlusion angle γi in the i-th direction;

[0048] The cosine squared average value is calculated as the SVF based on the horizon occlusion angle in each direction. The calculation formula includes:

[0049] ;

[0050] Where N refers to the total number of directions.

[0051] In some embodiments, to provide a feasible method for calculating the direct correction value of the direct component based on direct normal radiation, total horizontal surface radiation before correction, albedo parameter, SVF, Vf, and Minnaert non-Lambertian correction factor k, the calculation of the direct correction value of the direct component based on direct normal radiation, total horizontal surface radiation before correction, albedo parameter, SVF, Vf, and Minnaert non-Lambertian correction factor k includes:

[0052] The reflection compensation term is calculated based on the albedo parameter, SVF, and the total horizontal radiation before correction. The calculation formula includes:

[0053] ;

[0054] Where Iref refers to the reflection compensation term, and ρ refers to the albedo parameter;

[0055] The Minnaert correction factor is calculated using the following formula:

[0056] ;

[0057] Where M refers to the Minnaert correction factor, k refers to the Minnaert non-Lambert correction coefficient k, and C refers to the illumination factor obtained by Hillshade or slope incidence cosine normalization.

[0058] The direct radiation correction value for the direct radiation component is calculated using the following formula:

[0059] ;

[0060] Among them, Idir corr It refers to the direct correction value of the direct component, while DNI refers to direct normal radiation.

[0061] In some embodiments, to provide a feasible method for calculating the scattering correction value of the scattering component based on scattered radiation and SVF, the calculation formula for calculating the scattering correction value of the scattering component based on scattered radiation and SVF includes:

[0062] ;

[0063] Among them, Idif corr It refers to the scattering correction value of the scattering component, and DHI refers to the scattered radiation.

[0064] In some embodiments, to provide a feasible method for adding the direct illumination correction value and the scattering correction value, and combining them with a shadow mask and a horizontal plane geometric correction coefficient to calculate the optimized corrected total horizontal plane radiation, the calculation formula for adding the direct illumination correction value and the scattering correction value, and combining them with a shadow mask and a horizontal plane geometric correction coefficient to calculate the optimized corrected total horizontal plane radiation includes:

[0065] ;

[0066] Where GHIcorr refers to the optimized and corrected total horizontal radiation, and Idir corr This refers to the direct correction value of the direct component, Idif corr This refers to the scattering correction value for the scattering component, and shadow_mask refers to the shadow mask. θz refers to the solar zenith angle, and θi refers to the angle of incidence of sunlight relative to the plane to be corrected.

[0067] In some embodiments, for ease of viewing and analysis, after the optimized and corrected total horizontal radiation is calculated, when outputting the results, a non-negative constraint is applied to the optimized and corrected total horizontal radiation corresponding to the valid area of ​​the input data, and a uniform invalid value is assigned to the invalid area of ​​the input data, the invalid area of ​​Hillshade, or any raster that is missing a necessary factor, and the results are output in raster format.

[0068] The beneficial effects of this invention are that, in the scheme of this invention, the direct and scattered components are processed with different topographic factors. The direct component is corrected using Vf, albedo parameters, and Minnaert non-Lambertian correction coefficient k, while the scattered component is corrected using SVF. This improves the accuracy of the total horizontal radiation correction and can directly output the optimized and corrected total horizontal radiation. It can be directly used for photovoltaic resource assessment, meteorological models, spatial interpolation, and eco-hydrological simulation. It also has engineering grid calculation capabilities and can run automatically based on DEM, radiation grid, and related parameters, making it suitable for batch calculations over large areas. Attached Figure Description

[0069] Figure 1 This is a schematic flowchart illustrating the irradiance optimization and correction method based on terrain and mountain shadow discrimination and complex terrain in an embodiment of the present invention. Detailed Implementation

[0070] The technical solution of the present invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0071] like Figure 1 As shown in the embodiment of the present invention, the irradiance optimization and correction method based on terrain and mountain shadow discrimination and complex terrain includes:

[0072] Acquire the Digital Elevation Model (DEM), the latitude and longitude of the target area, the time information required for calculation, the Hillshade output path, the SVF (Sky View Factor) raster, the Vf (Vicinity Openness factor) raster, the albedo parameter, the Minnaert non-Lambertian correction factor k, and the radiative data.

[0073] Generate and use Hillshade raster as the output spatial reference, and perform coordinate system, resolution, width and height uniformity processing on SVF raster, Vf raster and radiometric data;

[0074] Calculate the corresponding solar altitude angle, solar zenith angle, and solar azimuth angle based on the required time information and the latitude and longitude of the target area.

[0075] Based on the digital elevation model, the east-west and north-south elevation gradients are calculated using third-order neighborhood difference to obtain slope and aspect data. Combined with the solar zenith angle and solar azimuth angle, the slope incident light relationship of each grid is calculated. The Hillshade grid is filled according to the Hillshade output path and a shadow mask is generated.

[0076] For any target grid corresponding to the target region in the Hillshade grid, the horizon occlusion angle of the target grid in multiple directions is calculated, and the cosine square average value is calculated based on the horizon occlusion angle in each direction and used as the SVF to fill the SVF grid.

[0077] For any target grid corresponding to the target region in the Hillshade grid, the Vf grid is filled by multi-directional horizon angle calculation under finite search distance and finite number of directions.

[0078] The direct normal radiation, scattered radiation, and total horizontal radiation before correction were obtained from the radiation data.

[0079] The direct correction value of the direct component is calculated based on the direct normal radiation, the total horizontal radiation before correction, the albedo parameter, SVF, Vf, and the Minnaert non-Lambertian correction factor k.

[0080] The scattering correction value of the scattering component is calculated based on the scattered radiation and SVF.

[0081] The direct radiation correction value and the scattering correction value are added together, and combined with the shadow mask and the horizontal plane geometric correction coefficient, to calculate the optimized and corrected total horizontal plane radiation.

[0082] Understandably, in the above embodiments, in complex terrain areas, the mountain shadows (forming Hillshade grids) and lighting conditions (solar altitude angle, solar zenith angle, and solar azimuth angle) are first identified. Then, based on data availability, the direct normal radiation and scattered radiation are obtained. The direct normal radiation is calculated by introducing terrain visibility (SVF and Vf), albedo parameters, and Minnaert non-Lambert correction coefficient k. The scattered radiation is corrected by introducing sky visibility factor (SVF). Finally, the optimized and corrected total horizontal radiation is obtained by combining the shadow mask and the horizontal plane geometric correction coefficient, which improves the accuracy of optimizing and correcting the total horizontal radiation.

[0083] Here, total horizontal radiation refers to the total solar radiation reaching the horizontal plane, composed of the contribution of the direct component projected onto the horizontal plane and the scattered component. Direct normal radiation refers to direct solar radiation perpendicular to the direction of sunlight; when synthesizing the direct contribution to the horizontal plane, it needs to be converted in conjunction with the solar altitude angle or zenith angle. Scattered radiation refers to the solar radiation component that reaches the horizontal plane after being scattered by the atmosphere; its influence by terrain is mainly manifested in the variation of the visible sky range. Hillshade refers to the topographic illumination grid calculated based on the solar azimuth angle, solar altitude angle, slope, and aspect, used to reflect the local slope illumination state of the target grid at a specific time. SVF, or Sky Visibility Factor, specifically refers to the proportion of the sky hemisphere visible to the target grid, with a value ranging from 0 to 1, used for terrain occlusion correction of scattered radiation. Vf, or the nearest terrain openness constraint factor, specifically refers to the factor used to characterize the openness of the target grid's nearest terrain and the impact of local terrain shading on the effective contribution of the direct radiation component. Its value ranges from 0 to 1, and this factor serves as an auxiliary constraint term for terrain correction of the direct radiation component. The introduction of Vf ensures that the direct radiation correction considers not only direct illumination from the direction of solar incidence but also the supplementary radiation effect of nearby terrain on the target point corresponding to the target grid under visible conditions, thus improving the accuracy of the final result.

[0084] When using the Hillshade raster as the output spatial reference, and performing coordinate system, resolution, width, and height consistency processing on SVF raster, Vf raster, and radiometric data, a bilinear resampling method can be used to complete raster registration, ensuring that each factor corresponds pixel by pixel.

[0085] Among them, the radiation data can be the total horizontal radiation before correction;

[0086] It can also be the total horizontal radiation and the direct normal radiation before correction;

[0087] It can also be direct normal radiation and scattered radiation.

[0088] It is understood that the radiation data in the above embodiments can be of various types, and there is no limitation on them here. As can be seen from the following description, as long as the direct normal radiation, scattered radiation and the total horizontal radiation before correction can be obtained from the radiation data, it is acceptable.

[0089] Correspondingly, the direct normal radiation, scattered radiation, and total horizontal radiation before correction can be obtained from radiation data, which may include:

[0090] When the radiation data is the total horizontal radiation before correction, the Erbs model is used to perform an empirical decomposition on the total horizontal radiation before correction to obtain the estimated values ​​of direct normal radiation and scattered radiation, which are used as direct normal radiation and scattered radiation, respectively.

[0091] When the radiation data is the total horizontal radiation and direct normal radiation before correction, or when the radiation data is direct normal radiation and scattered radiation, the scattered radiation is estimated based on the total horizontal radiation and direct normal radiation, or the total horizontal radiation before correction is calculated based on the direct normal radiation and scattered radiation. The calculation formula includes:

[0092] ;

[0093] Where GHIbase refers to the total horizontal radiation before correction, DNI refers to the direct normal radiation, h refers to the solar altitude angle, and DHI refers to the diffuse radiation.

[0094] It is understood that, according to the above embodiments, the irradiance optimization and correction method based on terrain and mountain shadow discrimination and complex terrain can be adapted to a variety of radiation data, enabling it to be applied under different data conditions and improving its versatility.

[0095] Here, when using the Erbs model to empirically decompose the total horizontal radiation before correction and obtain estimates of direct normal radiation and scattered radiation, the Erbs model is used to calculate the scattering ratio kd, and then... and The direct normal radiation DNI and the scattered radiation DHI were estimated.

[0096] In some embodiments, to provide a feasible method for calculating the solar altitude angle, solar zenith angle, and solar azimuth angle, the method for calculating the solar altitude angle, solar zenith angle, and solar azimuth angle based on the required time information and the latitude and longitude of the target area may include:

[0097] The corresponding solar declination and solar hour angle are calculated based on the required time information and the latitude and longitude of the target area.

[0098] The solar altitude angle is calculated based on the latitude of the target area, the corresponding solar declination, and the solar hour angle. The calculation formula includes:

[0099] ;

[0100] Where h refers to the solar altitude angle, φ refers to the latitude of the target area, δ is the corresponding solar declination, and ω is the corresponding solar hour angle;

[0101] The solar zenith angle is calculated from the solar altitude angle, and the calculation formula includes:

[0102] ;

[0103] Where θz refers to the solar zenith angle;

[0104] The solar altitude angle is calculated using the arctangent function based on the solar hour angle, the latitude of the target region, and the solar declination.

[0105] Understandably, this step provides a unified solar geometry for calculating the mountain shadow and separating the direct and scattered components in the total horizontal radiation before correction.

[0106] In some embodiments, to provide a feasible method for calculating the slope incident lighting relationship of each grid cell and filling the Hillshade grid according to the Hillshade output path, the calculation formula may include:

[0107] ;

[0108] Where H refers to the relative light intensity of the target grid under the current sun position, θz refers to the solar zenith angle, α refers to the slope aspect, β refers to the slope, and A refers to the solar azimuth angle.

[0109] When H is less than or equal to 0, it indicates that the target grid is under backlight or shadow conditions; when H is greater than 0, it indicates that the grid is under direct lighting conditions. Mountain shadows are reclassified based on occlusion (corresponding to backlight or shadow conditions) and non-occlusion (corresponding to direct lighting conditions).

[0110] The shadow mask refers to the identification of the illuminated area by assigning a value of 1 to areas where the Hillshade or normalized slope incident cosine is greater than zero, and assigning a value of 0 or invalid value to the remaining areas. The corresponding shadow mask is: marking the target grid with direct illumination conditions as 1, and marking the target grid under backlight or shadow conditions as 0 or invalid value.

[0111] In some embodiments, to provide a feasible method for calculating the horizon occlusion angles of any target raster corresponding to a target region in a Hillshade raster in multiple directions, and calculating the cosine squared average as the SVF based on the horizon occlusion angles in each direction, the method for calculating the horizon occlusion angles of any target raster corresponding to a target region in a Hillshade raster and calculating the cosine squared average as the SVF based on the horizon occlusion angles in each direction may include:

[0112] For any target grid p corresponding to the target region in the Hillshade grid, the elevation of adjacent grids is read step by step along the search distance d in the i-th azimuth direction, and the candidate elevation angle γi,d is calculated. The calculation formula may include:

[0113] ;

[0114] Where Zi,d refers to the elevation value of the target grid p when it searches to a distance d along the i-th azimuth direction, Zp refers to the elevation of the target grid p, and r refers to the grid resolution;

[0115] Take the maximum value among all candidate elevation angles in the i-th direction as the maximum horizon occlusion angle γi in the i-th direction;

[0116] The cosine squared average value is calculated based on the horizon occlusion angle in each direction as the SVF. The calculation formula may include:

[0117] ;

[0118] Where N refers to the total number of directions.

[0119] Here, a total of 32 directional data points can be selected. The maximum search distance can be set according to the resolution of the digital elevation model and the terrain scale of the target area, with a default value of 1000m. Since the SVF value is limited to between 0 and 1, a larger value indicates a higher degree of sky openness and weaker terrain obstruction of scattered radiation.

[0120] For any target grid cell corresponding to the target region in the Hillshade grid, the multi-directional horizon angle Vf is calculated using a finite search distance and a finite number of directions. The calculation of Vf employs a similar framework to that used for SVF, but its parameter range should differ. A shorter maximum search distance and fewer directions are preferred to characterize the openness of the target grid's neighboring terrain. In Vf, the horizon angle can be denoted as γi,Vf, which refers to the maximum horizon occlusion angle in the i-th direction within the neighboring search range. SVF emphasizes the openness of the sky hemisphere, while Vf emphasizes the constraints of neighboring terrain openness.

[0121] In some embodiments, to provide a feasible method for calculating the direct correction value of the direct component based on direct normal radiation, total horizontal surface radiation before correction, albedo parameter, SVF, Vf, and Minnaert non-Lambertian correction factor k, the calculation of the direct correction value of the direct component based on direct normal radiation, total horizontal surface radiation before correction, albedo parameter, SVF, Vf, and Minnaert non-Lambertian correction factor k may include:

[0122] The reflection compensation term is calculated based on the albedo parameter, SVF, and the total horizontal radiation before correction. The calculation formula includes:

[0123] ;

[0124] Where Iref refers to the reflection compensation term, and ρ refers to the albedo parameter;

[0125] The Minnaert correction factor is calculated using the following formula:

[0126] ;

[0127] Where M refers to the Minnaert correction factor, k refers to the Minnaert non-Lambert correction coefficient k, and C refers to the illumination factor obtained by Hillshade or slope incidence cosine normalization.

[0128] The direct radiation correction value for the direct radiation component is calculated using the following formula:

[0129] ;

[0130] Among them, Idir corr It refers to the direct correction value of the direct component, while DNI refers to direct normal radiation.

[0131] As can be seen, this formula simultaneously reflects the direct incidence terrain visibility constraint (corresponding to Vf), adjacent terrain reflection compensation (corresponding to Iref), and corrections for different surface radiation characteristics (corresponding to M). The specific method for obtaining C is as follows: if the values ​​in the Hillshade raster are stored as grayscale values ​​from 0 to 255, then C = clip(Hillshade / 255, 0, 1); if the original value of the slope incidence cosine H is directly output, then C = clip(H, 0, 1).

[0132] Here, when C≤0, the target grid does not have effective direct light reception conditions, so M can be set to 0, so that the direct light correction value of the subsequently calculated direct light component is 0.

[0133] The Minnaert non-Lambertian correction factor k is determined based on the surface type, specifically as follows: when the target raster's surface type is water or smooth, the k value ranges from 0.95 to 1, preferably 1; when the target raster's surface type is snow, the k value ranges from 0.8 to 1, preferably 0.9; when the target raster's surface type is bare soil or sand, the k value ranges from 0.7 to 0.85, preferably 0.8; when the target raster's surface type is grassland or shrubs, the k value ranges from 0.6 to 0.75, preferably 0.7; and when the target raster's surface type is forest or dense vegetation, the k value ranges from 0.5 to 0.7, preferably 0.6.

[0134] In some embodiments, to provide a feasible method for calculating the scattering correction value of the scattering component based on scattered radiation and SVF, the calculation formula for calculating the scattering correction value of the scattering component based on scattered radiation and SVF may include:

[0135] ;

[0136] Among them, Idif corr It refers to the scattering correction value of the scattering component, and DHI refers to the scattered radiation.

[0137] Here, the smaller the SVF, the more severely the sky hemisphere is blocked by the surrounding terrain, and the lower the proportion of scattered radiation retained.

[0138] In some embodiments, to provide a feasible method for adding the direct illumination correction value and the scattering correction value, and combining the shadow mask and the horizontal plane geometric correction coefficient to calculate the optimized corrected total horizontal plane radiation, the calculation formula may include:

[0139] ;

[0140] Where GHIcorr refers to the optimized and corrected total horizontal radiation, and Idir corr This refers to the direct correction value of the direct component, Idif corr This refers to the scattering correction value for the scattering component, and shadow_mask refers to the shadow mask. θz refers to the solar zenith angle, and θi refers to the angle of incidence of sunlight relative to the plane to be corrected.

[0141] Here, in the horizontal total radiation correction scenario of this application, since the plane to be corrected is horizontal, the incident angle θi is equal to the solar zenith angle θz, i.e., K=1. This term is reserved as an interface for subsequent extension to other incident geometries, such as non-horizontal planes and other incident geometries.

[0142] In some embodiments, for ease of viewing and analysis, after the optimized and corrected total horizontal radiation is calculated, when outputting the results, a non-negative constraint is applied to the optimized and corrected total horizontal radiation corresponding to the valid area of ​​the input data, and a uniform invalid value is assigned to the invalid area of ​​the input data, the invalid area of ​​Hillshade, or any raster that is missing a necessary factor, and the results are output in raster format.

[0143] Here, the valid input data area refers to the area where all necessary rasters involved in the calculation have valid values ​​at the same target cell and are not areas without data; the invalid input data area refers to the area where any necessary raster involved in the calculation has a null value or no data at the target cell; the invalid Hillshade area refers to the area where the Hillshade raster has a null value or no data at the target cell, or the corresponding original cell of the digital elevation model has invalid elevation. Necessary factors may include direct normal radiation, diffuse radiation (which can be the decomposed direct normal radiation and / or diffuse radiation), SVF, Vf, Hillshade, shadow mask, albedo parameter, Minnaert non-Lambertian correction factor k, and solar geometric parameters (including solar altitude angle, solar zenith angle, and solar azimuth angle). Invalid values ​​can be -9999, and the raster format can be GeoTIFF, etc.

Claims

1. A method for optimizing and correcting irradiance based on terrain and mountain shadow discrimination and complex terrain, characterized in that, include: Acquire digital elevation model, target area latitude and longitude, required calculation time information, Hillshade output path, SVF raster, Vf raster, albedo parameter, Minnaert non-Lambert correction factor k, and radiometric data; Generate and use Hillshade raster as the output spatial reference, and perform coordinate system, resolution, width and height uniformity processing on SVF raster, Vf raster and radiometric data; Calculate the corresponding solar altitude angle, solar zenith angle, and solar azimuth angle based on the required time information and the latitude and longitude of the target area. Based on the digital elevation model, the east-west and north-south elevation gradients are calculated using third-order neighborhood difference to obtain slope and aspect data. Combined with the solar zenith angle and solar azimuth angle, the slope incident light relationship of each grid is calculated. The Hillshade grid is filled according to the Hillshade output path and a shadow mask is generated. For any target grid corresponding to the target region in the Hillshade grid, the horizon occlusion angle of the target grid in multiple directions is calculated, and the cosine square average value is calculated based on the horizon occlusion angle in each direction and used as the SVF to fill the SVF grid. For any target grid corresponding to the target region in the Hillshade grid, the Vf grid is filled by multi-directional horizon angle calculation under finite search distance and finite number of directions. The direct normal radiation, scattered radiation, and total horizontal radiation before correction were obtained from the radiation data. The direct correction value of the direct component is calculated based on the direct normal radiation, the total horizontal radiation before correction, the albedo parameter, SVF, Vf, and the Minnaert non-Lambertian correction factor k. The scattering correction value of the scattering component is calculated based on the scattered radiation and SVF. The direct radiation correction value and the scattering correction value are added together, and combined with the shadow mask and the horizontal plane geometric correction coefficient, to calculate the optimized and corrected total horizontal plane radiation.

2. The irradiance optimization and correction method based on terrain and mountain shadow discrimination and complex terrain as described in claim 1, characterized in that, The radiation data mentioned is the total horizontal radiation before correction; Or the radiation data are the total horizontal radiation and direct normal radiation before correction; Or the radiation data may be direct normal radiation and scattered radiation.

3. The irradiance optimization and correction method based on terrain and mountain shadow discrimination and complex terrain as described in claim 2, characterized in that, The process of obtaining direct normal radiation, scattered radiation, and total horizontal radiation before correction based on radiation data includes: When the radiation data is the total horizontal radiation before correction, the Erbs model is used to perform empirical decomposition on the total horizontal radiation before correction to obtain the estimated values ​​of direct normal radiation and scattered radiation, which are respectively used as direct normal radiation and scattered radiation. When the radiation data is the total horizontal radiation and direct normal radiation before correction, or when the radiation data is direct normal radiation and scattered radiation, the scattered radiation is estimated based on the total horizontal radiation and direct normal radiation, or the total horizontal radiation before correction is calculated based on the direct normal radiation and scattered radiation. The calculation formula includes: ； Where GHIbase refers to the total horizontal radiation before correction, DNI refers to the direct normal radiation, h refers to the solar altitude angle, and DHI refers to the diffuse radiation.

4. The irradiance optimization and correction method based on terrain and mountain shadow discrimination and complex terrain as described in claim 1, characterized in that, The method for calculating the solar altitude angle, solar zenith angle, and solar azimuth angle based on the required time information and the latitude and longitude of the target area includes: The corresponding solar declination and solar hour angle are calculated based on the required time information and the latitude and longitude of the target area. The solar altitude angle is calculated based on the latitude of the target area, the corresponding solar declination, and the solar hour angle. The calculation formula includes: ； Where h refers to the solar altitude angle, φ refers to the latitude of the target area, δ is the corresponding solar declination, and ω is the corresponding solar hour angle; The solar zenith angle is calculated from the solar altitude angle, and the calculation formula includes: ； Where θz refers to the solar zenith angle; The solar altitude angle is calculated using the arctangent function based on the solar hour angle, the latitude of the target region, and the solar declination.

5. The irradiance optimization and correction method based on terrain and mountain shadow discrimination and complex terrain as described in claim 1, characterized in that, The slope incident light relationship of each grid cell is calculated and filled into the Hillshade grid according to the Hillshade output path. The calculation formula includes: ； Where H refers to the relative light intensity of the target grid under the current sun position, θz refers to the solar zenith angle, α refers to the slope aspect, β refers to the slope, and A refers to the solar azimuth angle.

6. The irradiance optimization and correction method based on terrain and mountain shadow discrimination and complex terrain as described in claim 1, characterized in that, For any target grid cell corresponding to the target region in the Hillshade grid, the horizon occlusion angle of the target grid cell in multiple directions is calculated, and the cosine squared average value is calculated as the SVF based on the horizon occlusion angle in each direction, including: For any target grid p corresponding to the target region in the Hillshade grid, the elevation of adjacent grids is read step by step along the search distance d in the i-th azimuth direction, and the candidate elevation angle γi,d is calculated. The calculation formula includes: ； Where Zi,d refers to the elevation value of the target grid p when it searches to a distance d along the i-th azimuth direction, Zp refers to the elevation of the target grid p, and r refers to the grid resolution; Take the maximum value among all candidate elevation angles in the i-th direction as the maximum horizon occlusion angle γi in the i-th direction; The cosine squared average value is calculated as the SVF based on the horizon occlusion angle in each direction. The calculation formula includes: ； Where N refers to the total number of directions.

7. The irradiance optimization and correction method based on terrain and mountain shadow discrimination and complex terrain as described in claim 1, characterized in that, The direct correction value for the direct component, calculated based on direct normal radiation, total horizontal radiation before correction, albedo parameter, SVF, Vf, and Minnaert non-Lambertian correction factor k, includes: The reflection compensation term is calculated based on the albedo parameter, SVF, and the total horizontal radiation before correction. The calculation formula includes: ； Where Iref refers to the reflection compensation term, and ρ refers to the albedo parameter; The Minnaert correction factor is calculated using the following formula: ； Where M refers to the Minnaert correction factor, k refers to the Minnaert non-Lambert correction coefficient k, and C refers to the illumination factor obtained by Hillshade or slope incidence cosine normalization. The direct radiation correction value for the direct radiation component is calculated using the following formula: ； Among them, Idir corr It refers to the direct correction value of the direct component, while DNI refers to direct normal radiation.

8. The irradiance optimization and correction method based on terrain and mountain shadow discrimination and complex terrain as described in claim 1, characterized in that, The scattering correction value for the scattering component is calculated based on the scattered radiation and SVF, and the calculation formula includes: ； Among them, Idif corr It refers to the scattering correction value of the scattering component, and DHI refers to the scattered radiation.

9. The irradiance optimization and correction method based on terrain and mountain shadow discrimination and complex terrain as described in claim 1, characterized in that, The direct radiation correction value and the scattering correction value are added together, and combined with the shadow mask and the horizontal plane geometric correction coefficient, to calculate the optimized and corrected total horizontal plane radiation. The calculation formula is as follows: include: ； Where GHIcorr refers to the optimized and corrected total horizontal radiation, and Idir corr This refers to the direct correction value of the direct component, Idif corr This refers to the scattering correction value for the scattering component, and shadow_mask refers to the shadow mask. θz refers to the solar zenith angle, and θi refers to the angle of incidence of sunlight relative to the plane to be corrected.

10. The irradiance optimization and correction method based on terrain and mountain shadow discrimination and complex terrain as described in any one of claims 1-9, characterized in that, After the optimized and corrected total horizontal radiation is obtained from the calculation, when outputting the results, a non-negative constraint is applied to the optimized and corrected total horizontal radiation corresponding to the valid area of ​​the input data. A uniform invalid value is assigned to the invalid area of ​​the input data, the invalid area of ​​Hillshade, or any raster that is missing a necessary factor, and the results are output in raster format.

Images