A 2D triangulated network constrained neural radiation field DSM generation method

Abstract

The present invention provides a method for generating a neural radiation field (DSM) constrained by a two-dimensional triangulated mesh, comprising the following steps: Step 1: implicit volume density and radiation field modeling; Step 2: extracting density field gradients and surface normals; Step 3: constructing a two-dimensional triangulated mesh; Step 4: implementing normal consistency constraints based on the two-dimensional triangulated mesh; Step 5: calculating uncertainty-aware rendering loss; and Step 6: performing end-to-end joint optimization. By jointly optimizing geometric constraints and rendering loss in an end-to-end manner, the method achieves a lower root mean square error (RMS) of three-dimensional reconstruction than traditional NeRF methods, providing a high-precision DSM generation solution for remote sensing applications such as urban planning and disaster monitoring.

Description

Technical Field

[0001] The present invention belongs to the technical field of three-dimensional reconstruction, and in particular relates to a method for generating a neural radiation field DSM constrained by a two-dimensional triangulated network. Background Art

[0002] With the rapid development of remote sensing technology and deep learning algorithms, three-dimensional reconstruction based on satellite images has become one of the important means of obtaining geographic information. As an effective form of expression of three-dimensional terrain data, digital terrain models have wide application value in urban planning, disaster monitoring, environmental assessment and other fields. Traditional DSM generation methods mainly rely on image matching and multi-view geometry techniques, including image stitching based on stereo pairs, dense matching, and shape from structure. However, these methods often face challenges such as height discontinuity, occlusion and shadow, and lack of texture in urban areas with dense high-rise buildings or large-scale terrain changes, making it difficult for DSM accuracy to meet the needs of actual applications.

[0003] Traditional image matching algorithms typically match based on pixel grayscale and feature point descriptors (such as SIFT and SURF). While they perform reasonably well in areas with rich terrain texture, matching accuracy plummets in flat areas like building rooftops due to the lack of distinct edges and corners, leading to depth recovery errors. Furthermore, due to the significant height differences between high-rise buildings, traditional matching methods fail to capture sudden height changes, easily introducing geometric artifacts that result in discontinuous DSM surfaces or irregular stepped edges. Furthermore, buildings and their surroundings cast large shadows, making it difficult for passive optical sensors to obtain sufficient reflectance information. Image matching algorithms either mismatch or skip over shadowed areas, resulting in missing terrain data or accumulating errors in these areas.

[0004] The recent emergence of neural rendering and neural radiance field technology has provided new insights for the three-dimensional reconstruction of complex scenes. NeRF treats two-dimensional images from different perspectives as projections of the scene's radiance field. By optimizing implicit voxel functions, it obtains continuous radiance and density distributions, thereby generating high-quality view synthesis and depth maps. The satellite imagery-based SatNeRF method introduces a terrain height prior into the NeRF framework, enabling the reconstruction process to better account for terrain undulations and perspective diversity. However, SatNeRF itself still has some difficulties matching discontinuous building heights, missing texture in shadow areas, and flat areas, and there is still considerable room for improvement in detail accuracy.

[0005] To further improve the accuracy of DSM reconstruction from urban satellite imagery, it is necessary to introduce more robust geometric constraints or prior information into the NeRF optimization process to address the aforementioned issues. Two-dimensional triangulation (TMN), as an efficient meshing technique, can leverage known ground or building outline data to generate high-quality triangulated meshes that meet input constraint boundaries. Introducing triangulated mesh vertex and face information into the NeRF depth optimization process provides geometric constraints for the scene point cloud, guiding the model to achieve reasonable depth estimates in both highly abrupt and flat areas, while also leveraging the mesh structure to enhance geometric reconstruction in shadowed areas. Summary of the Invention

[0006] Purpose of the invention: The technical problem to be solved by the present invention is to address the deficiencies of the existing technology and provide a method for generating a neural radiation field DSM constrained by a two-dimensional triangulated network, comprising the following steps:

