A fast surface reconstruction method for remote sensing scenes

Abstract

The application discloses a fast surface reconstruction method for a remote sensing scene, relates to the technical field of remote sensing image processing, and aims at the problem that the existing surface reconstruction method needs additional mask input to ensure the quality of surface reconstruction, thereby causing low reconstruction efficiency. The application takes neural rendering as a basis, uses a signed distance field (SDF) as a logical representation of a 3D scene to ensure detailed surfaces, and adopts a voxel grid as an underlying data structure to ensure efficiency. Therefore, the application does not need additional mask input, can realize fast reconstruction of a three-dimensional curved surface representation of a scene, and can synthesize high-quality views from any view angle. In addition, the surface reconstruction result of the application can represent accurate geometric shapes.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image processing technology, specifically to a rapid surface reconstruction method for remote sensing scenes. Background Technology

[0002] 3D spatial modeling of natural scenes, and scene view synthesis based on spatial modeling priors, have been important directions for information technology efforts related to human interactive experience. For decades, researchers have been trying to transform real-world natural scenes into effective digital assets, but this is no easy task, especially for remote sensing images.

[0003] With the continuous advancement of drone and satellite technologies, the availability of airborne and satellite imagery has significantly improved, leading to a series of photogrammetric pipelines such as Pix4D, MicMac, and SURE. Although their imaging models differ, these traditional pipelines generally follow a common methodology. Given camera parameters, they first perform stereo matching, computing a dense point cloud from the input image. Subsequently, surface reconstruction is performed based on this point cloud. However, over-reliance on the generated point cloud often results in surface defects or noise issues.

[0004] Compared to traditional methods based on explicit representations (such as voxels, point clouds, and meshes), recent research in the field of neural rendering has focused on using implicit functions to represent scenes, known as Neural Radiation Fields (NeRF). NeRF can directly learn continuous 3D geometry from 2D images and has achieved excellent performance in view synthesis tasks. Currently, some researchers have attempted to apply NeRF technology to airborne and satellite imagery. While these attempts demonstrate the potential advantages of NeRF in remote sensing, they generally face a common problem: learning all geometric and color details based on a complex MLP framework leads to excessively long training times. On a single GPU device, it typically takes several hours to optimize each scene, or even a whole day for pre-training.

[0005] Existing NeRF training acceleration methods employ voxel meshes and small MLPs to store scene data, or use symbolic distance fields (SDF) to model the scene. However, their 3D surface reconstruction results cannot represent accurate geometry, exhibit significant noise and holes, and require additional mask inputs to ensure the quality of surface reconstruction, leading to low reconstruction efficiency. Summary of the Invention

[0006] The purpose of this invention is to address the problem that existing surface reconstruction methods require additional mask input to ensure the quality of surface reconstruction, which leads to low reconstruction efficiency. This invention proposes a fast surface reconstruction method for remote sensing scenarios.

[0007] The technical solution adopted by the present invention to solve the above-mentioned technical problems is as follows:

[0008] A fast surface reconstruction method for remote sensing scenes includes the following steps:

[0009] Step 1: Determine the reconstruction area. Then, acquire remote sensing images of the reconstruction area and process them using the Structure in Motion (SIM) algorithm to obtain camera intrinsic and extrinsic parameters and camera position. Finally, process the remote sensing images using a multi-view stereo matching algorithm to obtain the reconstructed point cloud p. rec ;

[0010] Step 2: Using the reconstructed point cloud p rec The bounding box of the reconstructed region is constructed based on the camera's shooting position, and two voxel meshes V are initialized within the reconstructed region. sdf With V feat The voxel grid V sdf The voxel grid V is used to store the symbolic distance function (SDF) values. feat Used to store feature values;

[0011] Step 3: For the voxel mesh V sdf Set the SDF value of the middle layer of the voxel mesh to 0, and set the upper layer of the two adjacent layers of the middle layer to positive and the lower layer to negative.

[0012] Step 4: Select any pixel P in the remote sensing image. ixel A ray p(t,v) is constructed using the camera's intrinsic and extrinsic parameters, and the ray p(t,v) is used to represent the pixel P. ixel After penetration, the ray p(t,v) is sampled to obtain pixel P. ixel The corresponding set of sampling points P = {p(t)} i ,v)||i=0,1,2,…,N-1}, where N is the total number of sampling points, i represents the i-th sampling point, and the ray p(t,v) is a three-dimensional vector function, i.e. p(t,v)=o+tv, where o represents the spatial coordinates of the camera, v represents the direction of the ray, and t represents the depth of the ray;

