Building global sunshine analysis method based on live-action three-dimensional model

Abstract

The invention discloses a building global sunshine analysis method based on a live-action three-dimensional model, and the method comprises the steps: calculating a sunshine vector according to a time window, the number of analysis days, a moment list and the latitude and longitude of a location, and carrying out the loading, normal vector correction, roof recognition and adaptive subdivision of a building model. And for the subdivided building model, constructing a ray parallel to the sunlight direction by taking the center point of the surface as a starting point, judging whether the surface is shielded or not through ray-grid intersection detection, and calculating the sunshine duration of the triangular surface. And further calculating direct radiation, scattered radiation and ground reflection radiation of the surface, and estimating the total irradiance of the triangular surface. And calculating a mean value of all results on the time sequence to obtain a comprehensive sunlight analysis result. And matching an analysis result with a model surface index, exporting a data table, and embedding a surface patch color into a model for storage.

Description

Technical Field

[0001] The present invention relates to the fields of geographic information systems and building environment analysis, and in particular to a global sunlight analysis method for buildings based on a real-scene three-dimensional model. Background Art

[0002] Sunlight analysis refers to the process of analyzing the hours of sunlight and their impact on individual buildings or complexes using manual or computer analysis software. With the acceleration of urbanization and the rise of high-density buildings in my country, ensuring the quality of sunlight has become a major concern for residents. Furthermore, with the rise of distributed photovoltaic power generation, sunlight has become a key factor affecting building environmental quality, living comfort, and energy efficiency.

[0003] In order to conduct sunlight analysis on buildings, there are mainly the following methods: traditional manual sunlight analysis and calculation method, two-dimensional CAD sunlight time calculation method, and three-dimensional sunlight time calculation method.

[0004] Traditional manual methods for calculating sunlight exposure rely on a stick-shadow solar chart. This method uses the fundamental relationship between a stick and its shadow to describe the sun's motion, using the trajectory of the stick's shadow in the sun to simulate the trajectory of the object's projection. The stick-shadow solar chart is superimposed on the design plan to analyze the time period each building is exposed to sunlight. This method is cumbersome, time-consuming, and labor-intensive, and only yields an approximate estimate of sunlight duration.

[0005] The 2D CAD method for calculating sunlight duration, which uses CAD software, is the current mainstream method for sunlight analysis. Its principle is also derived from the stick shadow sun chart. However, this method is based on 2D construction drawings and does not consider the impact of terrain undulations and building roof forms. Alternatively, using Arc GIS to convert a building base map into an elevation raster and then utilizing tools such as a raster calculator, slope and aspect, and hillshade can also serve as a 2D sunlight analysis tool. However, this method also suffers from shortcomings such as an inability to account for detailed roof structures and complex occlusion relationships.

[0006] Three-dimensional sunshine duration calculation methods use 3D design software combined with ray tracing to simulate building sunlight conditions. While this method is more advanced, faster, and more accurate, it also has limitations. Using SketchUp combined with ShadowAnalysis for sunlight analysis can produce intuitive sunshine hour charts, but the result is a two-dimensional raster image from a fixed camera perspective, preventing a three-dimensional global sunlight analysis. Simulating sunlight using ray tracing renderers in 3D modeling software such as 3ds Max, Blender, and Maya also faces the aforementioned issues, resulting in only a single-view two-dimensional raster image and loss of building facet attribute information. Chinese Patent No. 2023115096538 discloses a sunlight analysis method for real-world 3D urban scenes. This method allows for viewing results from different perspectives, but essentially relies on real-time rendering from different perspectives, making it difficult to statically store the results. Furthermore, this method uses stretched results from 2D building shapefile data for analysis, without considering detailed roof structures. Using Rhino in conjunction with Grasshopper and Ladybug can produce 3D global solar radiation analysis results for buildings. However, this requires the building to be enclosed, making it difficult to process low-precision and coarse real-world 3D data. Furthermore, visual programming requires a certain level of familiarity. Existing methods fail to simultaneously account for building details, terrain undulations, and the generation of global solar radiation analysis results. Summary of the Invention

