Neural radiation field DSM generation method based on two-dimensional triangulation network constraint

Abstract

The invention provides a two-dimensional triangulation network constrained neural radiation field DSM generation method. The method comprises the following steps: step 1, implicit volume density and radiation field modeling; 2, extracting a density field gradient and a surface normal; step 3, constructing a two-dimensional triangulation network; step 4, realizing normal consistency constraint based on the two-dimensional triangulation network; step 5, calculating uncertainty perception rendering loss; and step 6, carrying out end-to-end joint optimization. According to the method, geometric constraint and rendering loss are optimized in an end-to-end combined mode, the root mean square error of three-dimensional reconstruction is reduced compared with a traditional NeRF method, and a high-precision DSM generation scheme is provided for remote sensing application such as urban planning and disaster monitoring.

Description

Technical Field

[0001] The present invention belongs to the technical field of 3D reconstruction, and particularly relates to a method for generating a DSM of a neural radiance field constrained by a 2D triangular mesh. Background Art

[0002] With the rapid development of remote sensing technology and deep learning algorithms, 3D reconstruction based on satellite images has become one of the important means for obtaining geographical information. As an effective expression form of 3D terrain data, the digital surface model (DSM) has wide application value in fields such as urban planning, disaster monitoring, and environmental assessment. Traditional DSM generation methods mainly rely on image matching and multi-view geometry techniques, including methods such as image stitching based on stereo image pairs, dense matching, and shape from structure. However, these methods often face challenges such as height discontinuity, occlusion shadows, and lack of texture in cases where high-rise buildings are dense in cities or large-scale terrain changes occur, resulting in the DSM accuracy being difficult to meet the actual application requirements.

[0003] Traditional image matching algorithms usually perform matching based on pixel grayscale and feature point descriptors (such as SIFT, SURF). They perform well in areas with rich ground object textures, but in flat areas such as building roofs, due to the lack of obvious edges and corners, the matching accuracy drops sharply, resulting in depth recovery errors. At the same time, the height difference between high-rise buildings is significant, and traditional matching methods are insufficient in capturing the discontinuous height information, and are prone to introducing geometric artifacts, causing the DSM surface to be discontinuous or with irregular stepped edges. In addition, large-area shadow regions will be generated by buildings themselves and their surrounding environments, and passive optical sensors are difficult to obtain sufficient reflection information. The image matching algorithm either mismatches or directly skips the shadow areas, resulting in the lack of terrain data or error accumulation in these areas.

[0004] In recent years, the emergence of neural rendering and neural radiance field technologies has provided new ideas for 3D reconstruction of complex scenes. NeRF regards 2D images from different perspectives as projections of the scene radiance field, and obtains continuous radiance and density distributions by optimizing the implicit voxel function, thereby generating high-quality view synthesis and depth maps. The SatNeRF method based on satellite images makes the reconstruction process better take into account terrain undulations and view diversity by introducing a terrain height prior in the NeRF framework. However, SatNeRF itself still has some difficulties in matching building height discontinuity, lack of texture in shadow areas, and flat areas, and there is still much room for improvement in enhancing the detail accuracy.

[0005] To further improve the accuracy of DSM reconstruction from urban satellite images, more robust geometric constraints or prior information need to be introduced during the optimization of NeRF to address the above problems. As an efficient mesh generation technique, the 2D triangular mesh can utilize the known ground or building contour line data to generate high-quality triangular meshes that satisfy the input constraint boundaries. Introducing the vertex and face information of the triangular mesh into the NeRF depth optimization process can provide geometric constraints for the scene point cloud, guiding the model to adopt reasonable depth estimation in both highly variable and flat regions, while enhancing the geometric reconstruction effect in the shadow areas using the mesh structure. Summary of the Invention

[0006] Object of the Invention: The technical problem to be solved by the present invention is to provide a method for generating DSM of neural radiance fields constrained by a 2D triangular mesh, in view of the deficiencies of the prior art, which includes the following steps: Step 1, implicit volumetric density and radiance field modeling; Step 2, extracting the density field gradient and surface normal; Step 3, constructing a 2D triangular mesh; Step 4, implementing normal consistency constraints based on the 2D triangular mesh; Step 5, calculating the uncertainty-aware rendering loss ; Step 6, performing end-to-end joint optimization.