[0013] Step 5: Filter the sampling points in the sampling point set P, such that the sampling points t in the sampling point set P are... i Satisfy P 1 ={p(t) i )|t near ≤t i ≤t far}, where t near With t far These represent the upper and lower bounds of the bounding box of the reconstructed region, respectively. Then, for the set of sampling points P... 1 Sampling point p(t)i ), such that 0 < Ψ s (f(p(t i Since )))<1, we obtain the set of sampling points P. 2 , among which, Ψ s This represents the cumulative opacity, and f(·) represents the SDF value calculation function;

[0014] Step Six: Targeting P 2 For each sampling point, its position p and camera pose v are obtained. Then, the SDF value of the sampling point is obtained using trilinear interpolation. and eigenvalues Next, calculate the normal vector value of the sampling point location p. Then, p, v, The input is fed into MLP neural network 1 to obtain the initial color c of the sampling point. 0 Finally, p, v, c 0 , The input is fed into MLP neural network 2 to obtain the final color value c of the sampling point, where,

[0015] Step 7: Repeat step 6 to obtain P. 2 The final color value c corresponding to all sampling points in the sample, {c(t) i ,v)|i=0,1,2,…,K}, where K≤N;

[0016] Step 8: For the sampling point set P 2 Weight estimation is performed to obtain the sampling point set P. 2 The color weight w(t) of each sampling point i Then, based on the color weight w(t) of each sampling point... i The pixel P is obtained by weighted summation of the final color value c and the color value c. ixel Corresponding reconstructed color C volume_rendering ;

[0017] Step 9: Utilize the reconstructed color C volume_rendering With c 0 Construct the color loss function L rgb This leads to the total loss function L. total L total Represented as:

[0018]

[0019] Among them, L z Indicates the Z-axis limit, L pc L represents the point cloud constraint. ws Indicates weight limit, V represents sdfThe regularization parameter in the text, V represents feat The regularization parameter in the text, V represents sdf The gradient regularization parameter, λ z L represents z The weight, λ pc L represents pc The weight, λ ws L represents ws The weight, λ tv Represents the regularization parameter and The weight, λ grad Represents the gradient regularization parameter The weight.

[0020] Step 10: Repeat steps 4 through 9 until the total loss function L... total When the preset requirements are met, the optimized V is obtained. sdf This means completing surface reconstruction.

[0021] Furthermore, the bounding box of the reconstructed region includes airborne image bounding boxes and satellite image bounding boxes.

[0022] Furthermore, the steps for obtaining the airborne image bounding box are as follows:

[0023] First, reconstruct the point cloud p rec Outlier filtering is performed to obtain the ROI of the point cloud region. points Next, the geocentric coordinate system is transformed into a station-centered coordinate system. Then, the camera's shooting position is projected onto the xy plane of the station-centered coordinate system to obtain the minimum bounding rectangle, which is defined as the camera's location region (ROI). cam Then, rotate the entire space along the z-axis to make the edge of the smallest bounding rectangle parallel to the x-axis or y-axis. Finally, based on the point cloud region ROI... points and camera location area ROI cam And by using linear interpolation, the final accurate point cloud region ROI is obtained. final ROI final Represented as:

[0024] ROI final =min(ROI) cam ROI points )+λ|ROI points -ROI cam |

[0025] Among them, ROI cam The Region of Interest (ROI) represents the area containing all camera locations. points This represents the region containing all filter points, where λ represents the ROI.points and ROI cam The proportionality coefficient between them is λ = 0.5.

[0026] Furthermore, the size of each voxel in the airborne image bounding box is represented as follows:

[0027]

[0028] Among them, D x D y D z The bounding box represents the length along the x, y, and z axes, and M represents the total number of voxels.

[0029] Furthermore, the steps for obtaining the satellite image bounding box are as follows:

[0030] The projection from geographic coordinates to image coordinates is obtained using the RPC model. The projection from geographic coordinates to image coordinates is represented as follows:

[0031] P:P(lat,lon,alt)=(row,line)

[0032] Where (lat,lon,alt) represent longitude, latitude, and altitude in geographic coordinates, respectively, and (row,line) represent row and column in image coordinates;

[0033] The inverse transformation of the projection from geographic coordinates to image coordinates is expressed as:

[0034] Pinv:Pinv(row,line,alt)=(lat,lon)

[0035] Obtain the maximum altitude from remote sensing images (alt) min with minimum value alt max After that, based on alt min Minimum value alt max Using the inverse transformation formula with (row, line), we can obtain the latitude and longitude coordinates (lat, lone, alt) corresponding to the extreme altitude. min ) and (lat,lon,alt) max This process yields a set of geographic points, i.e., the bounding box of the satellite image, which is represented as:

[0036] P = (lat,lon,alt) min )∪(lat,lon,alt max ).

[0037] Furthermore, the set of sampling points P 2 Each sampling point is represented as:

[0038]

[0039]

[0040] Where Ω represents the interior space of the boundary region, S represents the region boundary, x represents the coordinates of the previous sampling point on the ray, y represents the coordinates of the sampling point on the boundary S, and I Ω f(x) represents the sign indicator function that determines whether the sampling point x is in space Ω, f(x) represents the SDF value corresponding to the sampling point x, and T and Ψ s This represents the cumulative opacity along the light ray.

[0041] Furthermore, the color weight w(t) i ) is represented as:

[0042]

[0043]

[0044]

[0045] Where σ represents the volume density of the sampling points, a i This represents the discrete opacity of the sampling points. This represents the product of the normal to the sampling point and the ray direction.

[0046] Furthermore, the reconstructed color C volume_rendering Represented as:

[0047]

[0048] Furthermore, the L z Represented as:

[0049]

[0050]

[0051] Where l, i, j, k represent the positional information in voxel space. N x N y N z V sdf The number of voxels in the x, y, and z axes, N = C × N x ×N y ×N z C represents the number of channels, C=1;

[0052] L pc Represented as:

[0053]

[0054] Where, p rec This represents the SDF value of the reconstructed sparse point f(p) at point p;

[0055] L ws Represented as:

[0056]

[0057] Where m represents the number of rays sampled in a batch, n represents the number of sampling points along the ray, and w i Indicates the color weight of the sampling points;

[0058] as well as The steps to obtain it are as follows:

[0059]

[0060]

[0061] Among them, L smooth With L tv For smoothness regularization;

[0062] L tv Applying this to Vsdf and Vfeat yields... and

[0063] Using V sdf Gradient mesh is obtained Then L smooth Applied to get

[0064] Furthermore, the L rgb Represented as:

[0065]

[0066] in, This represents the final true color value.

[0067] The beneficial effects of this invention are:

[0068] This application is based on neural rendering, utilizing signed distance fields (SDF) as the logical representation of 3D scenes to ensure detailed surfaces, while employing voxel meshes as the underlying data structure to ensure efficiency. Therefore, this application does not require additional input masks, enabling rapid reconstruction of 3D surface representations of scenes and the synthesis of high-quality views from arbitrary perspectives. Furthermore, the surface reconstruction results of this application can represent precise geometric shapes. Attached Figure Description

[0069] Figure 1 This is a flowchart illustrating the process of reconstructing airborne and satellite imagery for the framework of this application;

[0070] Figure 2 A schematic diagram of the bounding box initialization strategy for airborne and satellite imagery;

[0071] Figure 3 SDF distribution diagrams corresponding to spherical initialization and planar initialization;

[0072] Figure 4 A curve of sampling points;

[0073] Figure 5 Before and after restoration images showing the elimination of floating objects in space;

[0074] Figure 6 This is a schematic diagram of the reconstruction results on the LENIV-NVS dataset;

[0075] Figure 7 A visualization comparing the MAE (Mathematical Analysis) between the true depth and the estimated depth;

[0076] Figure 8 Comparison of surface reconstruction results on the BlendedMVS dataset;

[0077] Figure 9 This is a result of surface reconstruction from a satellite image;

[0078] Figure 10 A visualization comparing MAE and MRA between the actual altitude and the estimated altitude. Detailed Implementation

[0079] It should be noted that, where there is no conflict, the various embodiments disclosed in this application can be combined with each other.

[0080] Specific implementation method one: Refer to Figure 1 This embodiment specifically describes a rapid surface reconstruction method for remote sensing scenes, comprising:

[0081] Step 1: Use the structure-on-motion algorithm to recover the camera intrinsic and extrinsic parameters and camera position for each image from the remote sensing images, and then use the multi-view stereo matching algorithm to reconstruct the point cloud p from the remote sensing images. rec ;

[0082] Step 2: Using the reconstructed point cloud p rec The camera's shooting position is used to construct the boundary of the reconstructed region, and boundary initialization is completed. Two voxel meshes V are initialized within the reconstructed region. sdf With V featThese are used to store the symbolic distance function (SDF) value and the eigenvalue, respectively;

[0083] Step 3: Apply Voxel mesh V sdf Perform SDF initialization by setting the SDF value of the middle layer of the voxel mesh to 0, the upper part of the middle layer to positive, and the lower part to negative;

