A double-layer monte carlo optimized satellite photogrammetry terrain modeling method

Abstract

The application discloses a kind of double-layer Monte Carlo optimization satellite photogrammetry terrain modeling methods, it is related to remote sensing and photogrammetry technical field, first, satellite laser height measurement data of study area is collected, and original data is filtered, and the stereo image of study area is obtained to carry out preliminary control point layout and aerial triangulation encryption;Second, the first layer Monte Carlo sampling system is constructed, and high confidence sample points are screened;Third, the second layer Monte Carlo offset experiment is established, and the multi-direction disturbance experiment is generated by analyzing the geometric response relationship between laser points and stereo images Multiple sets of random offset samples;Finally, the optimal offset parameter is determined by error response analysis, and the final high-precision digital elevation model is output;The method makes full use of laser point reference, combined with Monte Carlo optimization and spatial offset sensitivity analysis, realizes the high-precision, robust construction of DEM, and significantly improves the reliability and precision of digital elevation model.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of remote sensing and photogrammetry, in particular to a satellite photogrammetry terrain modeling method based on double-layer Monte Carlo optimization. BACKGROUND

[0002] DEM (Digital Elevation Model) is an important basic data of geographic information system, and has a wide range of applications in the fields of terrain analysis, ecological environment monitoring, hydrological simulation, land use assessment and geological disaster prediction. The precision of DEM directly determines the reliability of subsequent spatial analysis and environmental modeling, so high-precision terrain modeling has always been the focus of research in photogrammetry and remote sensing. The traditional photogrammetry method based on satellite stereo image pairs relies on ground control points GCP for aerial triangulation and stereoscopic orientation, and the modeling accuracy depends largely on the spatial distribution and height accuracy of the control points. However, in mountainous, uninhabited or sparsely controlled areas, it is difficult to obtain high-quality control points, resulting in photogrammetry results being easily affected by terrain undulations, attitude errors and image pair matching errors, resulting in local distortion and error accumulation, thus limiting the overall accuracy of DEM.

[0003] The development of satellite laser altimetry technology provides high-precision ground elevation observation data on a global scale. For example, the ICESAT-2 laser ranging accuracy can reach centimeter level, which can provide reliable elevation control reference for stereo mapping in uncontrolled or sparsely controlled areas. However, ICESat-2 laser data is affected by factors such as surface reflectivity, atmospheric scattering, spot energy distribution and satellite attitude error, and the laser points have certain random errors and plane offsets on the ground.

[0004] Existing researches mostly use static threshold or empirical rules to screen and remove abnormal points from laser altimetry data. This kind of method can reduce the influence of abnormal data to a certain extent, but it is difficult to fully characterize the statistical characteristics of errors, and it cannot dynamically reflect the propagation law of control point errors in stereo modeling, and it also lacks a systematic correction mechanism for spatial offset of laser points, so that precision fluctuations and local error accumulation may still occur in complex terrain conditions. SUMMARY

[0005] In order to solve the above technical problems, the present application provides a satellite photogrammetry terrain modeling method based on double-layer Monte Carlo optimization, comprising the following steps:

[0006] S1, collecting satellite laser altimetry data in the study area, screening the original data according to the characteristics of satellite laser points, and generating a set of screened laser altimetry points;

[0007] S2, obtain stereo image pairs of remote sensing images of the study area, perform initial control point layout and aerial triangulation based on the laser height measurement point set obtained in the previous step, establish a stereo geometric model and generate a preliminary digital elevation model;

[0008] S3, based on the laser height measurement point set obtained in step S1, a first layer Monte Carlo sampling system is constructed, a fixed proportion of sample points are randomly selected in each round of experiment for re-estimation of geometric model parameters, the remaining samples are used for independent accuracy evaluation, and a preliminary optimization result is obtained through adjustment optimization;

[0009] S4, error variance and stability analysis is performed on the results of multiple rounds of experiments in the previous step, and a sample point set with a confidence level greater than 0.8 is selected as the reference point of the stereo image, and a first layer Monte Carlo optimized digital elevation model is generated based on this;