[0007] Step 1 includes: using a deep neural network to implicitly model the 3D scene in the satellite image. The input of the deep neural network includes: the spatial sampling point coordinates x, the camera ray direction d0, the sun illumination direction s, and the temporal embedding t0; through the forward propagation of the deep neural network, the following output parameters are obtained: volumetric density , albedo color , shadow scalar , uncertainty coefficient ; It is represented by the following formula: (1), where represents the deep neural network, represents the forward propagation, ; ; ; represents the real three-dimensional space, that is, x is the sampling point coordinates in the three-dimensional space; represents the unit sphere; represents the volumetric density of the spatial sampling point, predicted by the network, reflecting the material thickness or existence probability at that point; Denotes the albedo color of the spatial sampling point in the observation direction d0, output by the network; ∈[0,1] represents the weight indicating whether the spatial sampling point is in the shadow, with the value range in [0,1], where 0 represents full direct illumination and 1 represents full shadow; For each ray r(t)=o+tD, the volume rendering formula is used to calculate the final color: where o represents the starting point of the ray, i.e., the position of the camera optical center; D represents the unit direction vector, which is the sampling direction starting from the camera; t is the distance parameter along the direction D, usually integrated between the near clipping plane t n and the far clipping plane t f ; Calculate the color value of each ray r(t) : (2), where T(t) represents the transfer function; is the color of the sky ambient light; dt represents the integration variable; t n is the near clipping plane; t f is the far clipping plane; the calculation formula of T(t) is: (3), where exp is the natural exponential function; represents the volume density distribution on the path parameter u∈[t n ,t]; du represents the integration variable.

[0008] Step 2 includes: calculating the volume density with respect to the gradient of the spatial sampling point coordinate x, and obtaining the normal vector of each surface point through backpropagation. The calculation formula of the gradient g(x) is: (4), where represents the gradient calculation for the point of the spatial sampling point coordinate x; By performing automatic differentiation on the volume density of the spatial sampling point , the normal vector is obtained: (5), where >0 is a constant to prevent the denominator from being zero; the negative sign ensures that the normal points outside the object; Through Step 2, a numerical gradient operation is performed on the volume density output by the deep neural network, and a continuous set of surface point positions in the scene is extracted.

[0009] Step 3 includes: Step 3-1, extracting corner points; The Harris operator is used to perform response detection on the satellite image data of the open-source data DFC2019 (Data Fusion Contest 2019). The threshold is set to capture more potential corners, and candidate corner points are obtained. In the Harris operator corner detection, the response threshold T Harris is set to 1% - 5% of the maximum response value. For example, let T Harris = 0.02×H max , where H max is the maximum value of the Harris response in the entire image; or a constant threshold of 0.08 - 0.20 is directly set after normalizing the response. For example, take 0.15; For all detected candidate corner points, they are grouped and clustered according to their spatial positions in the image coordinate system. In each cluster, several corner points with the highest response values are selected, and finally a set of feature corner points with uniform distribution and high response intensity is obtained as candidate nodes for constructing the triangular network; Step 3 - 2, straight-line segment detection and cleaning; The LSD (Line Segment Detector) algorithm is used to extract the most prominent line segments in the image, and the following cleaning is performed on the initial line segment set: Merge: If the included angle between two line segments is less than 10° and the end-point distance is less than 2 pixels, they are merged into one continuous line segment; Extension: If the end-point distance between two line segments is less than 12 pixels, it is judged whether the two line segments approach along the extended direction. If so, intersection points are supplemented along the extended direction of the two line segments; Alignment: When the Euclidean distance between the end points where two line segments should intersect in the pixel coordinate system does not exceed 2 - 3 pixels, it is determined that the end points have offsets, and the offset end points are adsorbed and adjusted so that all line segments intersect correctly at the true edge; Finally, a set of corner points and a set of cleaned line segments are obtained; Step 3 - 3, initialize an empty triangular network; Using all the cleaned line segments as hard constraint edges, the Constrained Delaunay Triangulation (CDT) algorithm with constraint conditions is called on the image plane. First, a non-overlapping triangulation is generated to ensure that all line segments exactly serve as the boundaries of the triangular network; Step 3 - 4, corner point insertion and quality control; Based on the triangular network obtained by initialization, the corner points extracted in Step 3 - 1 are sequentially tried to be inserted into the grid. For each corner point p, first find the triangle where the corner point p falls, and then check the following two conditions: (6), (7), where n(T) represents the number of pixels within triangle T, and is the threshold; is the i-th interior angle of triangle T; when i = 1, 2, 3, it corresponds to the interior angle values of the three vertices respectively; If both conditions are satisfied, the corner point p is added as a constraint point to the grid; Step 3 - 5, constraint update and iteration; Adjust the triangles that conflict with the existing line segments due to the inserted corner points, delete or subdivide the new triangles that do not satisfy formulas (6) and (7) until all corner points are successfully inserted or removed; finally, a 2D triangular mesh is obtained.