[0007] Purpose of the invention: The purpose of the present invention is to provide a global sunlight analysis method for buildings based on a real-scene three-dimensional model, which takes into account the detailed structure of the building and has adjustable analysis granularity.

[0008] Technical solution: The method for analyzing global sunlight exposure of a building based on a real-scene three-dimensional model of the present invention comprises the following steps:

[0009] (1) Set the time window, analysis days and sampling time list, sample sunshine information, determine the sunshine analysis time series set, and determine the sunlight unit direction vector based on the longitude and latitude of the building target location;

[0010] (2) Load the target building mesh model, identify and separate the roof and facade panels;

[0011] (3) Using different resolutions to triangulate the roof and facade slices, respectively. The resolution of the roof slice is greater than that of the facade slice. After triangulation, the roof and facade slices are combined into a model.

[0012] (4) Based on the sunlight unit direction vector obtained in step (1) and the surface center point obtained by triangulating the surface in step (3), construct an emission ray starting from the surface center point and along the sunlight unit direction vector; perform batch emission ray-grid intersection detection, using "0" to mark the obscured surface and "1" to mark the unobscured surface;

[0013] (5) Considering the occlusion of the patch, calculate the total irradiance of the patch;

[0014] (6) Calculate the sum of the shielding marks and the average irradiance at all times of the day to obtain the sunshine hours and average irradiance of the day; take the average of the analysis results of different days to obtain the average sunshine hours and average irradiance;

[0015] (7) Output the sunshine analysis result table, map the sunshine analysis results to the color spectrum, embed the model as the color information of the surface, and export the model containing the sunshine analysis results.

[0016] Furthermore, in step (1), the longitude and latitude of the building target geographical location are combined to determine the unit direction vector of sunlight, as follows:

[0017] Determine the sunlight azimuth and altitude angle based on the sampling time and the longitude and latitude of the building's target location, and calculate the sunlight unit direction vector S:

[0018]

[0019] Among them, α is the solar altitude angle, which ranges from θ is the solar azimuth angle, ranging from [0, 2π); S x 、S y and S z They represent the components of the sunlight unit direction vector in the X direction (east-west direction), Y direction (north-south direction), and Z direction (vertical direction).

[0020] Furthermore, step (2) is as follows:

[0021] Calculate the angle θ between the normal vector of the triangle and the vertical direction:

[0022]

[0023] Wherein, the range of θ is [0°, 180°]; is the normal vector of the triangle; is the vertical unit vector;

[0024] Set a roof inclination angle threshold and mark the triangles whose calculation results are less than the roof inclination angle threshold as candidate roof faces; the others are facade faces;

[0025] Taking the center point of each candidate roof triangle as the starting point and the vertical upward direction as the ray direction, a ray list is constructed; ray and grid intersection detection is performed, the number of hits for each ray is counted, and the triangles corresponding to rays with a hit number greater than 0 are deleted from the candidate roof patches to obtain the final roof patch.

[0026] Furthermore, in step (2), before identifying and separating the roof face and the facade face, the normal vector of the triangle face is corrected, and the normal vector of the triangle face with the tilt angle θ∈[150°, 180°] is reversed and adjusted to face the front side upward.

[0027] Furthermore, in step (3), the triangular surface is subdivided as follows:

[0028] Assume that the three vertices of a triangle in three-dimensional space are A, B, and C. Calculate the longest side of the triangle, mark the two vertices where the longest side is located as P1 and P2, and the third vertex as P3; calculate the X-axis direction vector of the target two-dimensional coordinate system

[0029]

[0030] Triangle unit normal vector

[0031]

[0032] Calculate the Y-axis direction vector of the target two-dimensional coordinate system

[0033]

[0034] Determine the coordinate transformation relationship from the three-dimensional space projection to the target two-dimensional plane as follows:

[0035]

[0036] Among them, P is the coordinate of a point in the original three-dimensional coordinate system, P project is the coordinate of the point after projection onto the two-dimensional plane;