[0010] S5, using the sample point set selected in the previous step, a second layer Monte Carlo offset experiment is established, the offset range and step are determined by analyzing the geometric response relationship between the laser points and the stereo image, and are used for spatial position error sensitivity evaluation;

[0011] S6, under the offset parameter conditions determined in the previous step, a multi-direction perturbation experiment is performed with the laser point spot center as the center, and a plurality of groups of random offset sample coordinate sets are generated, each group of samples corresponds to a spatial geometric hypothesis;

[0012] S7, input each group of offset samples generated in the previous step into the stereo image geometric model, respectively perform forward intersection and parameter re-estimation, form a plurality of groups of digital elevation model samples, determine the optimal offset parameter through error response analysis, and output the final digital elevation model.

[0013] The further defined technical solution of the application is:

[0014] Further, in step S1, the ICESat-2 satellite ATL08 laser height measurement data of the study area is obtained, first, the original photon point cloud is converted into a standardized format, and abnormal points affected by atmospheric clouds, terrain shielding or photon noise are removed according to laser point attribute information, the laser point attribute information includes signal confidence level mark, ground object type identification, reflection energy distribution, echo shape and detection channel; Then select the top 30 ground laser points with the highest quality from the removed laser points for subsequent experiments.

[0015] The satellite photogrammetry terrain modeling method of double-layer Monte Carlo optimization as described above, in step S2, the high-resolution stereo image pair remote sensing image of the study area is obtained, first, the image is preprocessed, and the attitude deviation and light influence are eliminated by using the radiation calibration and geometric correction method; combined with the laser points selected in step S1, the image is punctured and aerial triangulation is performed, half of the number of laser points are used as control points to participate in the aerial triangulation encryption, and the remaining half of the number of points are used as independent check points; the geometric constraint relationship between the ground points and the image points is established by using the collinearity equation, and the initial solution of the exterior orientation elements is performed by the minimum projection residual method, and a preliminary digital elevation model is generated to provide a benchmark for subsequent optimization.

[0016] The satellite photogrammetry terrain modeling method of double-layer Monte Carlo optimization as described above, in step S3, a first layer Monte Carlo experiment is designed, half of the number of laser points selected from step S1 are randomly extracted as control points to participate in the geometric parameter reestimation and adjustment optimization, and the remaining half of the number of points are used as check points and do not participate in the adjustment, a total of 300 random extractions are performed, and the error of each point is calculated after each extraction and adjustment is completed, and the error calculation formula is as follows:

[0017] ;

[0018] Wherein, Xi represents the observation value of the laser point in the x direction, xi represents the true value of the laser point in the x direction; Yi represents the observation value of the laser point in the y direction, yi represents the true value of the laser point in the y direction; Zi represents the observation value of the laser point in the z direction, and zi represents the true value of the laser point in the z direction.

[0019] The satellite photogrammetry terrain modeling method of double-layer Monte Carlo optimization as described above, in step S3, the adjustment optimization method adopted is the minimum projection residual optimization, which specifically includes the following steps:

[0020] S3.1, the geometric relationship between the ground points and the image coordinates is established by using the collinearity equation, and the formula is as follows:

[0021] ;

[0022] Wherein, represents the ground point; represents the image coordinate; is the exterior orientation element of the satellite sensor, wherein represents the spatial position of the satellite sensor in the ground coordinate system, represents the attitude angle of the satellite sensor, that is, the rotation angle around the X, Y and Z axes respectively;

[0023] S3.2, for each laser point i, the projection residual is defined, and the formula is as follows:

[0024] ;

[0025] wherein, and denote the projection residual of each laser point i, and denote the real coordinates of laser points on the image, and are the theoretical coordinates of laser points on the image;

[0026] S3.3, taking the residual square sum as the optimization objective function, the formula is as follows:

[0027] ;

[0028] wherein, denotes the projection residual square sum objective function;

[0029] S3.4, the iterative least square method is used for parameter solving, the formula is as follows:

[0030] ;