[0007] Step 1: implicit volume density and radiation field modeling;

[0008] Step 2, extract the density field gradient and surface normal;

[0009] Step 3, construct a two-dimensional triangulated network;

[0010] Step 4: Implement normal consistency constraint based on 2D triangulated network;

[0011] Step 5: Calculate uncertainty-aware rendering loss ;

[0012] Step 6: Perform end-to-end joint optimization.

[0013] Step 1 includes: using a deep neural network to implicitly model the three-dimensional scene in the satellite image. The input of the deep neural network includes: spatial sampling point coordinates x, camera light direction d0, sun radiation direction s and time series 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 ;

[0014] It is expressed by the following formula:

[0015] (1),

[0016] in represents a deep neural network, represents forward propagation, ; ; ; Represents real three-dimensional space, that is, x is the coordinate of the sampling point in three-dimensional space; represents the unit sphere;

[0017] Represents the volume density of a spatial sampling point, predicted by the network, reflecting the material thickness or existence probability at that point;

[0018] Represents the albedo color of the spatial sampling point in the viewing direction d0, which is output by the network;

[0019] ∈[0,1] represents the weight of whether the spatial sampling point is in the shadow, and the value range is [0,1], where 0 represents complete direct lighting and 1 represents complete shadow;

[0020] For each ray r(t)=o+tD, the final color is calculated using the volume rendering formula:

[0021] Where o represents the starting point of the ray, that is, the position of the camera's optical center; D represents the unit direction vector, which is the sampling direction from the camera; t is the distance parameter along the direction D, usually at the near clipping plane t n and the far clipping plane t f Integrate between

[0022] Calculate the color value of each ray r(t) :

[0023] (2),

[0024] Where T(t) represents the transfer function; is the color of the sky ambient light; dt represents the integral variable; t n is the near clipping plane; t f is the far clipping plane; the calculation formula of T(t) is:

[0025] (3),

[0026] Where exp is the natural exponential function; Indicates that the path parameter u∈[t n ,t]; du represents the integral variable.

[0027] Step 2 includes: calculating the volume density The gradient relative to the spatial sampling point coordinate x is used to obtain the normal vector of each surface point through back propagation. The gradient g(x) is calculated as follows:

[0028] (4),

[0029] in Indicates that the gradient calculation is performed on the point with the spatial sampling point coordinate x;

[0030] By the volume density of spatial sampling points Perform automatic differentiation to obtain the normal vector :

[0031] (5),

[0032] in >0 is a constant to prevent the denominator from being zero; the negative sign ensures that the normal points outside the object;

[0033] Through step 2, a numerical gradient operation is performed on the volume density output by the deep neural network to extract a continuous set of surface point positions in the scene.

[0034] Step 3 includes:

[0035] Step 3-1, extract corner points;

[0036] Harris operator is used to detect the response of open source data DFC2019 (Data Fusion Contest 2019) satellite image data. The threshold is set to capture more potential corners and obtain candidate corner points. In the Harris operator corner point detection, the response threshold T is set to 0. Harris Set it to 1% to 5% of the maximum response value, for example, let T Harris =0.02×H max , where H max The maximum value of the Harris response in the entire image; or after normalizing the response, directly set a constant threshold of 0.08 to 0.20, for example, 0.15;

[0037] For all detected candidate corner points, group and cluster them according to their spatial positions in the image coordinate system. In each cluster, select several corner points with the highest response values. Finally, a set of feature corner points with uniform distribution and high response intensity are obtained as candidate nodes for constructing the triangulated network.

[0038] Step 3-2, straight line segment detection and cleaning;

[0039] Use the Line Segment Detector (LSD) algorithm to extract the most significant line segments in the image and clean up the initial line segment set as follows:

[0040] Merge: If the angle between two line segments is less than 10° and the distance between their endpoints is less than 2 pixels, they are merged into a continuous line segment;