[0010] Step 4 includes: Through Step 3, the 2D triangular mesh grid data corresponding to each satellite image of DFC2019 is obtained. The 2D triangular mesh grid data includes all vertex coordinates and triangle indices; for each triangle, calculate the 2D plane center and construct a spatial index using the KD - Tree algorithm KD - Tree; For the set of continuous surface point positions in the scene extracted in Step 2, process them in batches. In each batch, after projecting the 3D coordinates to 2D, query the index of the nearest triangle center near each point through the KD - Tree, specifically including: After determining that a surface point is located inside triangle T, generate neighborhood points in the form of barycentric coordinates by first independently and uniformly sampling in the interval [0, 1] , if , then let , to ensure that is still within the barycentric coordinate domain; then let , at this time satisfies ; where represents the three sets of coefficients that satisfy the barycentric coordinate constraint; the neighborhood sampling point is , where are the 2D coordinates of the corresponding vertices of triangle T; For each surface point, check the queried triangles in turn: Use the vertex coordinates to perform a geometric judgment of the current surface point inside the triangle. Once it is confirmed that the point is located inside a triangle, immediately randomly generate a neighborhood sample point inside the triangle in the form of barycentric coordinates. The height of the sample point is aligned with the original surface point, and it is ensured that the sampling displacement does not exceed the preset threshold.

[0011] If no mesh containing the current surface point P1 can be found in the queried triangle list, neighborhood points are generated for point P1 by means of a small random perturbation to ensure that each surface point can obtain at least one geometric constraint sample (neighborhood points generated within the triangular mesh or randomly perturbed neighborhood points for providing geometric consistency constraints); The density gradient Y1 of the surface point P1 and the density gradient Z1 of the sampling points within the corresponding two-dimensional triangular mesh of the surface point P1 are calculated respectively, and the normal vectors Y2 and Z2 are obtained by normalizing Y1 and Z1 respectively. The average difference between Y2 and Z2 is used as the consistency metric, and the average difference is incorporated into the overall loss function to achieve a smooth constraint on the normal vector extraction result in the network; after the normal vector is extracted, the normal vector difference between the target point and the neighborhood points is calculated, and the consistency of the normal vector is constrained through the consistency loss function.

[0012] In step 4, the form of the consistency loss function is: (8), where s is the set of positions of the sampled surface points; represents the normal vector at the spatial sampling point coordinate x; represents the coordinates of the sampling points within the two-dimensional triangular mesh corresponding to the surface points; represents the coordinate and the normal vector at that place.

[0013] In step 5, the uncertainty-aware rendering loss is calculated using the following formula : (9), where E r represents the expectation over all rays r; represents the volume rendering weight at the k-th sampling position along the ray , which is calculated from the previous cumulative transmittance and the current density; C k represents the predicted color of the k-th sampling point obtained by rendering; represents the ground truth color corresponding to the k-th sampling point; β k represents the uncertainty coefficient at the k-th sampling point, which is output by the network and used to measure the uncertainty degree of the color prediction at that point.