[0031] wherein, denotes the correction amount of exterior orientation elements, and denote the parameter estimation value of the kth and k+1th iteration respectively, A denotes the partial derivative matrix, and L denotes the constant term vector.

[0032] A satellite photogrammetry terrain modeling method of double-layer Monte Carlo optimization as described above, in step S4, the laser point confidence calculation formula is as follows:

[0033] ;

[0034] ;

[0035] ;

[0036] wherein, denotes the average error of the laser point, denotes the error of the laser point after the jth Monte Carlo test, and m denotes the number of Monte Carlo tests, denotes the standard deviation of the point, denotes the confidence index of the point.

[0037] A satellite photogrammetry terrain modeling method of double-layer Monte Carlo optimization as described above, in step S5, the laser point offset range is determined by the following formula:

[0038] ;

[0039] wherein, represents an experience coefficient, represents a laser point spot diameter.

[0040] The double-layer Monte Carlo optimization satellite photogrammetry terrain modeling method as described above, in step S5, the offset amount of the laser point is determined according to the set maximum offset range M, and the calculation method is as follows:

[0041] ;

[0042] wherein, represents a laser point random offset amount, represents a laser point offset range, t is a grading coefficient, q represents an integer, and 1≤q≤10.

[0043] The double-layer Monte Carlo optimization satellite photogrammetry terrain modeling method as described above, in step S6, each offset sample coordinate is expressed as:

[0044] ;

[0045] wherein, represents a laser point offset sample coordinate, represents an original laser point coordinate; controls the offset direction, and corresponds to the X-axis and Y-axis directions respectively; and represents a random disturbance amount, and is subject to a normal distribution respectively:

[0046] ;

[0047] wherein, U(a,b) represents a uniformly distributed random number in the interval [a,b], and respectively represent the maximum disturbance amount in the X and Y directions.

[0048] The double-layer Monte Carlo optimization satellite photogrammetry terrain modeling method as described above, in step S7, an error response analysis is carried out by using an average absolute error analysis method, and the average absolute error calculation formula is as follows:

[0049] ;

[0050] wherein, MAE represents a laser point average absolute error, Z represents an elevation value of a laser point on a digital elevation model, z represents a true value of a laser point on a digital elevation model, and n represents the number of inspection points.

[0051] The beneficial effects of the present application are:

[0052] (1) In the present application, by establishing a double-layer Monte Carlo optimization system, the geometric response of laser points and stereo images is simulated in multiple directions and multiple scales, the systematic analysis of control point selection, spatial offset sensitivity and error response is realized, thereby effectively improving the accuracy and reliability of the digital elevation model;

[0053] (2) In the present application, not only can a stable and reliable elevation model be obtained under different laser point quality and distribution conditions, but also the adaptability of DEM to complex terrain and point cloud noise is enhanced by quantifying error response, and it has high repeatability and wide applicability, which can provide accurate and reliable data basis for various topographic mapping and geographic information analysis. BRIEF DESCRIPTION OF DRAWINGS

[0054] Figure 1 It is the overall flowchart of the present application;

[0055] Figure 2 It is the high-resolution seven satellite stereo image pair remote sensing image diagram of the sample area of the embodiment of the present application;

[0056] Figure 3 It is the digital elevation model DEM generated by ICESAT-2 after double-layer Monte Carlo optimization generated by the embodiment of the present application, wherein (a) is the DEM after the first layer Monte Carlo optimization, (b) is the DEM after the second layer Monte Carlo optimization, and the point number on the diagram is the check point, and the larger the shape, the greater the check point error;

[0057] Figure 4 It is the spot diagram of ICESAT-2 satellite laser point in the embodiment of the present application, wherein R is the spot radius, d is the spatial resolution of high-resolution seven stereo image pairs, and is 0.8 meters;

[0058] Figure 5 It is the offset diagram in the coordinate offset correction in the second layer Monte Carlo experiment in the embodiment of the present application, the green arrow indicates that the laser point coordinates are offset along the x direction with the spot center as the center, the red arrow indicates that the laser point coordinates are offset along the x direction with the spot center as the center, and the blue arrow indicates that the laser point coordinates are offset along the xy direction with the spot center as the center;