[0041] Extension: If the distance between the endpoints of two line segments is less than 12 pixels, determine whether the two line segments are approaching along the extension direction. If so, add the intersection point along the extension direction of the two line segments;

[0042] Alignment: When the Euclidean distance between the endpoints of two line segments that should intersect in the pixel coordinate system does not exceed 2 to 3 pixels, the endpoints are considered offset and the offset endpoints are adjusted so that all line segments intersect correctly at the true edge.

[0043] Finally, we get a set of corner points and a set of cleaned line segments;

[0044] Step 3-3, initialize the empty triangulated network;

[0045] Using all cleaned line segments as hard constraint edges, the Constrained Delaunay Triangulation (CDT) algorithm with constraints is called on the image plane. First, a non-overlapping triangulation is generated to ensure that all line segments are exactly as the triangulation boundaries.

[0046] Step 3-4, corner point insertion and quality control;

[0047] Based on the initialized triangulated mesh, try to insert the corner points extracted in step 3-1 into the mesh in sequence. For each corner point p, first find the triangle where the corner point p falls, and then check the following two conditions:

[0048] (6),

[0049] (7),

[0050] Where n(T) represents the number of pixels in triangle T, and is the threshold; is the i-th interior angle of triangle T; i=1, 2, 3 correspond to the interior angle values ​​of the three vertices respectively;

[0051] If both conditions are met, the corner point p is added to the grid as a constraint point;

[0052] Steps 3-5, constraint update and iteration;

[0053] Adjust the triangles that conflict with 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 eliminated; finally, a two-dimensional triangulated network is obtained.

[0054] Step 4 includes: obtaining the two-dimensional triangulated mesh data corresponding to each satellite image of DFC2019 through step 3, wherein the two-dimensional triangulated mesh data includes all vertex coordinates and triangle indexes; for each triangle, calculating the two-dimensional plane center and constructing a spatial index using the K-dimensional tree algorithm KD-Tree;

[0055] For the continuous surface point position set in the scene extracted in step 2, it is processed in batches. In each batch, the three-dimensional coordinates are projected to two dimensions, and the nearest triangle center index near each point is queried through the KD-Tree. Specifically, the following steps are performed: After being inside triangle T, the barycentric coordinates Generate neighborhood points by sampling independently and uniformly in the interval [0,1] ,if , then let , to ensure Still in the center of gravity coordinate domain; let ,at this time satisfy ;in Represents three sets of coefficients that satisfy the barycentric coordinate constraints; the neighborhood sampling points are ,in are the two-dimensional coordinates of the corresponding vertices of triangle T;

[0056] For each surface point, check the queried triangles in turn: use the vertex coordinates to perform the current surface point on the two-dimensional plane Geometric judgment within the triangle, once the point is confirmed Located inside a triangle, a neighborhood sample point is randomly generated inside the triangle in the form of barycentric coordinates. The height of the sample point is aligned with the original surface point, and the sampling displacement is guaranteed not to exceed the preset threshold.

[0057] If the mesh containing the current surface point P1 cannot be found in the queried triangle list, a small random perturbation method is used to generate neighboring points for point P1 to ensure that each surface point can obtain at least one geometric constraint sample (neighboring points generated in the triangulation network or neighborhood points generated by random perturbation to provide geometric consistency constraints);

[0058] The density gradient Y1 of the surface point P1 and the density gradient Z1 of the sampling point in the two-dimensional triangulation network corresponding to the surface point P1 are calculated respectively, and Y1 and Z1 are normalized to obtain the normal vectors Y2 and Z2 respectively. The average difference between Y2 and Z2 is used as the consistency measure, and the average difference is incorporated into the overall loss function to achieve smooth constraints on the normal extraction results in the network; after the normal vector is extracted, the normal vector difference between the target point and the neighboring points is calculated, and the consistency of the normal vector is constrained by the consistency loss function.

[0059] In step 4, the consistency loss function is in the form of:

[0060] (8),