[0014] Step 6 includes: combining the loss functions to form an end-to-end optimization objective, and the final total loss function L is: (10), where is the balance coefficient.

[0015] The present invention also provides an electronic device, including a processor and a memory. The memory stores program codes. When the program codes are executed by the processor, the processor is caused to execute the steps of the method described above.

[0016] The present invention also provides a storage medium storing a computer program or instruction. When the computer program or instruction runs on a computer, the steps of the method described above are executed.

[0017] Beneficial effects: While retaining the advantages of high-quality view synthesis of NeRF, the present invention effectively solves three core problems in urban satellite images with the help of two-dimensional triangular mesh geometric priors: Reconstruction of regions with sudden height changes: By constraining the normal consistency with a triangular mesh, the stepped artifacts at the building edges are significantly suppressed, improving the restoration accuracy of the height jumps at the vertical boundaries of high-rise buildings. Filling of shadow missing regions: Using the triangular mesh topology to guide the geometric reasoning in the shadow regions and combining the uncertainty-aware loss to dynamically balance the illumination differences, filling the terrain data missing caused by shadows in traditional methods. Optimization of flat texture regions: In regions lacking feature points such as rooftops and squares, by adaptively densifying the grid of the triangular mesh and combining the density field gradient constraints, reducing the reconstruction error of flat surfaces. Overall accuracy improvement: By jointly optimizing the geometric constraints and rendering loss end-to-end, the root mean square error of the 3D reconstruction is reduced compared with the traditional NeRF method, providing a high-precision DSM generation solution for remote sensing applications such as urban planning and disaster monitoring. Description of the Drawings

[0018] Figure 1 It is the overall flowchart of the method of the present invention.

[0019] Figure 2 It is the effect diagram of the reconstruction model and the detailed comparison diagram. Detailed Embodiments

[0020] The following further detailed description of the present invention is made in conjunction with the drawings and specific embodiments. The above and / or other advantages of the present invention will become clearer.

[0021] As Figure 1 shown, the embodiment of the present invention provides a method for generating a DSM of a neural radiance field constrained by a two-dimensional triangular mesh, including the following steps: Step 1, implicit volume density and radiance field modeling; Step 2, obtaining the normal vector from the density field gradient; Step 3, constructing a two-dimensional triangular mesh; Step 4, normal vector consistency constraint based on the two-dimensional triangular mesh; Step 5, calculate the uncertainty-aware rendering loss; Step 6, perform end-to-end joint optimization.

[0022] Step 1 includes: using a deep neural network to implicitly model the 3D scene in the satellite image, where the input of the deep neural network includes: the spatial sampling point coordinates x, the camera ray direction d0, the sun illumination direction s, and the temporal embedding t0; through the forward propagation of the deep neural network, the following output parameters are obtained: the volume density , the albedo color , the shadow scalar , the uncertainty coefficient ; It is represented by the following formula: (1), where represents the deep neural network, represents the forward propagation, ; ; ; represents the real three-dimensional space, that is, x is the sampling point coordinates in the three-dimensional space; represents the unit sphere; represents the volume density of the spatial sampling point, predicted by the network, reflecting the material thickness or existence probability at that point; represents the albedo color of the spatial sampling point in the observation direction d0, output by the network; ∈[0,1] represents the weight of whether the spatial sampling point is in the shadow, with a value range in [0,1], where 0 represents full direct sunlight and 1 represents full shadow; For each ray r(t)=o+tD, the volume rendering formula is used to calculate the final color: where o represents the starting point of the ray, that is, the camera optical center position; D represents the unit direction vector, which is the sampling direction starting from the camera; t is the distance parameter along the direction D, usually integrated between the near clipping plane t n and the far clipping plane t f ; Calculate the color value of each ray r(t) : (2), where T(t) represents the transfer function; is the color of the sky ambient light; dt represents the integration variable; t n is the near clipping plane; t f is the far clipping plane; The calculation formula of T(t) is: (3), where exp is the natural exponential function; represents the volume density distribution on the path parameter u ∈ [t n , t]; du represents the integration variable.