[0059] Figure 6 It is the comparison diagram of the best point position after offset correction of the present application and the point position matched with the stereo image of ICESAT-2 satellite laser point in the embodiment of the present application, wherein (a) is the forward-looking image of high-resolution seven stereo image pairs, and (b) is the rearview image of high-resolution seven stereo image pairs. DETAILED DESCRIPTION

[0060] The satellite photogrammetry terrain modeling method provided in the embodiment of the present application is a double-layer Monte Carlo optimization method, which comprises the following steps: Figure 1As shown, it includes the following steps:

[0061] S1. Obtain the ICESat-2 satellite ATL08 laser altimetry data of the study area. First, the original photon point cloud is standardized and converted. Based on the laser point attribute information (including signal confidence marker, ground object type identifier, reflected energy distribution, echo morphology, and detection channel), abnormal points affected by atmospheric clouds, terrain obstruction, or photon noise are removed. Then, the top 30 ground laser points with the highest quality are selected from the removed laser points for subsequent experiments.

[0062] S2, obtain as follows Figure 2 The Gaofen-7 stereo image pairs of the study area shown were first preprocessed, and radiometric calibration and geometric correction methods were used to eliminate attitude deviations and illumination effects. Combined with the 30 laser points selected in step S1, point-spotting and aerial triangulation were performed on the image, with 15 points used as control points for aerial triangulation refinement and the other 15 points used as independent check points.

[0063] The geometric constraint relationship between ground points and image points is established by using collinearity equations, and the exterior orientation elements are initially solved by the minimum projection residual method to generate a preliminary DEM (Digital Elevation Model), which provides a benchmark for subsequent optimization.

[0064] S3. Based on the 30 high-quality laser points selected in step S1, design the first layer of Monte Carlo experiments. Randomly select 50%, i.e., 15 laser points, as control points to participate in geometric parameter reestimation and adjustment optimization. The remaining 15 points are used as check points and do not participate in adjustment. A total of 300 random samplings are performed. After each sampling and adjustment, the error of each point is calculated. The error calculation formula is as follows:

[0065] ;

[0066] Where Xi represents the observed value of the laser point in the x-direction, xi represents the actual value of the laser point in the x-direction; Yi represents the observed value of the laser point in the y-direction, yi represents the actual value of the laser point in the y-direction; Zi represents the observed value of the laser point in the z-direction, zi represents the actual value of the laser point in the z-direction.

[0067] The adjustment optimization method is minimum projection residual optimization, which includes the following steps:

[0068] S3.1. Establish the geometric relationship between ground points and image coordinates using the collinearity equation, as shown in the following formula:

[0069] ;

[0070] in, Indicates a ground point; represents the image coordinates; is the exterior orientation elements of satellite sensor, wherein represents the spatial position of satellite sensor in the ground coordinate system, represents the attitude angle of satellite sensor (i.e. the rotation angle around X, Y, Z axis respectively).

[0071] S3.2, for each laser point i, define the projection residual, the formula is as follows:

[0072] ;

[0073] wherein, and represent the projection residual of each laser point i, and represent the real coordinates of the laser point on the image, and are the theoretical coordinates of the laser point on the image.

[0074] S3.3, again taking the residual sum of squares as the optimization objective function, the formula is as follows:

[0075] ;

[0076] wherein, represents the projection residual sum of squares objective function.

[0077] S3.4, the iterative least squares method is used to solve the parameters, the formula is as follows:

[0078] ;

[0079] wherein, represents the correction amount of exterior orientation elements, and respectively represent the parameter estimation value of the kth and k+1th iteration, A represents the partial derivative matrix, and L represents the constant term vector.

[0080] S4, after 300 iterations, the average error of each point is calculated according to the error of each point calculated in step S3, then the standard deviation of each point is calculated, finally the confidence index of each point is calculated, and each point is arranged in descending order according to the confidence, and the laser points with confidence greater than 0.8 are selected as control points to prick on the stereo image pair, then the first layer of Monte Carlo optimized digital elevation model is generated as shown in Figure 3 (a), the confidence index of each point is arranged as shown in Table 1, and the greater the confidence index, the higher the accuracy of the point, and the confidence of each point is calculated according to the following formula:

[0081] ;

[0082] ;

[0083] ;

[0084] wherein, represents the average error of the laser point, represents the error of the laser point after the jth Monte Carlo test, m represents the number of Monte Carlo tests, represents the standard deviation of the point, represents the confidence index of the point.

[0085] Table 1 Confidence index ranking table of each point

[0086]

[0087] S5, since the ICESat-2 satellite laser altimetry system adopts a photon counting ranging principle, the spot diameter formed by the laser beam on the ground is about 17.5 meters, and the spot diagram is as shown in Figure 4 .

[0088] First, the offset range is calculated by the following formula:

[0089] ;

[0090] wherein, represents an empirical coefficient, represents the laser spot diameter.

[0091] Secondly, the random offset of the laser point is determined by the following formula:

[0092] ;

[0093] wherein, represents the random offset of the laser point, represents the laser point offset range, t is the grading coefficient, q represents an integer, and 1≤q≤10.

[0094] Then, according to the results calculated above, the points after the first layer Monte Carlo optimization in step S4 are randomly disturbed in the x direction, y direction or x, y double direction with the spot center as the center within the range of ±0.8 meters or its multiples for 100 times, and all offsets are limited within the radius of 8.75 meters of the spot range, and the offset diagram is as shown in Figure 5 . After each disturbance, the coordinates after the offset are brought into the stereo image pair forward intersection model to generate the corresponding digital elevation model DEM, and the point position matched with the ICESAT-2 satellite laser point and the best point position after the offset correction of the embodiment are compared as Figure 6The (a) figure is a front view image of the high-resolution seven stereoscopic pair, and the (b) figure is a rear view image of the high-resolution seven stereoscopic pair.

[0095] S6, under the offset parameter condition determined in the previous step, a multi-direction perturbation experiment is performed with the center of the laser point spot circle, a plurality of sets of random offset sample coordinate sets are generated, each set of samples corresponds to a spatial geometric hypothesis; each offset sample coordinate is represented as:

[0096] ;

[0097] wherein, represents the laser point offset sample coordinate, represents the original laser point coordinate; control the offset direction, respectively corresponding to the X-axis and Y-axis directions; and represent the random perturbation, respectively obeying the normal distribution:

[0098] ;

[0099] wherein, U(a,b) represents a uniformly distributed random number in the interval [a,b], and respectively represent the maximum perturbation in the X and Y directions.

[0100] S7, for the 100 DEMs generated in step S5, the accuracy is evaluated by using the true elevation values of the independent check points, the average absolute error of each DEM is calculated according to the elevation values of the independent check points on the DEM, and the one with the highest accuracy is selected as the final modeling result, as shown in Figure 3 (b) shows that the average absolute error of each DEM is calculated as follows:

[0101] ;

[0102] wherein, MAE represents the average absolute error of the laser point, Z represents the elevation value of the laser point on the digital elevation model, z represents the true value of the laser point on the digital elevation model, and n represents the number of check points.

[0103] This embodiment constructs three kinds of digital elevation models (DEM) respectively:

[0104] (1) the initial DEM generated based on the static threshold screening control point (i.e. DEM before Monte Carlo optimization);

[0105] (2) DEM generated based on the first layer Monte Carlo random experiment optimization control point set;

[0106] (3) the final DEM based on the second layer Monte Carlo offset correction optimization.

[0107] The DEM precision evaluation method generated for steps S2 and S4 is evaluated by using the mean absolute error. The results of the comparison and analysis of the precision of the three DEMs are shown in Table 2. The DEM precision established by the control point set before the Monte Carlo optimization is 6.38 meters. After the first layer of Monte Carlo optimization, the precision is significantly improved to 4.165 meters. After the second layer of Monte Carlo optimization, the precision is improved to 3.043 meters. That is, the method of this embodiment has a higher improvement.