[0084] Step 4: Randomly select a remote sensing image, and randomly choose a pixel P in the image. ixel Using the camera's intrinsic and extrinsic parameters, a ray p(t,v) is constructed for pixel P. ixel The ray p(t,v) is then penetrated, and N points are sampled to obtain pixel P. ixel The corresponding set of sampling points P = {p(t)} i ,v)||i=0,1,2,…,N-1}, where the ray p(t,v) is a three-dimensional vector function, i.e. p(t,v)=o+tv, where o represents the spatial coordinates of the camera, v represents the direction of the ray, and t represents the depth of the ray;

[0085] Step 5: Filter the sampling points so that the set of sampling points satisfies P. 1 ={p(t) i )||t near ≤t i ≤t far}, where t near With t far These represent the upper and lower boundaries of the reconstructed region, respectively, and the sampling point set P is taken. 1 Point p(t) in i Make it satisfy 0 < Ψ s (f(p(t i )))<1,Ψ s Representing the cumulative opacity, the second filtering is completed, resulting in the set of sampling points P. 2 ;

[0086] Step Six: Randomly obtain P 2 Given a sampling point p with camera pose v, calculate the SDF value and features of this point using trilinear interpolation. And calculate the normal vector value of p. p, v, The input is fed into an MLP1 neural network to obtain the initial color c of the sampling point. 0 p, v, c 0 , The input is fed into an MLP2 neural network to obtain the final color c of the sampling point. This process is repeated to obtain the color value {c(t)} of the sampling point corresponding to the ray. i ,v)|i=0,1,2,…,K}, where K≤N;

[0087] Step 7: For the sampling point set P 2 Weight estimation is performed to obtain the sampling point set P. 2 The corresponding color weight {w(t)} i The color values ​​are weighted and summed according to their color weights (i = 0, 1, 2, ..., K), and the sum is used as the pixel P. ixel Corresponding reconstructed color C volume_rendering ;

[0088] Step 8: Calculate the Z-axis limit L z Point cloud limitation L pc Weight Limitation L ws , calculate SDF voxel mesh V sdf Feature Voxel Mesh V feat Regularization parameters in and SDF voxel grid V sdf gradient regularization parameter in

[0089] Step 9: Utilize the reconstructed color C volume_rendering With c 0 Construct the color loss function L rgb Construct the total loss function Where λ represents the weight of the corresponding loss term, and λ represents the ROI. points and ROI cam The scaling factor between two regions, with a value between 0 and 1. For example, when the result is 0, the ROI is... final =min(ROI) cam ROI points When the result is 1, the ROI is... final =max(ROI) cam ROI points To take into account the area of ​​both regions, this paper uses a value of 0.5.

[0090] Step 10: Repeat steps 4 through 9 until the total loss function L... total When the preset requirements are met, the optimized V is obtained. sdf This means completing surface reconstruction.

[0091] Given a set of multi-view images of a remote sensing scene as input (airborne or satellite images), this application aims to rapidly reconstruct a 3D surface representation of the scene without requiring additional input masks, and to synthesize high-quality views from any perspective. This invention is based on neural rendering, utilizing signed distance fields (SDF) as the logical representation of the 3D scene to ensure detailed surfaces, while employing voxel meshes as the underlying data structure to ensure efficiency. This invention proposes an SDF distributed initialization strategy and several geometric constraint mechanisms customized for scene priors. Furthermore, to adapt to remote sensing scenes, this application proposes a method for merging pinhole camera-based and rational polynomial camera (RPC)-based rendering into a single framework using different spatial initialization strategies. This application utilizes SDF as the logical representation of the 3D scene to ensure surface detail. Simultaneously, it employs voxel meshes as the underlying data structure to ensure efficiency. To simultaneously process airborne and satellite images, we merge pinhole camera-based and rational polynomial camera-based rendering into a single framework and employ different spatial initialization strategies. In addition, we introduce an SDF distributed initialization strategy and multiple geometric constraint prior mechanisms to enhance the adaptability of this invention in remote sensing environments.

[0092] Step 1 implementation: This invention uses the COLMAP algorithm to recover the camera pose and reconstruct the point cloud.

[0093] The specific implementation of step two: The main step is to establish the boundary region containing the space and adopt the SDF plane initialization strategy to facilitate the modeling of the space using voxel meshes.

[0094] 1. Region boundary initialization

[0095] Assuming the bounding box is aligned with Cartesian coordinates, D x D y D z M is the length of the bounding box along the x, y, and z axes, and M is the total number of voxels, with each voxel having the following size:

[0096]

[0097] Specific steps for initializing airborne image scene bounding boxes using sparse point clouds:

[0098] 1) Apply statistical filtering to the reconstructed points to filter out most outliers.

[0099] 2) Convert ECEF coordinates to ENU coordinates using a 4×4 matrix, as ENU coordinates have the advantage of automatically aligning the top of the scene with the z-axis.

[0100] 3) The minimum bounding rectangle (MER) is calculated by projecting the camera position onto the xy plane.

[0101] 4) Rotate the entire space along the z-axis so that the calculated edge of MER is parallel to the x-axis or y-axis.

[0102] 5) Determine the precise regions of interest in the point cloud. We now have two regions of interest (ROIs): ROI cam , representing the region containing all camera locations, ROI points This represents the region containing all filter points. The final region is calculated using linear interpolation:

[0103] ROI final =min(ROI) ca, ROI points )+λ|ROI points -ROI cam |#(2)

[0104] λ is a value between 0 and 1. When λ = 0, the region ROI is... final =min(ROI) cam ROI points When λ = 1, the region ROI final =max(ROI) cam ROI points ).

[0105] 6) Select ROI final The points in the bounding box are calculated.

[0106] Satellite imagery is typically described by an RPC model. Each RPC model defines a bidirectional projection between geographic coordinates and image coordinates. We call the projection from geographic coordinates to image coordinates...

[0107] P:P(lat,lon,alt)=(row,line)#(3)

[0108] Its inverse transformation is

[0109] Pinv:Pinv(row,line,alt)=(lat,lon)#(4)

[0110] We first apply Pinv to alt. min and alt max The appropriate bounding box of the scene is calculated for all image coordinates {(row,line)} within the defined elevation range. This involves using two rectangular planes to delineate the sampling range of the scene, and this process generates a set of geographic points.

[0111] P = (lat,lon,alt) min )∪(lat,lon,alt max )#(5)

[0112] This process is as follows Figure 2 As shown in part a.

[0113] 2. SDF plane initialization strategy:

[0114] This invention uses SDF (Surface Depositional Focal Array) to model the geometry of a 3D scene; therefore, given a bounding box, it is necessary to specify the initial distribution of the SDF within it. Unlike previous methods that support spherical initialization, this invention proposes a planar SDF initialization strategy. This strategy addresses the convergence difficulties and shape unpredictability associated with spherical initialization, and is more suitable for fitting remote sensing data due to its closer alignment with terrain object morphology. This process achieves the following results: Figure 3 As shown.

[0115] Step 3 Implementation: The purpose of this step is to optimize the implicit representation of the 3D scene. The entire optimization process is as follows: Figure 2 As shown. Since this is the core part of the invention, for the sake of clarity, we will describe the specific implementation process of step three in four parts:

[0116] 1. Scene Representation Methods

[0117] This invention uses a signed distance field (SDF) f to model the entire space (including occupied and unoccupied areas). Assume S represents the boundary of the occupied space, i.e., the Earth's surface. I represents the space occupied by objects such as trees and buildings. Ω (·) is an exponential function of the space Ω. This distance field represents the spatial position. Mapped to the scale value as follows:

[0118]

[0119] This invention uses SDF functions to represent volume rendering. Given a pixel, assume the ray from the camera to this pixel is {p(t) = o + tv | t ≥ 0}, where o is the camera position and v is the unit direction vector of the ray. Therefore, the weighted sum of the volume rendering colors can be expressed as:

[0120]

[0121] Among them, T i a represents discrete cumulative transmittance. i Represents discrete opacity, c i This indicates a color estimate.

[0122] This invention will accumulate transmission T i Modeled as a function of SDF, where Ψ s It's the Sigmoid function, where s controls Ψ. sTrainable parameters of the slope:

[0123]

[0124] Using Equation 8, the volume density σ and discrete opacity a i This can be represented using SDF as follows:

[0125]

[0126]

[0127] in, This indicates that at position p(t) i The product of the normal at point 0 and the direction of the ray.

[0128] 2. Adaptive point cloud filtering:

[0129] This invention proposes a novel point cloud filtering method, which first filters point light sampling points {p(t)}. i )|t near ≤t i ≤t far For points outside the specified range, the second step is to calculate the cumulative transmittance in Formula 8, filtering out points whose transmittance is outside the range shown in Formula 11.

[0130] 0<Ψ s (f(p(t i )))<1#(11)