[0037] Project the three vertices of the triangle (A, B, and C) to obtain points A2, B2, and C2. Calculate the X and Y coordinate ranges of the projected triangle. Using the resolution of the triangulated surface as the spacing, generate X and Y coordinate sequences, combine them into grid points, and use the barycentric coordinate formula to determine whether the grid point is inside the triangle.

[0038] Add A2, B2, C2, and the mesh points inside the triangle to the vertex list. Calculate the intersection of the three sides of the projected triangle with parallel lines parallel to the Y axis and spaced at the resolution of the triangulated face, and add the intersections to the vertex list. Perform triangulation on the vertices in the vertex list to generate the subdivided faces.

[0039] Furthermore, in step (3), the barycentric coordinate formula is used to determine whether the grid point is inside the triangle, as follows:

[0040]

[0041] u=1-vw

[0042] Among them, u, v, and w are the centroid coordinates of point P on the two-dimensional plane for triangles A2, B2, and C2 respectively. P is the coordinate of a point on the two-dimensional plane. If u, v, w ≥ 0 and u + v + w = ​​1, then P is inside or on the edge of triangles A2, B2, and C2.

[0043] Furthermore, in step (3), the points in the vertex list are converted back to the original three-dimensional coordinate system as follows:

[0044] P3=u·A+v·B+w·C

[0045] Wherein, P3 is the coordinate converted back to the three-dimensional coordinate system.

[0046] Furthermore, in step (4), the surface occlusion judgment method is as follows: based on the sunlight unit direction vector in step (1) and the surface center point obtained by triangulating the surface in step (3), a ray list is constructed with each triangulated surface center point as the starting point and the sunlight direction as the ray direction; batch ray and grid intersection detection is performed, the number of hits for each ray is counted, and the value corresponding to the ray with a hit number greater than 0 is set to "0", indicating that the surface is occluded.

[0047] Furthermore, in step (5), the total irradiance of the patch I total for

[0048] I total =I direct +I diffuse +I reflected

[0049] Among them, I total is the total irradiance of the patch, which means the total electromagnetic radiation power received by the patch per unit area; I direct I is the solar direct radiation power received by the patch per unit area; diffuse I is the atmospheric scattered radiation power received by the patch per unit area; reflectedis the ground reflected radiation power received by the patch per unit area;

[0050] According to the unobstructed surface, calculate I direct :

[0051]

[0052] Where DNI is the irradiance of sunlight directly from the solar disk to the surface orthogonal to the light path; is the patch normal vector of the unoccluded face; is the unit vector of sunlight direction. It is the cosine of the surface normal vector and the direction of sunlight. A negative value means facing away from the sun.

[0053] Based on all faces, calculate I diffuse :

[0054] I diffuse =DHI·F sky

[0055] F sky =(1+cosθ) / 2

[0056] Where DHI is the atmospheric scattered radiation irradiance received on the horizontal plane; F sky is the sky visibility factor, which represents the visibility between the patch and the sky; θ is the difference between the patch normal vector and the vertical unit vector The angle between

[0057] Calculate I for all faces reflected :

[0058] I reflected =ρ·GHI·F ground

[0059] F ground =(1-cosθ) / 2

[0060] Where GHI is the total horizontal irradiance of the sun’s DHI and DNI reaching the horizontal surface; ρ is the ground reflectivity; F ground is the ground visibility factor, which represents the visibility between the patch and the ground; θ is the difference between the patch normal vector and the vertical unit vector Angle.

[0061] Furthermore, in step (6), the output sunshine analysis result table includes the index information of the triangle face "face_index"; the average sunshine hours of the triangle face "illumination", in units of h, with a minimum value of 0 and a maximum value of 13; the average irradiance of the triangle face "irradiance", in units of W / m 2 .

[0062] Beneficial effects: Compared with the prior art, the present invention has the following significant advantages: 1. The present invention takes into account the detailed structure of the building and the analysis granularity is adjustable, realizing different analysis granularities for the roof plane and the facade; 2. The present invention generates global three-dimensional static sunshine analysis results, realizing efficient storage and browsing of the results; 3. The present invention further calculates the irradiance while taking into account the sunshine hours, supporting scenarios such as photovoltaic potential assessment; 4. It has low data requirements and supports rough grid data, topological structure errors, and normal vector errors; 5. The present invention is easy to operate, the program is lightweight, and the threshold for use is low. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 It is a schematic diagram of the process of the present invention;