[0108] Table 2 Quality evaluation table of terrain modeling results

[0109]

[0110] The method of this embodiment first collects satellite laser altimetry data in the study area, filters the original data, and obtains stereo image pairs in the study area to preliminarily arrange control points and perform aerial triangulation encryption. Secondly, a first layer Monte Carlo sampling system is constructed, and high-confidence sample points are selected. Thirdly, a second layer Monte Carlo offset experiment is established, and a multi-direction perturbation experiment is performed to generate multiple groups of random offset samples by analyzing the geometric response relationship between the laser points and the stereo images. Finally, the optimal offset parameters are determined through error response analysis, and the final high-precision digital elevation model is output. This method fully utilizes the laser point reference, combines Monte Carlo optimization and spatial offset sensitivity analysis, realizes the high-precision and robust construction of DEM, and significantly improves the reliability and precision of the digital elevation model.

[0111] The Monte Carlo method is a numerical statistical optimization technique based on random sampling, which has significant advantages in error analysis and uncertainty evaluation. The method of this embodiment can effectively reduce the influence of control point errors on aerial triangulation adjustment solution through a double-layer random optimization mechanism, optimize the systematic correction mechanism of laser point spatial offset, and significantly improve the precision and stability of DEM.

[0112] In addition to the above embodiments, the present application can have other implementation manners. Any technical solution formed by equivalent substitution or equivalent transformation falls within the protection scope required by the present application.

Claims

1. A dual-layer Monte Carlo optimized satellite photogrammetry terrain modeling method, characterized by: The method comprises the following steps: S1, collecting satellite laser altimetry data in the research area, screening the original data according to the characteristics of the satellite laser points, and generating a screened laser altimetry point set; S2, obtaining stereo image remote sensing images in the research area, laying out initial control points and aerotriangulation based on the laser altimetry point set obtained in the previous step, establishing a stereo geometric model and generating a preliminary digital elevation model; S3, based on the laser altimetry point set obtained in step S1, a first layer Monte Carlo sampling system is constructed, a fixed proportion of sample points are randomly selected in each experiment for geometric model parameter reestimation, and the remaining samples are used for independent accuracy evaluation, and a preliminary optimization result is obtained through adjustment optimization; a first layer Monte Carlo experiment is designed, half of the laser points selected from step S1 are randomly selected as control points to participate in geometric parameter reestimation and adjustment optimization, and the remaining half of the points are not involved in adjustment as check points, a total of 300 random extractions are performed, the error of each point is calculated after each extraction, and the error calculation formula is as follows: ； Wherein, Xi represents the observation value of the laser point in the x direction, xi represents the true value of the laser point in the x direction; Yi represents the observation value of the laser point in the y direction, yi represents the true value of the laser point in the y direction; Zi represents the observation value of the laser point in the z direction, and zi represents the true value of the laser point in the z direction; In step S3, the adjustment optimization method adopted is the minimum projection residual optimization, which comprises the following steps: S3.1, the geometric relationship between the ground point and the image coordinate is established by using the collinear equation, and the formula is as follows: ； wherein denotes a ground point; denotes an image coordinate; is an exterior orientation element of the satellite sensor, wherein denotes a spatial position of the satellite sensor in the ground coordinate system, denotes an attitude angle of the satellite sensor, i.e. the rotation angle around the X, Y, Z axes, respectively; S3.2, for each laser point i, the projection residual is defined, and the formula is as follows: ； wherein, and denotes the projection residual of each laser point i, and denotes the real coordinates of the laser points on the image, and are the theoretical coordinates of the laser points on the image; S3.3, taking the sum of squares of residuals as the optimization objective function, the formula is as follows: ； wherein, denotes the projection residual sum of squares objective function; S3.4, the iterative least square method is used for parameter solving, and the formula is as follows: ； wherein denotes a correction of the exterior orientation elements, and denotes the parameter estimates of the kth and k+1th iteration, respectively, A denotes a matrix of partial derivatives, and L denotes a vector of constant terms. S4, error variance and stability analysis is performed on the results of multiple experiments in the previous step, a sample point set with a confidence level greater than 0.8 is selected as the reference point of the stereo image, and a first layer Monte Carlo optimized digital elevation model is generated based on this; S5, using the sample point set selected in the previous step, a second layer Monte Carlo offset experiment is established, the offset range and step are determined by analyzing the geometric response relationship between the laser point and the stereo image, and are used for spatial position error sensitivity evaluation; S6, under the offset parameter conditions determined in the previous step, a multi-direction perturbation experiment is performed with the laser point spot center as the center, a plurality of sets of random offset sample coordinate sets are generated, and each set of samples corresponds to a spatial geometric hypothesis; S7, input each group of offset samples generated in the previous step into the stereo image geometric model, perform forward intersection and parameter reestimation respectively, form a plurality of digital elevation model samples, determine the optimal offset parameter through error response analysis, and output the final digital elevation model.