[0131] Formula 8 models transmittance as a sigmoid function of SDF value, with parameter s -1 Control the steepness of its slope near zero. As the network converges, s -1 As the value approaches 0, the slope of the sigmoid function gradually increases near the zero point. This causes the number of sampling points satisfying Equation 11 to decrease over time, a process as follows: Figure 4 As shown.

[0132] 3. Hybrid network architecture:

[0133] To accelerate training, this invention uses a voxel-based architecture to encode the SDF function, employing two voxel networks V. feat and V sdf This is used to store the corresponding features and SDF values ​​for a given space. This voxel-based architecture allows for efficient querying of features or SDF at any location p via interpolation:

[0134] feat(p) = interp(p,V) feat )#(12)

[0135] sdf(p)=interp(p,V sdf )#(13)

[0136] To ensure consistent geometry while simultaneously recovering accurate details, we employed a two-stage training scheme. In the coarse stage, the color of the sampled points is predicted by a 3-layer MLP:

[0137] c(p,v)=MLP c (p,v,n,feat(p))#(14)

[0138] Where v is the viewing direction and n is the normal calculated based on the gradient of the SDF value.

[0139] In the refinement stage, the color of the sampling points is predicted using two 4-layer MLPs combined in a residual manner. Specifically, First use V sdf Calculate the hierarchical geometric features to predict color c 0 (p,v). Subsequently, MLP f c 0 Using (p,v) as input, predict the final color with fine details:

[0140] c(p,v)=MLP f (p,v,n,feat(p),c 0 (p,v))

[0141]

[0142]

[0143]

[0144] in and These are the normals and SDF values ​​calculated from adjacent positions of p, with the superscript i indicating the distance between the position and p.

[0145] 4. Geometric constraints based on scene priors

[0146] The constraints provided by color loss alone are loose geometric constraints and cannot guarantee the consistency of all views in the direction of geometric optimization. In order to promote the correct fitting of the scene by the SDF and derive a high-quality reconstructed surface, this invention introduces several regularization constraints.

[0147] Z-axis constraint:

[0148] After data preprocessing, the ground plane normal of the scene is aligned with the positive Z-axis of the coordinate system. Based on this, for any line parallel to the Z-axis, we expect that the line intersects the boundary of the occupied space S at only one point, and that the SDF value should gradually decrease as the line moves in the negative Z-axis direction. To describe this expectation, we propose a Z-axis constraint L. z :

[0149]

[0150] in, N x N y N z They are V sdf Number of voxels in the x, y, and z axes:

[0151]

[0152] At the same time, N = C × N x ×N y ×N z C represents the number of channels, for V sdf C = 1. L z Eliminating uncontrollable floating objects distributed in space facilitates the restoration of the building's vertical structure, such as... Figure 5 As shown.

[0153] Point constraints:

[0154] The SFM algorithm can reconstruct sparse points in a scene during camera parameter estimation. Although the points reconstructed by SFM are inevitably affected by noise, they still possess a certain level of accuracy. Furthermore, considering our end-to-end training framework, noise can be corrected to some extent during optimization. Therefore, we use these reconstructed sparse points p rec Additional geometric information is used to further constrain the SDF field:

[0155]

[0156] For satellite imagery, this invention employs bundle adjustment to perform relative calibration on the RPC model of all input images before retraining begins. Then, we use the points generated during the calibration process to constrain the SDF field.

[0157] Weight constraints:

[0158] For remote sensing scenes where camera light rays are emitted from top to bottom, they will inevitably hit objects such as the ground or buildings. This invention uses the weighted sum of light rays to characterize this, assuming that a batch of rays must penetrate the object:

[0159]

[0160] Where m represents the number of rays sampled in a batch, and n represents the number of sampling points along the ray. L ws It helps to eliminate holes distributed on the surface, and also helps to fit the correct geometry.

[0161] Loss function:

[0162] This invention trains vox2surf by minimizing the difference between rendered colors and real colors, i.e., color constraint:

[0163]

[0164] To maintain the continuity and smoothness of the surface, this invention employs a smoothness regularization term L. tv and L smooth Used to ensure consistency between adjacent voxel meshes:

[0165]

[0166]

[0167] Among them, L smooth With L tv The smoothness regularization term is used to ensure consistency between adjacent voxel meshes.

[0168] Where N represents the total number of grid cells multiplied by the number of channels, and conv(V) represents applying a 3D convolution operation on grid cells V with a stride of 1, using a Gaussian kernel of size 3. We will L tv Applied to SDF mesh V sdf and feature grid V feat and represent them as and In addition, we will L smooth Applied to from V sdf Calculated gradient grid and represent it as The final loss function is:

[0169]