[0061] Where s is the set of sampled surface point positions; Represents the normal vector at the spatial sampling point coordinate x; Indicates the coordinates of the sampling point in the two-dimensional triangulation corresponding to the surface point; Representing coordinates The normal vector at .

[0062] In step 5, the uncertainty-aware rendering loss is calculated using the following formula: :

[0063] (9),

[0064] Among them E r represents the expectation of all rays r; Indicates along the ray The volume rendering weight at the kth sampling position is calculated by the previous cumulative transmittance and the current density; C k Indicates 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 kth sampling point, which is output by the network and is used to measure the uncertainty of the color prediction at that point.

[0065] Step 6 includes: combining the various loss functions to form an end-to-end optimization objective. The final total loss function L is:

[0066] (10),

[0067] in is the balance coefficient.

[0068] The present invention also provides an electronic device, comprising a processor and a memory, wherein the memory stores program code, and when the program code is executed by the processor, the processor executes the steps of the method.

[0069] The present invention also provides a storage medium storing a computer program or instruction, which executes the steps of the method when the computer program or instruction is run on a computer.

[0070] Beneficial Effects: While retaining the advantages of NeRF high-quality view synthesis, this invention effectively solves three core problems in urban satellite imagery by leveraging 2D triangulated mesh geometry priors:

[0071] Reconstruction of height-change regions: By constraining normal consistency through triangulation, the stair-step artifacts at building edges are significantly suppressed, improving the accuracy of height-change reconstruction at the vertical boundaries of high-rise buildings.

[0072] Filling in missing shadow areas: Using triangulated mesh topology to guide geometric reasoning of shadow areas, combined with uncertainty-aware loss to dynamically balance illumination differences, filling in the gaps in terrain data caused by shadows in traditional methods.

[0073] Optimization of flat textured areas: In areas lacking feature points, such as rooftops and squares, the reconstruction error of flat surfaces is reduced by adaptively refining the triangulated mesh and combining it with density field gradient constraints.

[0074] Overall accuracy improvement: Through end-to-end joint optimization of geometric constraints and rendering loss, the root mean square error of 3D reconstruction is reduced compared to traditional NeRF methods, providing a high-precision DSM generation solution for remote sensing applications such as urban planning and disaster monitoring. BRIEF DESCRIPTION OF THE DRAWINGS

[0075] Figure 1 4 is an overall flow chart of the method of the present invention.

[0076] Figure 2 The following are the renderings of the reconstructed model and the detailed comparison chart. DETAILED DESCRIPTION

[0077] The present invention will be further described below in conjunction with the accompanying drawings and specific embodiments, and the above and / or other advantages of the present invention will become more apparent.

[0078] like Figure 1 As shown, an embodiment of the present invention provides a method for generating a neural radiation field DSM constrained by a two-dimensional triangulated network, comprising the following steps:

[0079] Step 1: implicit volume density and radiation field modeling;

[0080] Step 2, obtain the normal vector from the density field gradient;

[0081] Step 3, construct a two-dimensional triangulated network;

[0082] Step 4: Normal vector consistency constraint based on 2D triangulated network;

[0083] Step 5: Calculate the uncertainty-aware rendering loss;

[0084] Step 6: Perform end-to-end joint optimization.

[0085] Step 1 includes: using a deep neural network to implicitly model the three-dimensional scene in the satellite image. The input of the deep neural network includes: spatial sampling point coordinates x, camera light direction d0, sun radiation direction s and time series 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 ;

[0086] It is expressed by the following formula:

[0087] (1),

[0088] in represents a deep neural network, represents forward propagation, ; ; ; Represents real three-dimensional space, that is, x is the coordinate of the sampling point in three-dimensional space; represents the unit sphere;

[0089] Represents the volume density of a spatial sampling point, predicted by the network, reflecting the material thickness or existence probability at that point;

[0090] Represents the albedo color of the spatial sampling point in the viewing direction d0, which is output by the network;