[0064] Figure 2 A schematic diagram of the original triangle in a three-dimensional coordinate system;

[0065] Figure 3 Schematic diagram of a triangular face projected onto a two-dimensional plane;

[0066] Figure 4 A schematic diagram of all vertices used for subdivision;

[0067] Figure 5 This is a schematic diagram of the subdivision result after triangulation;

[0068] Figure 6 This is a schematic diagram of the result of converting the subdivision result back to a three-dimensional coordinate system;

[0069] Figure 7 This is the rendering of the original model;

[0070] Figure 8 This is the rendering of the model after subdivision;

[0071] Figure 9 This is a schematic diagram of the visualized irradiance results of the sunlight analysis of the test building according to the present invention. DETAILED DESCRIPTION

[0072] The present invention will be further described below with reference to the accompanying drawings.

[0073] The method for analyzing global sunlight exposure of a building based on a real-scene three-dimensional model of the present invention comprises the following steps:

[0074] (1) Set the time window, analysis days and sampling time list to sample sunshine information, determine the sunshine analysis time series set, and determine the unit direction vector of sunlight based on the longitude and latitude of the building target location.

[0075] Set the time window, analysis days, and moment list. The time window is the time range used for sunshine analysis, determined by the start time and end time (year, month, day, hour, minute, second). The analysis days are the days for sunshine analysis. The moment list is a list that stores a series of moments (hour, minute, second), indicating which moments of the day are selected for sunshine analysis. Based on the analysis days, sample dates at equal intervals within the time window to obtain the date information (year, month, day) for sunshine analysis. Concatenate the date information with the moment list to obtain the analysis time series (year, month, day, hour, minute, second). During sunshine analysis, sunshine analysis is performed on the times in the above analysis time series in sequence and the results are integrated.

[0076] Get the sunlight azimuth and altitude information based on time, longitude and latitude. Calculate the sunlight unit direction vector as follows:

[0077]

[0078] Among them, α is the solar altitude angle, which ranges from θ is the solar azimuth angle, ranging from [0, 2π); S x 、S y and S z They represent the components of the sunlight unit direction vector in the X direction (east-west direction), Y direction (north-south direction), and Z direction (vertical direction).

[0079] (2) Load the target building mesh model, identify and separate the roof and facade patches.

[0080] First, calculate the angle between the triangle normal vector and the vertical direction as follows:

[0081]

[0082] Where θ is the angle between the normal vector of the triangle surface and the vertical direction, ranging from [0°, 180°]; is the normal vector of the triangle; is the vertical unit vector.

[0083] If the input model has not been normal-checked and corrected, there may be inconsistencies in normal vectors, which can cause some roof triangles to face up incorrectly, resulting in significant deviations in subsequent solar radiation analysis. Normal correction can be performed to reverse the normal vectors of triangles with large tilt angles, adjusting them to face up. Here, θ∈[150°, 180°] is used. Normal correction can mitigate the effects of model normal errors and ensure that solar radiation analysis results reflect the maximum photovoltaic potential. However, this may cause some triangles that should be facing downward to face the wrong direction, introducing further deviations.

[0084] A roof inclination angle threshold (here: 50°) is set, and triangular faces with calculated results less than the threshold are marked as candidate roof faces. A ray list is constructed, starting from the center point of each candidate roof triangle and pointing vertically upward. Then, batch ray and mesh intersection detection is performed, and the number of hits for each ray is counted. To eliminate the influence of self-intersections of patches, a minimum distance threshold (here: 0.001m) between the hit point and the ray starting point is set to filter out hit points less than the threshold. Triangular faces corresponding to rays with a hit count greater than 0 are excluded from the candidate roof faces to obtain the final roof face.

[0085] (3) The roof and facade patches are triangulated using different resolutions. The resolution of the roof patch is greater than that of the facade patch. The triangulated roof and facade patches are then combined into a model.