[0170] Step 4 Implementation: Following the optimization in Step 3, the SDF network V sdf With feature network V featThe geometry and appearance of the target have been captured. At this point, to export the mesh model of the scene, simply perform uniform sampling in cubic space; then, use the SDF network to determine the SDF value of each sampling point's coordinates; finally, use the MarchingCubes algorithm to reconstruct the mesh model based on the coordinates and SDF values ​​of the sampling points. The entire process is as follows: Figure 1 As shown.

[0171] To generate a new view, simply provide the new viewpoint direction, and then perform forward inference between the SDF network and the color network. The result of the forward inference is the view synthesized from the new viewpoint. The entire process is as follows: Figure 1 As shown.

[0172] It should be noted that the specific embodiments are merely explanations and illustrations of the technical solution of the present invention and should not be used to limit the scope of protection. Any modifications made in accordance with the claims and specification of the present invention that are only partial should still fall within the protection scope of the present invention.

Claims

1. A rapid surface reconstruction method for remote sensing scenes, characterized in that... Includes the following steps: Step 1: Determine the reconstruction area. Then, acquire remote sensing images of the reconstruction area and process them using the Structure in Motion (SIM) algorithm to obtain camera intrinsic and extrinsic parameters and camera position. Finally, process the remote sensing images using a multi-view stereo matching algorithm to obtain the reconstructed point cloud p. rec ; Step 2: Using the reconstructed point cloud p rec The bounding box of the reconstructed region is constructed based on the camera's shooting position, and two voxel meshes V are initialized within the reconstructed region. sdf With V feat The voxel grid V sdf The voxel grid V is used to store the symbolic distance function (SDF) values. feat Used to store feature values; Step 3: For the voxel mesh V sdf Set the SDF value of the middle layer of the voxel mesh to 0, and set the upper layer of the two adjacent layers of the middle layer to positive and the lower layer to negative. Step 4: Select any pixel P in the remote sensing image. ixel A ray p(t, v) is constructed using the camera's intrinsic and extrinsic parameters, and the ray p(t, v) is used to represent the pixel P. ixel After penetration, the ray p(t, v) is sampled to obtain pixel P. ixel The corresponding set of sampling points P = {p(t)} i ,v)|i=0,1,2,...,N-1},where N is the total number of sampling points, i represents the i-th sampling point, and the ray p(t,v) is a three-dimensional vector function, i.e. p(t,v)=o+tv, where o represents the spatial coordinates of the camera, v represents the direction of the ray, and t represents the depth of the ray; Step 5: Filter the sampling points in the sampling point set P, such that the sampling points t in the sampling point set P are... i Satisfy P 1 ={p(t) i )|t near ≤t i ≤t far }, where t near With t far These represent the upper and lower bounds of the bounding box of the reconstructed region, respectively. Then, for the set of sampling points P... 1 Sampling point p(t) i ), such that 0 < Ψ s (f(p(t i Since )))<1, we obtain the set of sampling points P. 2 , among which, Ψ s This represents the cumulative opacity, and f(·) represents the SDF value calculation function; Step Six: Targeting P 2 For each sampling point, its position p and camera pose v are obtained. Then, the SDF value of the sampling point is obtained using trilinear interpolation. and eigenvalues Next, calculate the normal vector value of the sampling point location p. Then, p, v, The input is fed into MLP neural network 1 to obtain the initial color c of the sampling point. 0 Finally, p, v, c 0 , The input is fed into MLP neural network 2 to obtain the final color value c of the sampling point, where, Step 7: Repeat step 6 to obtain P. 2 The final color value c corresponding to all sampling points in the sample, {c(t) i ,v)|i=0,1,2,...,K},whereK≤N; Step 8: For the sampling point set P 2 Weight estimation is performed to obtain the sampling point set P. 2 The color weight w(t) of each sampling point i Then, based on the color weight w(t) of each sampling point... i The pixel P is obtained by weighted summation of the final color value c and the color value c. ixel Corresponding reconstructed color C volume_rendering ; Step 9: Utilize the reconstructed color C volum_erendering With c 0 Construct the color loss function L rgb This leads to the total loss function L. total L total Represented as: Among them, L z Indicates the Z-axis limit, L pc L represents the point cloud constraint. ws Indicates weight limit, V represents sdf The regularization parameter in the text, V represents feat The regularization parameter in the text, V represents sdf The gradient regularization parameter, λ z L represents z The weight, λ pc L represents pc The weight, λ ws L represents ws The weight, λ tv Represents the regularization parameter and The weight, λ grad Represents the gradient regularization parameter The weights; Step 10: Repeat steps 4 through 9 until the total loss function L... total When the preset requirements are met, the optimized V is obtained. sdf This means completing surface reconstruction.