[0091] ∈[0,1] represents the weight of whether the spatial sampling point is in the shadow, and the value range is [0,1], where 0 represents complete direct lighting and 1 represents complete shadow;

[0092] For each ray r(t)=o+tD, the final color is calculated using the volume rendering formula:

[0093] Where o represents the starting point of the ray, that is, the position of the camera's optical center; D represents the unit direction vector, which is the sampling direction from the camera; t is the distance parameter along the direction D, usually at the near clipping plane t n and the far clipping plane t f Integrate between

[0094] Calculate the color value of each ray r(t) :

[0095] (2),

[0096] Where T(t) represents the transfer function; is the color of the sky ambient light; dt represents the integral variable; t n is the near clipping plane; t f is the far clipping plane; the calculation formula of T(t) is:

[0097] (3),

[0098] Where exp is the natural exponential function; Indicates that the path parameter u∈[t n ,t]; du represents the integral variable.

[0099] Step 2 includes: calculating the volume density The gradient relative to the spatial sampling point coordinate x is used to obtain the normal vector of each surface point through back propagation. The gradient g(x) is calculated as follows:

[0100] (4),

[0101] in Indicates that the gradient calculation is performed on the point with the spatial sampling point coordinate x;

[0102] By the volume density of spatial sampling points Perform automatic differentiation to obtain the normal vector :

[0103] (5),

[0104] in >0 is a constant to prevent the denominator from being zero; the negative sign ensures that the normal points outside the object;

[0105] Through step 2, a numerical gradient operation is performed on the volume density output by the deep neural network to extract a continuous set of surface point positions in the scene.

[0106] Step 3 includes:

[0107] Step 3-1, extract corner points;

[0108] Harris operator is used to detect the response of open source data DFC2019 (Data Fusion Contest 2019) satellite image data. The threshold is set to capture more potential corners and obtain candidate corner points. In the Harris operator corner point detection, the response threshold T is set to 0.Harris Set it to 1% to 5% of the maximum response value, for example, let T Harris =0.02×H max , where H max The maximum value of the Harris response in the entire image; or after normalizing the response, directly set a constant threshold of 0.08 to 0.20, for example, 0.15;

[0109] For all detected candidate corner points, group and cluster them according to their spatial positions in the image coordinate system. In each cluster, select several corner points with the highest response values. Finally, a set of feature corner points with uniform distribution and high response intensity are obtained as candidate nodes for constructing the triangulated network.

[0110] Step 3-2, straight line segment detection and cleaning;

[0111] Use the Line Segment Detector (LSD) algorithm to extract the most significant line segments in the image and clean up the initial line segment set as follows:

[0112] Merge: If the angle between two line segments is less than 10° and the distance between their endpoints is less than 2 pixels, they are merged into a continuous line segment;

[0113] Extension: If the distance between the endpoints of two line segments is less than 12 pixels, determine whether the two line segments are approaching along the extension direction. If so, add the intersection point along the extension direction of the two line segments;

[0114] Alignment: When the Euclidean distance between the endpoints of two line segments that should intersect in the pixel coordinate system does not exceed 2 to 3 pixels, the endpoints are considered offset and the offset endpoints are adjusted so that all line segments intersect correctly at the true edge.

[0115] Finally, we get a set of corner points and a set of cleaned line segments;

[0116] Step 3-3, initialize the empty triangulated network;

[0117] Using all cleaned line segments as hard constraint edges, the Constrained Delaunay Triangulation (CDT) algorithm with constraints is called on the image plane. First, a non-overlapping triangulation is generated to ensure that all line segments are exactly as the triangulation boundaries.

[0118] Step 3-4, corner point insertion and quality control;

[0119] Based on the initialized triangulated mesh, try to insert the corner points extracted in step 3-1 into the mesh in sequence. For each corner point p, first find the triangle where the corner point p falls, and then check the following two conditions:

[0120] (6),

[0121] (7),