[0086] Set the subdivision resolution of the roof surface to 0.05m and the subdivision resolution of the facade to 1m, and subdivide the triangular surface. Assume that the three vertices of the triangle in the three-dimensional space are A, B, and C. Figure 2 Calculate the longest side of the triangle and mark the two vertices where the longest side is located as P1 and P2, and the third vertex as P3.

[0087] Calculate the X-axis direction vector of the target two-dimensional coordinate system Here’s how:

[0088]

[0089] Calculate the triangle unit normal vector Here’s how:

[0090]

[0091] Calculate the Y-axis direction vector of the target two-dimensional coordinate system Here’s how:

[0092]

[0093] Determine the coordinate transformation relationship from the three-dimensional space projection to the target two-dimensional plane as follows:

[0094]

[0095] Among them, P is the coordinate of a point in the original three-dimensional coordinate system, P project are the coordinates of the point after projection onto the two-dimensional plane.

[0096] Project the three vertices of the triangle A, B, and C to obtain points A2, B2, and C2, such as Figure 3Calculate the X and Y ranges of the projected triangle, generate X and Y coordinate sequences based on the specified resolution, and combine them into grid points. Use the barycentric coordinate formula to determine whether a grid point is inside the triangle. The method is as follows:

[0097]

[0098] u=1-vw

[0099] Where P is the coordinate of a point on the two-dimensional plane. If u, v, w ≥ 0 and u + v + w = ​​1, then P is inside or on the edge of triangles A2, B2, and C2.

[0100] Initialize an empty vertex list to store the subdivided vertices. Add A2, B2, C2, and the mesh points inside the triangle to the vertex list. Calculate the Y coordinates of the three sides of the triangle at the X coordinate sequence generated above with the specified resolution as the spacing, that is, the intersection of the three sides of the triangle with the parallel lines parallel to the Y axis and with the specified resolution as the spacing, and add the intersection points to the vertex list. The final content of the vertex list is as follows: Figure 4 Indicated by the red mark in the middle.

[0101] For the vertices in the vertex list, perform Delaunay triangulation to generate the subdivided surface. Delaunay triangulation has the characteristics of empty circle and maximizing the minimum angle, which can generate homogeneous sampling surfaces. This method is based on existing technology and will not be described in detail here. The result after triangulation is as follows Figure 5 shown.

[0102] Convert the points in the vertex list back to the original three-dimensional coordinate system as follows:

[0103] P3=u·A+v·B+w·C

[0104] Among them, u, v, and w are the centroid coordinates of point P on the two-dimensional plane for triangles A2, B2, and C2 respectively; P3 is the coordinate converted back to the three-dimensional coordinate system.

[0105] The result of converting the subdivision result back to the three-dimensional coordinate system is as follows Figure 6 The effect of using this adaptive subdivision method to subdivide a building model is shown in Figure 7 、 Figure 8 As shown, Figure 7 For the original model, Figure 8 The model after segmentation.

[0106] (4) Based on the sunlight unit direction vector of step (1) and the surface center point obtained by triangulating the surface in step (3), construct an emission ray starting from the surface center point and along the sunlight unit direction vector; perform batch emission ray-grid intersection detection, use "0" to mark the obscured surface and "1" to mark the unobscured surface.

[0107] Construct a ray list using the center point of each triangular face obtained in step (3) as the starting point and the direction of sunlight as the ray direction. Then perform batch ray and grid intersection detection and count the number of hits for each ray. To eliminate the influence of self-intersection of the facets, set a minimum distance threshold between the hit point and the ray starting point, which is 0.001m here, and filter out hit points smaller than the threshold. Initialize the occlusion result list, and the list size is the same as the number of faces. Set the value corresponding to the ray with a hit count greater than 0 to "0", indicating that the face is occluded.

[0108] (5) Considering the occlusion of the patch, calculate the total irradiance of the patch.