[0023] Step 2 includes: calculating the gradient of the volume density with respect to the spatial sampling point coordinates x, and obtaining the normal vector of each surface point through backpropagation. The gradient g(x) calculation formula is: (4), where represents the gradient calculation for the points of the spatial sampling point coordinates x; By performing automatic differentiation on the volume density of the spatial sampling points the normal vector is obtained: (5), where >0 is a constant to prevent the denominator from being zero; the negative sign ensures that the normal points to the outside of the object; Through Step 2, a numerical gradient operation is performed on the volume density output by the deep neural network, and a set of continuous surface point positions in the scene is extracted.

[0024] Step 3 includes: Step 3-1, extracting corner points; The Harris operator is used to perform response detection on the open-source data DFC2019 (Data Fusion Contest 2019) satellite image data, and the threshold is set to capture more potential corners to obtain candidate corner points. In the Harris operator corner detection, the response threshold T Harris is set to 1% - 5% of the maximum response value. For example, let T Harris = 0.02×H max , where H max is the maximum value of the Harris response in the entire image; or a constant threshold of 0.08 - 0.20 is directly set after normalizing the response. For example, take 0.15; For all detected candidate corner points, they are grouped and clustered according to their spatial positions in the image coordinate system. In each cluster, several corner points with the highest response values are selected, and finally a set of uniformly distributed and high-response-strength feature corner points is obtained as candidate nodes for constructing the triangular mesh; Step 3-2, straight-line segment detection and cleaning; Use the Line Segment Detector (LSD) algorithm to extract the most prominent line segments in the image, and clean the initial set of line segments as follows: Merge: If the angle between two line segments is less than 10° and the endpoint distance is less than 2 pixels, merge them into one continuous line segment; Extend: If the endpoint distance between two line segments is less than 12 pixels, determine whether the two line segments approach along the extension direction. If so, supplement the intersection point along the extension direction of the two line segments; Align: When the Euclidean distance between the endpoints where two line segments should intersect does not exceed 2 - 3 pixels in the pixel coordinate system, it is determined that the endpoints have offsets. Adjust the offset endpoints by adsorption so that all line segments intersect correctly at the true edge; Finally, obtain a set of corner points and a set of cleaned line segments; Step 3 - 3, Initialize an empty triangular mesh; Using all the cleaned line segments as hard constraint edges, call the Constrained Delaunay Triangulation (CDT) algorithm on the image plane. First, generate a non - overlapping triangulation that ensures all line segments exactly serve as the boundaries of the triangular mesh; Step 3 - 4, Corner point insertion and quality control; Based on the initialized triangular mesh, sequentially try to insert the corner points extracted in Step 3 - 1 into the mesh. For each corner point p, first find the triangle where the corner point p falls, and then check the following two conditions: (6), (7), where \(n(T)\) represents the number of pixels in triangle \(T\), and are thresholds; is the \(i\) - th interior angle of triangle \(T\); when \(i = 1,2,3\), it corresponds to the interior angle values of the three vertices respectively; is used to control the area (number of pixels) of the triangle, usually set in the range of 200 - 1000 pixels. For example, when the image resolution is 0.5 m / px - 1 m / px, can be taken as about 500 pixels; m / px is meters per pixel.

[0025] is used to exclude overly sharp or degenerate triangles, generally taken between \(\cos(2^{\circ})\approx0.9994\) and \(\cos(3^{\circ})\approx0.9986\). For example, it can be set that = 0.999.

[0026] The corner point p is added to the grid as a constraint point only if both conditions are satisfied. Step 3-5, constraint update and iteration; Adjust the triangles that conflict with the existing line segments due to the insertion of corner points, and delete or subdivide the new triangles that do not satisfy formulas (6) and (7) until all corner points are successfully inserted or removed; finally, a 2D triangular mesh is obtained.