[0122] Where n(T) represents the number of pixels in triangle T, and is the threshold; is the i-th interior angle of triangle T; i=1, 2, 3 correspond to the interior angle values ​​of the three vertices respectively;

[0123] Used to control the triangle area (number of pixels), usually set in the range of 200 to 1000 pixels. For example, when the image resolution is 0.5 m / px~1 m / px, Take about 500 pixels; m / px is meter / pixel.

[0124] It is used to exclude triangles that are too sharp or degenerate. Generally, it is set between cos(2°)≈0.9994 and cos(3°)≈0.9986. For example, = 0.999.

[0125] If both conditions are met, the corner point p is added to the grid as a constraint point;

[0126] Steps 3-5, constraint update and iteration;

[0127] Adjust the triangles that conflict with 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 eliminated; finally, a two-dimensional triangulated network is obtained.

[0128] Step 4 includes: obtaining the two-dimensional triangulated mesh data corresponding to each satellite image of DFC2019 through step 3, wherein the two-dimensional triangulated mesh data includes all vertex coordinates and triangle indexes; for each triangle, calculating the two-dimensional plane center and constructing a spatial index using the K-dimensional tree algorithm KD-Tree;

[0129] For the continuous surface point position set in the scene extracted in step 2, it is processed in batches. In each batch, the three-dimensional coordinates are projected to two dimensions, and the nearest triangle center index near each point is queried through the KD-Tree. Specifically, after determining that a surface point P1 (with the same meaning as x in the previous text) is located inside the triangle T, the center coordinates of the triangle are used to find the nearest triangle center index. Generate neighborhood points by sampling independently and uniformly in the interval [0,1] ,if , then let , to ensure Still in the center of gravity coordinate domain; let ,at this time satisfy ;in Represents three sets of coefficients that satisfy the barycentric coordinate constraints; the neighborhood sampling points are ,in are the two-dimensional coordinates of the corresponding vertices of triangle T;

[0130] For each surface point, the queried triangles are checked in turn: the vertex coordinates are used to perform geometric judgment on the two-dimensional plane to determine whether the current surface point P1 is inside the triangle. Once it is confirmed that point P1 is inside a triangle, a neighborhood sample point is randomly generated inside the triangle using barycentric coordinates. The height of the sample point is aligned with the original surface point, and the sampling displacement is guaranteed not to exceed the preset threshold.

[0131] If the mesh containing the current surface point P1 cannot be found in the queried triangle list, a small random perturbation method is used to generate neighboring points for point P1 to ensure that each surface point can obtain at least one geometric constraint sample (neighboring points generated in the triangulation network or neighborhood points generated by random perturbation to provide geometric consistency constraints);

[0132] The density gradient Y1 of the surface point P1 and the density gradient Z1 of the sampling point in the two-dimensional triangulation network corresponding to the surface point P1 are calculated respectively, and Y1 and Z1 are normalized to obtain the normal vectors Y2 and Z2 respectively. The average difference between Y2 and Z2 is used as the consistency measure, and the average difference is incorporated into the overall loss function to achieve smooth constraints on the normal extraction results in the network; after the normal vector is extracted, the normal vector difference between the target point and the neighboring points is calculated, and the consistency of the normal vector is constrained by the consistency loss function.

[0133] In step 4, the consistency loss function is in the form of:

[0134] (8),

[0135] Where s is the set of sampled surface point positions; Represents the normal vector at the spatial sampling point coordinate x; Indicates the coordinates of the sampling point in the two-dimensional triangulation corresponding to the surface point; Representing coordinates The normal vector at .

[0136] In step 5, the uncertainty-aware rendering loss is calculated using the following formula: :

[0137] (9),

[0138] Among them E r represents the expectation of all rays r; Indicates along the ray The volume rendering weight at the kth sampling position is calculated by the previous cumulative transmittance and the current density; C k Indicates 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 kth sampling point, which is output by the network and is used to measure the uncertainty of the color prediction at that point.