2. The rapid surface reconstruction method for remote sensing scenes according to claim 1, characterized in that... The bounding boxes of the reconstructed region include airborne image bounding boxes and satellite image bounding boxes.

3. The rapid surface reconstruction method for remote sensing scenes according to claim 2, characterized in that... The steps for obtaining the airborne image bounding box are as follows: First, reconstruct the point cloud p rec Outlier filtering is performed to obtain the ROI of the point cloud region. points Next, the geocentric coordinate system is transformed into a station-centered coordinate system. Then, the camera's shooting position is projected onto the xy plane of the station-centered coordinate system to obtain the minimum bounding rectangle, which is defined as the camera's location region (ROI). cam Then, rotate the entire space along the z-axis to make the edge of the smallest bounding rectangle parallel to the x-axis or y-axis. Finally, based on the point cloud region ROI... points and camera location area ROI cam And by using linear interpolation, the final accurate point cloud region ROI is obtained. final ROI final Represented as: KING final <min(ROI cam ,KING points )+λ|ROI points -KING cam | Among them, ROI cam The Region of Interest (ROI) represents the area containing all camera locations. points This represents the region containing all filter points, where λ represents the ROI. points and ROI cam The proportionality coefficient between them is λ = 0.

5.

4. The rapid surface reconstruction method for remote sensing scenes according to claim 3, characterized in that... The size of each voxel in the airborne image bounding box is represented as follows: Among them, D x D y D z The bounding box represents the length along the x, y, and z axes, and M represents the total number of voxels.

5. A rapid surface reconstruction method for remote sensing scenes according to claim 4, characterized in that... The steps for obtaining the satellite image bounding box are as follows: The projection from geographic coordinates to image coordinates is obtained using the RPC model. The projection from geographic coordinates to image coordinates is represented as follows: P:P(lat,lon,alt)=(row,line) Where (lat,lon,alt) represent longitude, latitude, and altitude in geographic coordinates, respectively, and (row,line) represent row and column in image coordinates; The inverse transformation of the projection from geographic coordinates to image coordinates is expressed as: Pinv:Pinv(row,line,alt)=(lat,lon) Obtain the maximum altitude from remote sensing images (alt) min with minimum value alt max After that, based on alt min Minimum value alt max Using the inverse transformation formula with (row, line), we can obtain the latitude and longitude coordinates (lat, lon, alar) corresponding to the extreme altitude. min ) and (lat,lon,alt) max This process yields a set of geographic points, i.e., the bounding box of the satellite image, which is represented as: P=(lat,lon,alt min )∪(lat,lon,alt max )。 6. A rapid surface reconstruction method for remote sensing scenes according to claim 5, characterized in that... The set of sampling points P 2 Each sampling point is represented as: Where Ω represents the interior space of the boundary region, S represents the region boundary, x represents the coordinates of the previous sampling point on the ray, y represents the coordinates of the sampling point on the boundary S, and I Ω f(x) represents the sign indicator function that determines whether the sampling point x is in space Ω, f(x) represents the SDF value corresponding to the sampling point x, and T and Ψ s This represents the cumulative opacity along the light ray.

7. A rapid surface reconstruction method for remote sensing scenes according to claim 6, characterized in that... The color weight w(t) i ) is represented as: Where σ represents the volume density of the sampling points, a i This represents the discrete opacity of the sampling points. This represents the product of the normal to the sampling point and the ray direction.

8. A rapid surface reconstruction method for remote sensing scenes according to claim 7, characterized in that... The reconstructed color C volume_rendering Represented as:

9. A rapid surface reconstruction method for remote sensing scenes according to claim 8, characterized in that... The L z Represented as: Where l, i, j, k represent the positional information in voxel space. N x N y N z V sdf The number of voxels in the x, y, and z axes, N = C × N x ×N y ×N z C represents the number of channels, C=1; L pc Represented as: Where, p rec This represents the SDF value of the reconstructed sparse point f(p) at point p; L ws Represented as: Where m represents the number of rays sampled in a batch, n represents the number of sampling points along the ray, and w i Indicates the color weight of the sampling points; as well as The steps to obtain it are as follows: Among them, L smooth With L tv For smoothness regularization; L tv Applied to V sdf and V feat ,get and Using V sdf Gradient mesh is obtained Then L smooth Applied to get 10. A rapid surface reconstruction method for remote sensing scenes according to claim 9, characterized in that... The L rgb Represented as: in, This represents the final true color value.

Images