[0109] Extract the unobstructed patches facing the sun, calculate the angle between the patch normal and the direction of sunlight, and calculate the product of the cosine of the angle and the DNI to obtain the patch's direct solar radiation component. Calculate the angle between the patch normal and the vertical direction, construct the sky visibility factor, and calculate the product of the sky visibility factor and the DHI to obtain the patch's diffuse radiation component. Construct a ground visibility factor that is complementary to the sky visibility factor, calculate the product of the DHI, ground reflectivity, and the ground visibility factor to obtain the patch's ground reflected radiation component. Calculate the sum of the patch's direct solar radiation component, diffuse radiation component, and ground reflected radiation component to obtain the patch's total irradiance.

[0110] Combined with the occlusion surface mark, calculate the irradiance of the patch: The irradiance calculation method of the patch is as follows:

[0111] I total =I direct +I diffuse +I reflected

[0112] Among them, I total is the total irradiance of the patch, which means the total electromagnetic radiation power received by the patch per unit area; I direct I is the solar direct radiation power received by the patch per unit area; diffuse I is the atmospheric scattered radiation power received by the patch per unit area; reflected is the ground reflected radiation power received by the patch per unit area.

[0113] For the occluded surface marked as "1", that is, the surface is not occluded, calculate I direct I direct The calculation method is as follows:

[0114]

[0115] Among them, DNI (Direct Normal Irradiance) is the irradiance of sunlight directly from the solar disk to the surface orthogonal to the light path; is the patch normal vector; is the unit vector of sunlight direction. It is the cosine of the surface normal vector and the direction of sunlight. A negative value means facing away from the sun.

[0116] Calculate I for all faces diffuse I diffuse The calculation method is as follows:

[0117] I diffuse =DHI·F sky

[0118] F sky =(1+cosθ) / 2

[0119] Among them, DHI (Diffuse Horizontal Irradiance) is the atmospheric scattered radiation irradiance received on the horizontal plane; F sky is the sky visibility factor, which represents the visibility between the patch and the sky; θ is the difference between the patch normal vector and the vertical unit vector Angle.

[0120] Calculate I for all faces reflected I reflected The calculation method is as follows:

[0121] I reflected =ρ·GHI·F ground

[0122] F ground =(1-cosθ) / 2

[0123] Where GHI (Global Horizontal Irradiance) is the total horizontal irradiance of the sun’s DHI and DNI reaching the horizontal surface; ρ is the ground reflectivity, which is taken as 0.2 here; F ground is the ground visibility factor, which represents the visibility between the patch and the ground; θ is the difference between the patch normal vector and the vertical unit vector Angle.

[0124] Calculate the I of the patch direct , I diffuse , I reflected The sum of the total irradiance I of the patch at a certain moment is obtained. total .

[0125] (6) Calculate the sum of the shielding marks and the average irradiance at all times of the day to obtain the sunshine hours and average irradiance of the day; take the average of the analysis results of different dates to obtain the average sunshine hours and average irradiance.

[0126] Calculate the sum of the occlusion marks at all times of the day to obtain the sunshine hours for all patches on that day. Calculate the mean irradiance at all times of the day to obtain the average irradiance for that day. Average the analysis results for different days to obtain the total average sunshine hours and average irradiance for the analysis window.

[0127] (7) Output the sunshine analysis result table, map the sunshine analysis results to the color spectrum, embed the model as the color information of the surface, and export the model containing the sunshine analysis results.

[0128] The analysis result table is constructed, including: "face_index": the index information of the triangle face, corresponding to the model triangle face after adaptive subdivision; "illumination": the average sunshine hours of the triangle face, in units of (h), with a minimum value of 0 and a maximum value of 13; "irradiance": the average irradiance of the triangle face, in units of (W / m 2 ). The output is a CSV file. The first 10 lines of the file are as follows:

[0129] face_index illumination irradiance 0 3 43.72594223030238 1 1 43.72594223030199 2 3 43.72594223030242 3 1 43.72594223030254 4 1 43.7259422 5 0 43.72594223030278 6 0 43.72594223030277 7 0 43.72594223030209 8 3 43.72594223030244 9 3 43.72594223030282