[0027] Step 4 includes: Through Step 3, the 2D triangular mesh grid data corresponding to each satellite image of DFC2019 is obtained. The 2D triangular mesh grid data includes all vertex coordinates and triangle indices; for each triangle, calculate the 2D plane center and construct a spatial index using the KD-Tree algorithm. For the set of continuous surface point positions in the scene extracted in Step 2, process them in batches. In each batch, after projecting the 3D coordinates to 2D, query the index of the nearest triangle center near each point through the KD-Tree. Specifically, after determining that a surface point P1 (with the same meaning as x appearing earlier) is inside a triangle T, generate neighborhood points in the form of barycentric coordinates by independently and uniformly sampling in the interval [0,1] . If , then let , to ensure that is still within the barycentric coordinate domain; then let . At this time, satisfies ; where represents the three sets of coefficients that satisfy the barycentric coordinate constraints; the neighborhood sampling point is , where are the 2D coordinates of the corresponding vertices of triangle T. For each surface point, check the queried triangles in turn: use the vertex coordinates to perform a geometric judgment of the current surface point P1 inside the triangle in the 2D plane. Once it is confirmed that point P1 is inside a triangle, immediately generate a neighborhood sample point randomly in the form of barycentric coordinates inside the triangle. The height of the sample point is aligned with the original surface point, and it is ensured that the sampling displacement does not exceed the preset threshold.

[0028] If no grid containing the current surface point P1 can be found in the queried triangle list, generate neighborhood points for point P1 in a way of small random perturbations to ensure that each surface point can obtain at least one geometric constraint sample (the neighborhood points generated inside the triangular mesh or the neighborhood points generated by random perturbations are used to provide geometric consistency constraints). Calculate the density gradient Y1 of the surface point P1 and the density gradient Z1 of the sampling points within the two-dimensional triangular mesh corresponding to the surface point P1 respectively, and normalize Y1 and Z1 respectively to obtain the normal vectors Y2 and Z2. Use the average difference between Y2 and Z2 as the consistency metric, and incorporate the average difference into the overall loss function to achieve smooth constraints on the normal extraction results in the network; after normal vector extraction, calculate the normal vector difference between the target point and the neighborhood points, and use the consistency loss function to constrain the consistency of the normal vectors.

[0029] In step 4, the form of the consistency loss function is: (8), where s is the set of positions of the sampled surface points; represents the normal vector at the spatial sampling point coordinate x; represents the coordinates of the sampling points within the two-dimensional triangular mesh corresponding to the surface points; represents the coordinate and the normal vector at that position.

[0030] In step 5, the uncertainty-aware rendering loss is calculated using the following formula : (9), where E r represents the expectation over all rays r; represents the volume rendering weight at the k-th sampling position along the ray , which is calculated from the previous cumulative transmittance and the current density; C k represents the predicted color of the k-th sampling point obtained by rendering; represents the ground truth color corresponding to the k-th sampling point; β k represents the uncertainty coefficient at the k-th sampling point, which is output by the network and used to measure the uncertainty degree of the color prediction at that point.

[0031] Step 6 includes: combining the various loss functions to form an end-to-end optimization objective, and the final total loss function L is: (10), where is the balance coefficient.

[0032] Such as Figure 2As shown, it is a comparison chart of the reconstruction effects of the original model and the improved model. In the first row, the "original image" is the satellite image input. In the second row, the "original model" shows that there are stepped errors at the edges of buildings with abrupt height changes and blurred textures in the shadow areas. In the third row, the "improved model" significantly improves the flatness of the vertical surfaces of buildings, the geometric consistency of the shadow areas, and the detail fidelity of the flat roofs through the constraint of the two-dimensional triangular network, verifying the reconstruction advantages of this method for complex scenes.

[0033] The present invention provides a method for generating a digital surface model (DSM) of a neural radiance field constrained by a two-dimensional triangular network. There are many methods and ways to specifically implement this technical solution. The above is only the preferred embodiment of the present invention. It should be noted that for those of ordinary skill in the art of this technology, without departing from the principle of the present invention, several improvements and refinements can be made, and these improvements and refinements should also be regarded as the protection scope of the present invention. Each component not clearly defined in this embodiment can be implemented by existing technologies.

Claims