[0139] Step 6 includes: combining the various loss functions to form an end-to-end optimization objective. The final total loss function L is:

[0140] (10),

[0141] in is the balance coefficient.

[0142] like Figure 2 The following figure compares the reconstruction results of the original model and the improved model. The first row, "Original Image," shows satellite imagery input. The second row, "Original Model," shows that traditional NeRF suffers from stair-step errors around highly abrupt building edges and blurry textures in shadowed areas. The third row, "Improved Model," employs 2D triangulated mesh constraints to significantly improve the flatness of vertical building surfaces, the geometric consistency of shadowed areas, and the fidelity of flat roof details, validating the advantages of this method for reconstructing complex scenes.

[0143] The present invention provides a method for generating a neural radiation field (DSM) constrained by a two-dimensional triangulated mesh. There are many methods and approaches for implementing this technical solution. The above is only a preferred embodiment of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and such improvements and modifications should also be considered within the scope of protection of the present invention. Any components not specified in this embodiment may be implemented using existing technologies.

Claims

1. A method for generating a neural radiation field DSM constrained by a two-dimensional triangulated network, characterized in that: The following steps are involved: Step 1: implicit volume density and radiation field modeling; Step 2: Obtain the normal vector from the density field gradient: Perform a numerical gradient operation on the volume density output by the deep neural network to extract a continuous set of surface point positions in the scene; Step 3, construct a two-dimensional triangulated network; Step 4: Normal vector consistency constraint based on 2D triangulated network; Step 4 includes: obtaining two-dimensional triangulated mesh data corresponding to each satellite image through step 3, wherein the two-dimensional triangulated mesh data includes all vertex coordinates and triangle indexes; calculating the two-dimensional plane center of each triangle, and constructing a spatial index using a K-dimensional tree algorithm (KD-Tree); For the continuous surface point position set in the scene extracted in step 2, it is processed in batches. In each batch, the three-dimensional coordinates are projected to two dimensions, and the nearest triangle center index near each point is queried through the KD-Tree. Specifically, the following steps are performed: After being inside triangle T, the barycentric coordinates Generate neighborhood points by sampling independently and uniformly in the interval [0,1] ,if , then let , to ensure Still in the center of gravity coordinate domain; let ,at this time satisfy ;in Represents three sets of coefficients that satisfy the barycentric coordinate constraints; the neighborhood sampling points are ,in are the two-dimensional coordinates of the corresponding vertices of triangle T; For each surface point, check the queried triangles in turn: use the vertex coordinates to perform the current surface point on the two-dimensional plane Geometric judgment within the triangle, once the point is confirmed If the triangle is inside, a random neighborhood sample point is immediately generated in the triangle using barycentric coordinates. The height of the sample point is aligned with the original surface point, and the sampling displacement is guaranteed not to exceed the preset threshold. If the mesh containing the current surface point P1 cannot be found in the queried triangle list, the neighborhood points are generated for point P1 by 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 sampling point in the 2D triangulated network corresponding to the surface point P1 are calculated respectively, and Y1 and Z1 are normalized to obtain the normal vectors Y2 and Z2 respectively. The average difference between Y2 and Z2 is used as the consistency measure, and the average difference is incorporated into the overall loss function to achieve smooth constraints on the normal extraction results in the network; after the normal vector is extracted, the normal vector difference between the target point and the neighboring points is calculated, and the consistency loss function is used to constrain the consistency of the normal vector; Step 5: Calculate the uncertainty-aware rendering loss; Step 6: Perform end-to-end joint optimization.

2. The method according to claim 1, characterized in that Step 1 includes: using a deep neural network to implicitly model the three-dimensional scene in the satellite image. The input of the deep neural network includes: spatial sampling point coordinates x, camera light direction d0, sun radiation direction s and time series 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 expressed by the following formula: （1）， in represents a deep neural network, represents forward propagation, ; ; ; represents real three-dimensional space; represents the unit sphere; represents the volume density of spatial sampling points; Represents the albedo color of the spatial sampling point in the observation direction d0; ∈[0,1] represents the weight of whether the spatial sampling point is in the shadow, and the value range is [0,1], where 0 represents complete direct lighting and 1 represents complete 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) : （2）， Where T(t) represents the transfer function; is the color of the sky ambient light; dt represents the integral 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; Indicates that the path parameter u∈[t n ,t]; du represents the integral variable.