[0130] Map the average sunshine hours or average irradiance to the color spectrum, embed it into the model as the surface color information, and export the OBJ model containing the sunshine analysis results.

[0131] The following is an embodiment of the present invention using Trimesh, Pyembree, and Scipy: (1) Triangulation

[0132]

[0133]

[0134]

[0135]

[0136] (2) Single time sunshine analysis

[0137]

[0138]

[0139]

Claims

1. A method for analyzing global sunlight exposure of a building based on a real-scene three-dimensional model, characterized in that: The following steps are involved: (1) Set the time window, analysis days and sampling time list, sample sunshine information, determine the sunshine analysis time series set, and determine the sunlight unit direction vector based on the longitude and latitude of the building target location; (2) Load the target building mesh model, identify and separate the roof and facade panels; (3) Using different resolutions to triangulate the roof and facade slices, respectively. The resolution of the roof slice is greater than that of the facade slice. After triangulation, the roof and facade slices are combined into a model. (4) Based on the sunlight unit direction vector obtained in step (1) and the surface center point obtained by triangulating the surface in step (3), construct an emission ray starting from the surface center point and along the sunlight unit direction vector; perform batch emission ray-grid intersection detection, using "0" to mark the obscured surface and "1" to mark the unobscured surface; (5) Considering the occlusion of the patch, calculate the total irradiance of the patch; (6) Calculate the sum of the shielding marks and the average irradiance at all times of the day to obtain the sunshine hours and average irradiance of the day; take the average of the analysis results of different days to obtain the average sunshine hours and average irradiance; (7) Output the sunshine analysis result table, map the sunshine analysis results to the color spectrum, embed the model as the color information of the surface, and export the model containing the sunshine analysis results.

2. The method for analyzing global sunlight exposure of a building based on a real-scene three-dimensional model according to claim 1, characterized in that: In step (1), the unit direction vector of sunlight is determined by combining the longitude and latitude of the building target's geographical location, as follows: Determine the sunlight azimuth and altitude angle based on the sampling time and the longitude and latitude of the building's target location, and calculate the sunlight unit direction vector S: Among them, α is the solar altitude angle, which ranges from θ is the solar azimuth angle, ranging from [0, 2π); S x 、S y and S z They represent the components of the sunlight unit direction vector in the X, Y, and Z directions respectively.

3. The method for analyzing global sunlight exposure of a building based on a real-scene three-dimensional model according to claim 2, characterized in that: Step (2) is as follows: Calculate the angle θ between the normal vector of the triangle and the vertical direction: Wherein, the range of θ is [0°, 180°]; is the normal vector of the triangle; is the vertical unit vector; Set a roof inclination angle threshold and mark the triangles whose calculation results are less than the roof inclination angle threshold as candidate roof faces; the others are facade faces; Taking the center point of each candidate roof triangle as the starting point and the vertical upward direction as the ray direction, a ray list is constructed; ray and grid intersection detection is performed, the number of hits for each ray is counted, and the triangles corresponding to rays with a hit number greater than 0 are deleted from the candidate roof patches to obtain the final roof patch.

4. The method for analyzing global sunlight exposure of a building based on a real-scene three-dimensional model according to claim 3, characterized in that: In step (2), before identifying and separating the roof facets and the facade facets, the triangle normal vector is corrected, and the triangle normal vector with the tilt angle θ∈[150°, 180°] is reversed and adjusted to face the front side upward.