1. A method for generating a DSM of a neural radiance field constrained by a two-dimensional triangular mesh, characterized in that It includes the following steps: Step 1, implicit volumetric density and radiation field modeling; Step 2, obtaining the normal vector from the density field gradient; Step 3, constructing a 2D triangular mesh; Step 4, normal vector consistency constraint based on the 2D triangular mesh; Step 5, calculating the uncertainty-aware rendering loss; Step 6, performing end-to-end joint optimization.

2. The method according to claim 1, wherein Step 1 includes: using a deep neural network to implicitly model the 3D scene in the satellite image, where the input of the deep neural network includes: spatial sampling point coordinates x, camera ray direction d0, solar illumination direction s, and temporal embedding t0; through the forward propagation of the deep neural network, the following output parameters are obtained: volume density , albedo color , shadow scalar , uncertainty coefficient ; It is represented by the following formula: （1）， Among them represents a deep neural network, represents forward propagation, ; ; ; represents a three-dimensional real space; represents a unit sphere; Represents the volumetric density of spatial sampling points; Represents the albedo color of the spatial sampling point in the observation direction d0; ∈[0,1] represents the weight indicating whether it is in the shadow at the spatial sampling point, and the value range is [0,1], where 0 represents full direct sunlight and 1 represents full shadow; For each ray r(t)=o+tD, the final color is calculated using the volume rendering formula: where o represents the starting point of the ray; D represents the unit direction vector; t is the distance parameter along the direction D; Calculate the color value of each ray r(t) :[[-END]] （2）， Among them, T(t) represents the transfer function; is the color of the sky ambient light; dt represents the integration variable; t n is the near clipping plane; t f is the far clipping plane; The calculation formula of T(t) is: （3）， where exp is the natural exponential function; denotes the volume density distribution over the path parameter u ∈ [t n , t]; du represents the integration variable.

3. The method according to claim 2, wherein Step 2 includes: calculating the bulk density the gradient with respect to the spatial sampling point coordinate x, and obtaining the normal vector of each surface point through backpropagation. The gradient g(x) calculation formula is: （4）， Among them means to perform gradient calculation on the points with the coordinates x of the spatial sampling points; By automatically differentiating the volumetric density of spatial sampling points a normal vector is obtained : （5）， Among them > 0 is a constant to prevent the denominator from being zero; Through Step 2, perform a numerical gradient operation on the volumetric density output by the deep neural network to extract the set of continuous surface point positions in the scene.

4. The method according to claim 3, characterized in that, Step 3 includes: Step 3-1, extracting corner points; Use the Harris operator to perform response detection on satellite image data, set the threshold, and obtain candidate corner points. In the Harris operator corner detection, set the response threshold T Harris to 1% - 5% of the maximum response value H max , or directly set a constant threshold after normalizing the response; For all detected candidate corner points, group and cluster them according to their spatial positions in the image coordinate system. Select the corner point with the highest response value in each cluster. Finally, obtain a set of uniformly distributed and highly responsive feature corner points as candidate nodes for constructing the triangular mesh; Step 3-2, straight line segment detection and cleaning; Use the LSD algorithm of the straight line segment detector to extract the most prominent line segments in the image, and perform the following cleaning on the initial line segment set: Merge: If the included angle between two line segments is less than 10° and the endpoint distance is less than 2 pixels, merge them into one continuous line segment; Extend: If the endpoint distance between two line segments is less than 12 pixels, determine whether the two line segments approach along the extended direction. If so, supplement the intersection point along the extended direction of the two line segments; Align: When the Euclidean distance between the endpoints where two line segments should intersect in the pixel coordinate system does not exceed 2 - 3 pixels, it is determined that the endpoints have offsets. Adjust the offset endpoints by adsorption so that all line segments intersect correctly at the real edge; Finally, obtain a set of corner points and a set of cleaned line segments; Step 3-3, initializing an empty triangular mesh; Using all the cleaned line segments as hard constraint edges, call the Constrained Delaunay Triangulation (CDT) algorithm with constraints on the image plane. First, generate a non-overlapping triangulation that ensures all line segments exactly serve as the boundaries of the triangular mesh; Step 3-4, corner point insertion and quality control; Based on the triangular mesh obtained by initialization, sequentially try to insert the corner points extracted in Step 3-1 into the mesh. For each corner point p, first find the triangle where the corner point p falls, and then check the following two conditions: （6）， （7）， where n(T) represents the number of pixels within triangle T, and is the threshold; is the i-th interior angle of triangle T; when i = 1, 2, 3, it corresponds to the interior angle values of the three vertices respectively; Only when both conditions are satisfied, add the corner point p as a constraint point to the mesh; Step 3-5, constraint update and iteration; Adjust the triangles that conflict with the existing line segments due to the insertion of corner points, delete or subdivide the new triangles that do not satisfy formulas (6) and (7) until all corner points are successfully inserted or removed; finally, obtain the 2D triangular mesh.