3. The method according to claim 2, characterized in that Step 2 includes: calculating the volume density The gradient relative to the spatial sampling point coordinate x is used to obtain the normal vector of each surface point through back propagation. The gradient g(x) is calculated as follows: （4）， in Indicates that the gradient calculation is performed on the point with the spatial sampling point coordinate x; By the volume density of spatial sampling points Perform automatic differentiation to obtain the normal vector : （5）， in >0 is a constant to prevent the denominator from being zero.

4. The method according to claim 3, characterized in that Step 3 includes: Step 3-1, extract corner points; Harris operator is used to detect the response of satellite image data, and the threshold is set to obtain candidate corner points. In the Harris operator corner point detection, the response threshold T is set to Harris Set to the maximum response value H max 1% to 5% of the response, or directly set the 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. In each cluster, select the corner point with the highest response value. Finally, a set of feature corner points with uniform distribution and high response intensity are obtained as candidate nodes for constructing the triangulated network. Step 3-2, straight line segment detection and cleaning; Use the Line Segment Detector (LSD) algorithm to extract the most significant line segments in the image and clean up the initial line segment set as follows: Merge: If the angle between two line segments is less than 10° and the distance between their endpoints is less than 2 pixels, they are merged into a continuous line segment; Extension: If the distance between the endpoints of two line segments is less than 12 pixels, determine whether the two line segments are approaching along the extension direction. If so, add the intersection point along the extension direction of the two line segments; Alignment: When the Euclidean distance between the endpoints of two line segments that should intersect in the pixel coordinate system does not exceed 2 to 3 pixels, the endpoints are considered offset and the offset endpoints are adjusted so that all line segments intersect correctly at the true edge. Finally, we get a set of corner points and a set of cleaned line segments; Step 3-3, initialize the empty triangulated network; Using all cleaned line segments as hard constraint edges, the Delaunay triangulation (CDT) algorithm with constraints is called on the image plane. First, a non-overlapping triangulation is generated to ensure that all line segments are exactly as the triangulation boundaries. Step 3-4, corner point insertion and quality control; Based on the initialized triangulated mesh, try to insert the corner points extracted in step 3-1 into the mesh in sequence. 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 is the threshold; is the i-th interior angle of triangle T; i=1, 2, 3 correspond to the interior angle values ​​of the three vertices respectively; If both conditions are met, the corner point p is added to the grid as a constraint point; Steps 3-5, constraint update and iteration; Adjust the triangles that conflict with 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 eliminated; finally, a two-dimensional triangulated network is obtained.

5. The method according to claim 4, characterized in that In step 4, the consistency loss function is in the form of: （8）， Where s is the set of sampled surface point positions; Represents the normal vector at the spatial sampling point coordinate x; Indicates the coordinates of the sampling point in the two-dimensional triangulation corresponding to the surface point; Representing coordinates The normal vector at .

6. The method according to claim 5, characterized in that In step 5, the uncertainty-aware rendering loss is calculated using the following formula: : （9）， Among them E r represents the expectation of all rays r; Indicates along the ray Volume rendering weight at the kth sampling position; C k Indicates 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 kth sampling point.

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

8. An electronic device, characterized in that: The method comprises a processor and a memory, wherein the memory stores program codes, and when the program codes are executed by the processor, the processor is caused to perform the steps of the method according to any one of claims 1 to 7.

9. A storage medium, characterized in that: A computer program or instruction is stored, and when the computer program or instruction is run on a computer, the steps of the method according to any one of claims 1 to 7 are executed.

Images