5. The method for analyzing global sunlight exposure of a building based on a real-scene three-dimensional model according to claim 4, characterized in that: In step (3), the triangular surface is subdivided as follows: Assume that the three vertices of a triangle in three-dimensional space are A, B, and C. Calculate the longest side of the triangle, mark the two vertices where the longest side is located as P1 and P2, and the third vertex as P3; calculate the X-axis direction vector of the target two-dimensional coordinate system Triangle unit normal vector Calculate the Y-axis direction vector of the target two-dimensional coordinate system Determine the coordinate transformation relationship from the three-dimensional space projection to the target two-dimensional plane as follows: Among them, P is the coordinate of a point in the original three-dimensional coordinate system, P project is the coordinate of the point after projection onto the two-dimensional plane; Project the three vertices of the triangle (A, B, and C) to obtain points A2, B2, and C2. Calculate the X and Y coordinate ranges of the projected triangle. Using the resolution of the triangulated surface as the spacing, generate X and Y coordinate sequences, combine them into grid points, and use the barycentric coordinate formula to determine whether the grid point is inside the triangle. Add A2, B2, C2, and the mesh points inside the triangle to the vertex list. Calculate the intersection of the three sides of the projected triangle with parallel lines parallel to the Y axis and spaced at the resolution of the triangulated face, and add the intersections to the vertex list. Perform triangulation on the vertices in the vertex list to generate the subdivided faces.

6. The method for analyzing global sunlight exposure of a building based on a real-scene three-dimensional model according to claim 5, characterized in that: In step (3), the barycentric coordinate formula is used to determine whether the grid point is inside the triangle. The method is as follows: u=1-vw Among them, u, v, and w are the centroid coordinates of point P on the two-dimensional plane for triangles A2, B2, and C2 respectively. P is the coordinate of a point on the two-dimensional plane. If u, v, w ≥ 0 and u + v + w = ​​1, then P is inside or on the edge of triangles A2, B2, and C2.

7. The method for analyzing global sunlight exposure of a building based on a real-scene three-dimensional model according to claim 6, characterized in that: In step (3), the points in the vertex list are converted back to the original three-dimensional coordinate system as follows: P3=u·A+v·B+w·C Wherein, P3 is the coordinate converted back to the three-dimensional coordinate system.

8. The method for analyzing global sunlight exposure of a building based on a real-scene three-dimensional model according to claim 7, characterized in that: In step (4), the surface occlusion judgment method is as follows: based on the sunlight unit direction vector in step (1) and the surface center point obtained by triangulating the surface in step (3), a ray list is constructed with each triangulated surface center point as the starting point and the sunlight direction as the ray direction; batch ray and grid intersection detection is performed, the number of hits for each ray is counted, and the value corresponding to the ray with a hit number greater than 0 is set to "0", indicating that the surface is occluded.

9. The method for analyzing global sunlight exposure of a building based on a real-scene three-dimensional model according to claim 8, characterized in that: In step (5), the total irradiance of the patch I total for I total =I direct +I diffuse +I reflected Among them, I total is the total irradiance of the patch, which means the total electromagnetic radiation power received by the patch per unit area; I direct I is the solar direct radiation power received by the patch per unit area; diffuse I is the atmospheric scattered radiation power received by the patch per unit area; reflected is the ground reflected radiation power received by the patch per unit area; According to the unobstructed surface, calculate I direct : Where DNI is the irradiance of sunlight directly from the solar disk to the surface orthogonal to the light path; is the patch normal vector of the unoccluded face; is the unit vector of sunlight direction. It is the cosine of the surface normal vector and the direction of sunlight. A negative value means facing away from the sun. Based on all faces, calculate I diffuse : I diffuse =DHI·F sky F sky =(1+cosθ) / 2 Where DHI is the atmospheric scattered radiation irradiance received on the horizontal plane; F sky is the sky visibility factor, which represents the visibility between the patch and the sky; θ is the difference between the patch normal vector and the vertical unit vector Angle; Calculate I for all faces reflected : I reflected =ρ·GHI·F ground F ground =(1-cosθ) / 2 Where GHI is the total horizontal irradiance of the sun’s DHI and DNI reaching the horizontal surface; ρ is the ground reflectivity; F ground is the ground visibility factor, which represents the visibility between the patch and the ground; θ is the difference between the patch normal vector and the vertical unit vector Angle.

10. The method for analyzing global sunlight exposure of a building based on a real-scene three-dimensional model according to claim 9, characterized in that: In step (6), the output sunshine analysis result table includes the index information of the triangle face "face_index"; the average sunshine hours of the triangle face "illumination", the unit is h, the minimum value is 0, and the maximum value is 13; the average irradiance of the triangle face "irradiance", the unit is W / m 2 .

Images