5. The method according to claim 4, characterized in that Step 4 includes: Through Step 3, the 2D triangular mesh grid data corresponding to each satellite image is obtained. The 2D triangular mesh grid data includes all vertex coordinates and triangle indices; for each triangle, calculate the 2D plane center and construct a spatial index using the KD-Tree algorithm; For the set of consecutive surface point positions in the scene obtained by step 2, process them in batches. In each batch, after projecting the three-dimensional coordinates into two dimensions, query the indices of the centers of the nearest triangles for each point through a KD-Tree. Specifically, it includes: After determining that a surface point is inside triangle T, generate neighborhood points in the form of barycentric coordinates . First, independently and uniformly sample in the interval [0, 1] . If , then let , to ensure that remains within the barycentric coordinate domain; then let . At this time, satisfies ; where represents the three sets of coefficients that satisfy the barycentric coordinate constraints; the neighborhood sampling points are , where are the two-dimensional coordinates of the corresponding vertices of triangle T. For each surface point, check the queried triangles in sequence: perform the geometric judgment of the current surface point within the triangle using the vertex coordinates in the two-dimensional plane Once it is confirmed that the point is located inside a triangle, immediately generate a neighborhood sample point randomly in the triangle in barycentric coordinates. The height of the sample point is aligned with the original surface point, and it is ensured that the sampling displacement does not exceed the preset threshold; If a grid containing the current surface point P1 cannot be found in the queried triangle list, neighborhood points are generated for point P1 by means of random perturbation to ensure that each surface point can obtain at least one geometric constraint sample; The density gradient Y1 of the surface point P1 and the density gradient Z1 of the in-sampled points within the two-dimensional triangular mesh corresponding to the surface point P1 are calculated respectively, and the normal vectors Y2 and Z2 are obtained by normalizing Y1 and Z1 respectively. The average difference between Y2 and Z2 is used as the consistency metric, and the average difference is incorporated into the overall loss function to achieve smooth constraint on the normal extraction result in the network; after the normal vector is extracted, the normal vector difference between the target point and the neighborhood points is calculated, and the consistency of the normal vector is constrained by the consistency loss function.

6. The method according to claim 5, wherein In step 4, the form of the consistency loss function is: （8）， where s is the set of sampled surface point positions; represents the normal vector at the spatial sampling point coordinate x; represents the coordinates of the sampling points within the two-dimensional triangular mesh corresponding to the surface points; represents the coordinate and the normal vector at that point.

7. The method according to claim 6, wherein In step 5, the uncertainty-aware rendering loss is calculated using the following formula :[[]]END]] （9）， Among them, E r represents the expectation for all rays r; represents along the ray the volume rendering weight at the k-th sampling position; C k represents the predicted color of the k-th sampling point obtained by rendering; represents the ground truth color corresponding to the k-th sampling point; β k represents the uncertainty coefficient at the k-th sampling point.

8. The method according to claim 7, wherein Step 6 includes: combining the loss functions to form an end-to-end optimization objective, and the final total loss function L is: （10）， wherein is the balance coefficient.

9. An electronic device, characterized in that, It includes a processor and a memory, and the memory stores program codes. When the program codes are executed by the processor, the processor is caused to execute the steps of the method according to any one of claims 1 to 8.

10. A storage medium, characterized in that, It stores a computer program or instruction. When the computer program or instruction runs on a computer, it executes the steps of the method according to any one of claims 1 to 8.

Images