2. The method of claim 1, wherein: In the step S1, the ICESat-2 satellite ATL08 laser altimetry data of the study area is acquired, first, the original photon point cloud is converted into a standardized format, and abnormal points affected by atmospheric clouds, terrain shielding or photon noise are removed according to laser point attribute information, the laser point attribute information includes signal confidence marker, ground object type identifier, reflection energy distribution, echo shape and detection channel; then the first 30 ground laser points with the highest quality are selected from the removed laser points for subsequent experiments.

3. The method of claim 1, wherein: In the step S2, the high-resolution seven stereo image pair remote sensing image of the study area is acquired, first, the image is preprocessed, and the radiation calibration and geometric correction method is used to eliminate the attitude deviation and light influence; combined with the selected laser points in the step S1, the image is punctured and aerial triangulation is performed, half of the number of laser points are used as control points for aerial triangulation encryption, and the remaining half of the number of points are used as independent check points; the geometric constraint relationship between the ground points and the image points is established by using the collinearity equation, and the initial solution of the exterior orientation elements is obtained by the minimum projection residual method, and a preliminary digital elevation model is generated to provide a benchmark for subsequent optimization.

4. The dual-layer Monte Carlo optimized satellite photogrammetry terrain modeling method of claim 1, wherein: In the step S4, the laser point confidence calculation formula is as follows: ； ； ； wherein, represents the average error of the laser spot, represents the error of the laser spot after the jth Monte Carlo trial, m represents the number of Monte Carlo trials, represents the standard deviation of the point, represents the confidence index of the point.

5. The dual-layer Monte Carlo optimized satellite photogrammetry terrain modeling method of claim 1, wherein: In the step S5, the laser point offset range is determined by the following formula: ； wherein represents an empirical coefficient, represents the laser spot diameter.

6. A dual-layer Monte Carlo optimized satellite photogrammetry terrain modeling method according to claim 5, characterized in that: In the step S5, the offset amount of the laser point is determined according to the set maximum offset range M, and the calculation method is: ； wherein, represents a random offset of the laser spot, represents a range of the laser spot offset, t is a grading coefficient, q represents an integer, and 1 ≤ q ≤ 10.

7. The dual-layer Monte Carlo optimized satellite photogrammetry terrain modeling method of claim 1, wherein: In the step S6, each offset sample coordinate is represented as: ； wherein, represents the laser point offset sample coordinate, represents the original laser point coordinate; control the offset direction, corresponding to the X-axis and Y-axis directions, respectively; and represents a random disturbance, respectively subject to a normal distribution: ； where U(a,b) denotes a uniformly distributed random number in the interval [a,b], and denote the maximum perturbation in the X and Y directions, respectively.

8. The dual-layer Monte Carlo optimized satellite photogrammetry terrain modeling method of claim 1, wherein: In the step S7, the error response analysis is performed by using the mean absolute error analysis method, and the mean absolute error calculation formula is: ； Wherein, MAE represents the average absolute error of the laser point, Z represents the elevation value of the laser point on the digital elevation model, z represents the true value of the laser point on the digital elevation model, and n represents the number of check points.

Images