Photogrammetry point cloud adjustment processing method and system based on satellite photon point cloud

Abstract

The invention relates to the technical field of photogrammetry, and discloses a photogrammetry point cloud adjustment processing method and system based on satellite photon point clouds, and the method comprises the steps: obtaining photon point cloud data and photogrammetry point cloud data; dividing the photon point cloud data into grids of a two-dimensional plane, and screening out a first control point; calculating the topographic relief similarity between the first control point and the photogrammetry point cloud data, and calculating the offset between the first control point and the photogrammetry point cloud data by using a simulated annealing algorithm; performing registration processing on the first control points according to the offset, and screening out second control points meeting a preset condition from the registered first control points; obtaining original RPC parameters of the photogrammetry point cloud data, and carrying out iterative adjustment optimization on the original RPC parameters to obtain RPC parameters meeting precision conditions; and according to the RPC parameter meeting the precision condition, calculating to obtain the adjusted photogrammetric point cloud. The method has the advantage of high universality, and the precision of photogrammetric point cloud is improved.

Description

Technical Field

[0001] The present invention relates to the technical field of photogrammetry, and more specifically, to a photogrammetry point cloud adjustment processing method and system based on satellite photon point cloud. Background Art

[0002] Improving the positioning accuracy of satellite photogrammetry in areas lacking or even without control points is of great significance in a variety of applications such as natural resource development and protection in mountainous areas, geological disaster monitoring and assessment, and infrastructure construction. However, due to the complex terrain of mountainous areas, traditional satellite photogrammetry methods often have difficulty in achieving ideal positioning accuracy in these areas, so it is urgent to develop new technical methods to solve this problem.

[0003] At present, some scholars have tried to jointly process satellite image data and photon radar data to establish the correspondence between laser points and images. For example, some studies have selected laser points in flat areas, and used the original or free network adjusted image positioning model to calculate the image projection points of the laser points as image corresponding points. Other scholars have filtered the photon point cloud data and selected points with satisfactory accuracy as control points for adjustment. However, these methods have some problems: the former inevitably has errors in the positioning model, resulting in inaccurate correspondence between laser points and images, especially in rugged mountainous and forested areas, where the deviation is greater. The latter is not universal for rugged mountainous and forested areas because the registration method is too single, and the accuracy of the registration results is not high enough. Summary of the invention

[0004] In order to overcome the defects of low universality and low precision in the prior art, the present invention proposes the following technical solutions: In a first aspect, the present invention proposes a photogrammetric point cloud adjustment processing method based on satellite photon point cloud, comprising: Acquire photon point cloud data and photogrammetric point cloud data.

[0005] The photon point cloud data is divided into a two-dimensional grid, and photon points meeting preset conditions are screened out from the grid as first control points.

[0006] Calculate the terrain relief similarity between the first control point and the photogrammetric point cloud data.

[0007] According to the terrain relief similarity, a simulated annealing algorithm is used to calculate the offset between the first control point and the photogrammetric point cloud data.

[0008] According to the offset, the first control points are registered, and second control points meeting preset conditions are selected from the registered first control points.

[0009] Get the raw RPC parameters of photogrammetric point cloud data.

[0010] According to the second control point, the original RPC parameters are iteratively adjusted and optimized to obtain RPC parameters that meet the accuracy conditions.

[0011] According to the RPC parameters that meet the accuracy conditions, the adjusted photogrammetric point cloud is calculated.

[0012] As a preferred technical solution, the photon point cloud data is divided into a two-dimensional grid, and photon points meeting preset conditions are screened out from the grid as control points, including: The photon point cloud data is divided into a two-dimensional grid with the y-axis as the horizontal coordinate and the z-axis as the vertical coordinate.

[0013] Count the number of photon points in each grid , the number of photon points Reaching the threshold The photon points in the grid are retained and the corresponding matrix is ​​constructed , where the retained grid is in the matrix The corresponding position is 1, and the rest are 0. The specific expression is as follows:

[0014]

[0015] in, i Indicates the vertical position of the grid and is used to traverse or reference different rows in the matrix. j Indicates the horizontal position of the grid and is used to traverse or reference different columns in the matrix.

[0016] According to the matrix Calculate the continuity score for each raster , whose expression is as follows:

[0017] Group the rasters in each column and keep the rasters with the largest continuity score and a continuity score greater than the threshold in each raster group. Photon points in a grid As a control point, its expression is as follows:

[0018] .

[0019] As a preferred technical solution, calculate the elevation difference between the control point and the photogrammetric point cloud data As terrain relief similarity, it includes: The photogrammetric points and control points The elevation is set to 0, the horizontal and vertical coordinates remain unchanged, and saved to the 2D photogrammetry point cloud carrier and 2D photon point cloud carrier The expression is as follows: ;

[0020] For two-dimensional photon point cloud carrier For each control point in , find the 2D photogrammetric point cloud vector The photogrammetric point with the smallest Euclidean distance is taken as the corresponding point cloud , whose expression is as follows:

[0021]

[0022] in, Represents the Euclidean distance between the calculated control points and the photogrammetric points.

[0023] The absolute value of the elevation difference between the control point and its corresponding photogrammetric point is calculated as the elevation difference, and its expression is as follows:

[0024] in, and are the elevations of the control point and its corresponding photogrammetric point, respectively.

[0025] As a preferred technical solution, according to the terrain undulation similarity, a simulated annealing algorithm is used to calculate the offset between the control point and the photogrammetric point cloud data, including: Change the spatial position of a control point within a preset search area.

[0026] Get the current temperature of the simulated annealing algorithm , whose expression is as follows:

[0027] in, n represents the number of iterations, is the cooling rate.

[0028] According to the current temperature, the control point is calculated x and y The offset in the direction is expressed as follows:

[0029]

[0030] in, for A random number.

[0031] Update the coordinates of the control point according to the offset of the control point in the x and y directions.

[0032] Calculate the elevation difference change of the control point based on the updated coordinates , whose expression is as follows:

[0033] in, Represents the elevation difference between the control point and the photogrammetric point. n Indicates the number of iterations.

[0034] Determine whether the elevation difference change meets the preset conditions. If not, iterate the above steps. If yes, output the current accumulated offsets in the x and y directions as the offsets of the optimal position of the control point on the horizontal plane according to the following formula:

[0035] .

[0036] As a preferred technical solution, after obtaining the offset of the optimal position of the control point on the horizontal plane, the method further includes: The calculation control point and the corresponding photogrammetric point are x , y and z Offset in direction , whose expression is as follows:

[0037]

[0038]

[0039] in, represents the total number of photogrammetric point clouds, No. i control points x coordinate, For the i The control points correspond to the photogrammetric points x coordinate, No. i control points y coordinate, For the i The control points correspond to the photogrammetric pointsy coordinate, No. i control points z coordinate, For the i The control points correspond to the photogrammetric points z coordinate.

[0040] As a preferred technical solution, according to the second control point, the original RPC parameters are iteratively adjusted and optimized to obtain RPC parameters that meet the accuracy conditions, including: Step a: Build the RFM model of the original RPC parameters.

[0041] Step b: Perform affine transformation compensation adjustment on the RFM model.

[0042] Step c: According to the second control point, the RFM model subjected to the adjustment process of affine transformation compensation is updated with the terrain-independent algorithm in turn, and the adjustment process of Fourier compensation is performed.

[0043] Step d: Determine whether the RPC parameters of the RFM model meet the preset accuracy conditions. If so, output the RPC parameters. If not, jump to step c.

[0044] As a preferred technical solution, according to the second control point, the RFM model adjusted by affine transformation compensation is updated with a terrain-independent algorithm using RPC parameters, including: According to the second control points, grid points are extracted at equal intervals in the image plane, and the surface is divided into several layers according to the elevation range.

[0045] According to the image space coordinates of each grid point, the object space coordinates of each grid point on each surface layer are calculated using the RFM model.

[0046] A virtual control grid is constructed based on the object coordinates of each grid point on each surface layer.

[0047] Based on the virtual control grid, the RPC parameters of the RFM model are updated using the least squares method.

[0048] As a preferred technical solution, after obtaining the adjusted photogrammetric point cloud coordinates, the method further includes: Calculate the average Euclidean distance between the preset verification point cloud and the adjusted photogrammetric point cloud. The expression is as follows:

[0049]

[0050]

[0051] in, is the average Euclidean distance in the x direction between the preset verification point cloud and the adjusted photogrammetric point cloud, is the average Euclidean distance in the y direction between the preset verification point cloud and the adjusted photogrammetric point cloud, is the average Euclidean distance in the z direction between the preset verification point cloud and the adjusted photogrammetric point cloud, The preset verification point cloud i Points x coordinate, The preset verification point cloud i Points y coordinate, The preset verification point cloud i Points z coordinate.

[0052] The accuracy of the photogrammetric point cloud after adjustment was evaluated based on the average Euclidean distance.

[0053] As a preferred technical solution, judging whether the elevation difference change meets a preset condition includes: Determine the current elevation difference change Is it less than the elevation difference change of the previous iteration? And the current temperature Is it less than the threshold, or , if so, then output the current accumulated offsets in the x and y directions as the offsets of the optimal position of the control point on the horizontal plane.

[0054] In a second aspect, the present invention further proposes a photogrammetric calibration system based on satellite photon point cloud, which is applied to the photogrammetric point cloud adjustment processing method based on satellite photon point cloud as described in any scheme of the first aspect, comprising: The first acquisition module acquires photon point cloud data and photogrammetry point cloud data.

[0055] The first screening module divides the photon point cloud data into a two-dimensional grid, and screens out photon points that meet preset conditions from the grid as first control points.

[0056] The first calculation module calculates the terrain relief similarity between the first control point and the photogrammetric point cloud data.

[0057] The second calculation module calculates the offset between the first control point and the photogrammetric point cloud data using a simulated annealing algorithm according to the terrain undulation similarity.

[0058] The second screening module performs registration processing on the first control points according to the offset, and screens out second control points that meet preset conditions from the registered first control points.

[0059] The second acquisition module acquires the original RPC parameters of the photogrammetric point cloud data.

[0060] The optimization module performs iterative adjustment optimization on the original RPC parameters according to the second control point to obtain RPC parameters that meet the accuracy condition.

[0061] The third calculation module calculates the adjusted photogrammetric point cloud according to the RPC parameters that meet the accuracy conditions.

[0062] The beneficial effects of the present invention include at least: The present invention effectively utilizes the high-precision characteristics of ICESat-2 satellite photon point cloud data by dividing the photon point cloud data into a two-dimensional plane grid and screening photon points that meet preset conditions as control points, providing reliable reference points for subsequent adjustment processing. The similarity of terrain undulation is used as an indicator and a simulated annealing algorithm is adopted to calculate the offset between the control point and the photogrammetric point to be processed. Compared with the traditional grid search method, the registration accuracy and universality are greatly improved, and it is particularly suitable for mountainous areas with complex terrain. The original RPC parameters of the photogrammetric point cloud data are obtained, and these parameters are iteratively adjusted and optimized using the screened high-quality second control points until the accuracy requirements are met. Finally, the adjusted photogrammetric point cloud is calculated using the optimized RPC parameters, which not only improves the reliability and accuracy of the control points, but also ensures the high accuracy of the RPC parameters, thereby significantly improving the overall accuracy of the final photogrammetric point cloud. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 This is a flow chart of the photogrammetry point cloud adjustment processing method based on satellite photon point cloud provided in Example 1.

[0064] Figure 2 This is a flow chart of the simulated annealing algorithm provided in Example 2.

[0065] Figure 3 The flowchart of iterative adjustment optimization of RPC parameters provided in Example 2.

[0066] Figure 4 This is a schematic diagram of a virtual control grid for calculating RPC parameters provided in Example 2.

[0067] Figure 5 Schematic diagram of the positional relationship between the photon point cloud and the photogrammetry point cloud of the actual measurement area provided in Example 3.

[0068] Figure 6A schematic diagram of selecting the number of photon point cloud control points provided in Example 3.

[0069] Figure 7 (a) is the elevation distribution diagram of the photon point cloud provided in Example 3, Figure 7 (b) is the elevation distribution map of the photogrammetric point cloud provided in Example 3.

[0070] Figure 8 This is the photon point cloud and photogrammetry point cloud registration data diagram provided in Example 3.

[0071] Fig. 9 This is a schematic diagram of the accuracy assessment results provided in Example 3.

[0072] Fig.10 This is a comparison view of the a window calibration results provided in Example 3.

[0073] Fig.11 This is the comparison view of the original data of window b provided in Example 3.

[0074] Fig.12 This is a comparison view of the mountain local calibration results of window a provided in Example 3.

[0075] Fig.13 This is a comparison view of the local original data of the mountain in window b provided in Example 3.

[0076] Fig.14 This is an architectural diagram of the photogrammetry calibration system based on satellite photon point cloud provided in Example 4. DETAILED DESCRIPTION

[0077] The following will describe the embodiments of the present invention with reference to the accompanying drawings and preferred technical solutions. Those skilled in the art can easily understand other advantages and effects of the present invention from the contents disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments, and the details in this specification can also be modified or changed in various ways based on different viewpoints and applications without departing from the spirit of the present invention. It should be understood that the preferred technical solutions are only for illustrating the present invention, not for limiting the scope of protection of the present invention.

[0078] It should be noted that the illustrations provided in the following embodiments are only schematic illustrations of the basic concept of the present invention, and thus the drawings only show components related to the present invention rather than being drawn according to the number, shape and size of components in actual implementation. In actual implementation, the type, quantity and proportion of each component may be changed arbitrarily, and the component layout may also be more complicated.

[0079] In the following description, numerous details are discussed to provide a more thorough explanation of the embodiments of the present invention. However, it is obvious to those skilled in the art that the embodiments of the present invention can be implemented without these specific details. In other embodiments, well-known structures and devices are shown in the form of block diagrams rather than in detail to avoid making the embodiments of the present invention difficult to understand.

[0080] Example 1 This embodiment proposes a photogrammetry point cloud adjustment processing method based on satellite photon point cloud, such as Figure 1 As shown, Figure 1 A schematic flow chart of a photogrammetric point cloud adjustment processing method based on satellite photon point cloud provided in this embodiment, the method comprises the following steps: S1: Acquire photon point cloud data and photogrammetry point cloud data.

[0081] S2: Divide the photon point cloud data into a two-dimensional grid, and select photon points that meet preset conditions from the grid as first control points.

[0082] S3: Calculate the terrain relief similarity between the first control point and the photogrammetric point cloud data.

[0083] S4: Calculate the offset between the first control point and the photogrammetric point cloud data using a simulated annealing algorithm according to the terrain undulation similarity.

[0084] S5: performing registration processing on the first control points according to the offset, and selecting second control points that meet preset conditions from the registered first control points.

[0085] S6: Get the original RPC parameters of the photogrammetric point cloud data.

[0086] S7: According to the second control point, the original RPC parameters are iteratively adjusted and optimized to obtain RPC parameters that meet the accuracy conditions.

[0087] S8: Calculate the adjusted photogrammetric point cloud based on the RPC parameters that meet the accuracy requirements.

[0088] It can be understood that the present invention effectively utilizes the high-precision characteristics of ICESat-2 satellite photon point cloud data by dividing the photon point cloud data into a two-dimensional plane grid and selecting photon points that meet the preset conditions as control points, providing reliable reference points for subsequent adjustment processing. Using the similarity of terrain undulation as an indicator and adopting the simulated annealing algorithm to calculate the offset between the control point and the photogrammetric point to be processed, compared with the traditional grid search method, the registration accuracy and universality are greatly improved, and it is particularly suitable for mountainous areas with complex terrain. By obtaining the original RPC parameters of the photogrammetric point cloud data, and using the selected high-quality second control points to iteratively adjust and optimize these parameters until the accuracy requirements are met. Finally, the adjusted photogrammetric point cloud is calculated using the optimized RPC parameters, which not only improves the reliability and accuracy of the control points, but also ensures the high accuracy of the RPC parameters, thereby significantly improving the overall accuracy of the final photogrammetric point cloud.

[0089] The present invention is applicable to various terrain conditions, especially in complex terrain areas with a lack of control points or no control points, and can be widely used in fields such as natural resource surveys in mountainous areas, geological disaster monitoring, and infrastructure construction. The entire processing flow has a high degree of automation, which reduces manual intervention, improves processing efficiency, and reduces the impact of human errors, ensuring the objectivity and reliability of the results. In summary, the present invention not only significantly improves the accuracy of satellite photogrammetry in areas without control points, but also has a high degree of universality, providing an effective technical means for solving photogrammetry problems in complex terrains such as mountainous areas.

[0090] Example 2 This embodiment makes improvements on the photogrammetry point cloud adjustment processing method based on satellite photon point cloud proposed in Embodiment 1.

[0091] In this embodiment, the photon point cloud data is divided into a two-dimensional grid, and photon points meeting preset conditions are screened out from the grid as control points, including: With the y-axis as the horizontal coordinate and the z-axis as the vertical coordinate, the two-dimensional plane of the photon point cloud data is divided into Gao Wei The grid is created and the grid side length is set according to the minimum and maximum values ​​of the y coordinates and the minimum and maximum values ​​of the z coordinates of the photon point cloud data.

[0092] Count the number of photon points in each grid , the number of photon points Reaching the threshold The photon points in the grid are retained and the corresponding matrix is ​​constructed , where the retained grid is in the matrix The corresponding position is 1, and the rest are 0. The specific expression is as follows:

[0093]

[0094] in, i Indicates the vertical position of the grid and is used to traverse or reference different rows in the matrix. j Indicates the horizontal position of the grid and is used to traverse or reference different columns in the matrix.

[0095] Through the matrix The distribution of the retained grids can be clearly seen, and the continuity score of each grid is given by judging whether there are retained grids in the eight directions around the retained grid. , one point is added for each retained grid around it, and its expression is as follows:

[0096] Group the rasters in each column and keep the rasters with the largest continuity score and a continuity score greater than the threshold in each raster group. Photon points in a grid As a control point, its expression is as follows:

[0097] .

[0098] In this embodiment, the sum of the absolute values ​​of the elevation differences of the nearest points of the plane coordinates is used. Evaluate the global terrain similarity. The more similar the terrain is, the greater the sum of the absolute values ​​of the elevation differences. The smaller.

[0099] The elevation difference between the control point and the photogrammetric point The calculation steps include: The photogrammetric points and control points The elevation is set to 0, the horizontal and vertical coordinates remain unchanged, and saved to the 2D photogrammetry point cloud carrier and 2D photon point cloud carrier The expression is as follows: ;

[0100] For two-dimensional photon point cloud carrier For each control point in , find the 2D photogrammetric point cloud vector The photogrammetric point with the smallest Euclidean distance is taken as the corresponding point cloud , whose expression is as follows:

[0101]

[0102] in, Represents the Euclidean distance between the calculated control points and the photogrammetric points.

[0103] The absolute value of the elevation difference between the control point and its corresponding photogrammetric point is calculated as the elevation difference, and its expression is as follows:

[0104] in, and are the elevations of the control point and its corresponding photogrammetric point, respectively.

[0105] After traversing the photon points, add up the absolute values ​​of all elevation differences to get the terrain similarity evaluation :

[0106] in, Represents the total number of photon point clouds.

[0107] In this embodiment, the simulated annealing algorithm is used to calculate the offset between the control point and the photogrammetric point to be processed. The simulated annealing algorithm makes up for the shortcomings of the grid search method, the ICP registration method, the excessive reliance on normal vectors, and the long calculation time. The basic logic of the simulated annealing algorithm is to change the spatial position of the control point to obtain different evaluations of the similarity of terrain relief, and finally take the position with the most similar terrain relief, that is, the smallest, to record the change of the control point. and ,like Figure 2 As shown, Figure 2 Flow chart of the simulated annealing algorithm provided for this embodiment: According to the terrain undulation similarity, the offset between the control point and the photogrammetric point cloud data is calculated using a simulated annealing algorithm, including: Change the spatial position of a control point within a preset search area.

[0108] Get the current temperature of the simulated annealing algorithm , whose expression is as follows:

[0109] in, n represents the number of iterations, is the cooling rate.

[0110] According to the current temperature, the control point is calculated x and y The offset in the direction is expressed as follows:

[0111]

[0112] in, for A random number.

[0113] Update the coordinates of the control point according to the offset of the control point in the x and y directions.

[0114] Calculate the elevation difference change of the control point based on the updated coordinates , whose expression is as follows:

[0115] in, Represents the elevation difference between the control point and the photogrammetric point. n Indicates the number of iterations.

[0116] Determine the current elevation difference change Is it less than the elevation difference change of the previous iteration? And the current temperature Is it less than the threshold, or , if so, then according to the following formula, output the current cumulative offset in the x and y directions as the offset of the optimal position of the control point on the horizontal plane:

[0117] .

[0118] In this embodiment, after obtaining the offset of the optimal position of the control point on the horizontal plane, the method further includes: The calculation control point and the corresponding photogrammetric point are x , y and z Offset in direction , whose expression is as follows:

[0119]

[0120]

[0121] in, represents the total number of photogrammetric point clouds, No. i control points x coordinate, For the i The control points correspond to the photogrammetric pointsx coordinate, No. i control points y coordinate, For the i The control points correspond to the photogrammetric points y coordinate, No. i control points z coordinate, For the i The control points correspond to the photogrammetric points z coordinate.

[0122] In this embodiment, according to the second control point, the original RPC parameters are iteratively adjusted and optimized to obtain RPC parameters that meet the accuracy conditions, such as Figure 3 As shown, including: Step a: Build the RFM model of the original RPC parameters.

[0123] As an exemplary explanation, the RFM (rational function model) model is further extended on the basis of the traditional direct linear transformation model (Direct Linear Transformation, DLT), and the transformation from the geodetic coordinates of the object space to the image point coordinates of the image space is described in the form of the ratio of two cubic polynomials, that is, the inverse solution form of the RFM model, as shown below:

[0124] in, is the normalized image point coordinate in the image plane coordinate system, is the normalized geodetic coordinate of the object space coordinate system. The purpose of the normalization parameter is to ensure that the values ​​of the image space coordinates and the object space coordinates are within the range of -1 to 1, avoiding the calculation error caused by the difference in the order of magnitude of the coordinates. Specifically:

[0125]

[0126]

[0127]

[0128]

[0129] in, Represents the image point coordinates in the image plane coordinate system. The object coordinates of ground points can be represented by geocentric rectangular coordinates, geodetic coordinates or map projection coordinates. is the normalization parameter of the image coordinate. is the normalized parameter of the object coordinate. , , , All represent polynomials. Here we introduce the cubic polynomial. For example, the form is as follows:

[0130] in, are the polynomial coefficients, i.e., the RPC parameters, , The constant term of is 1.

[0131] Step b: Perform affine transformation compensation adjustment on the RFM model.

[0132] Step c: According to the second control point, the RFM model subjected to the adjustment process of affine transformation compensation is updated with the terrain-independent algorithm in turn, and the adjustment process of Fourier compensation is performed.

[0133] As an exemplary explanation, in order to correct the positioning error of the RFM model, this embodiment adds a system error compensation term for the image coordinates without changing the original RPC parameters, and uses the compensation term coefficient as an adjustment parameter. The observation equation is shown in the following formula:

[0134] in is the calculated value of the image point coordinates obtained using the original RPC, for The error compensation term can generally be expressed as the image point observation value The more the error compensation term used in the adjustment conforms to the actual error distribution, the higher the adjustment accuracy.

[0135] This embodiment performs adjustment through affine transformation compensation term and Fourier compensation term, specifically: The first-order polynomial, i.e., the affine transformation formula, can effectively eliminate the image-side system error. The form of the affine transformation compensation term of this embodiment is shown in the following formula:

[0136] Theoretically, the Fourier series can approach any function that meets the continuity condition in a two-dimensional interval. Studies have shown that the two-dimensional Fourier polynomial has a good compensation effect on complex image distortion. The form of the Fourier compensation term of this embodiment is shown in the following formula:

[0137] in, and , and are the Fourier polynomial coefficients. and are the cosine component and the sine component respectively, as shown in the following formula:

[0138]

[0139] in, is the image width, is the image height, is the order. When , the Fourier polynomial is as follows:

[0140] Furthermore, if the error model The error compensation term in Included compensation coefficients, taking into account the unknown object coordinates of the same-name connection points extracted from the stereo image, The error equation of the linearized RFM adjustment is as follows:

[0141] It can be written into the error equation matrix form as:

[0142] in, is the residual vector of the image point coordinate observation. is the coefficient matrix, that is, the matrix of partial derivatives with respect to the unknowns. is the compensation factor The increment vector of . is the object space coordinate of the connection point The increment vector of . is a constant term, where is the observed value of the image point coordinates, To use the approximate value of the unknown number into the formula The calculated image point coordinates. is the weight matrix, which is set to the unit matrix when the accuracy of each observation point is equal For ground control points, the coordinate error equation of the control points can be listed. Specifically

[0143] Among them, , the weight of the control point coordinate observation It depends on the ratio of the image point coordinate variance to the corresponding coordinate variance, that is:

[0144] For the photon elevation point, the coordinate error equation is the same as the general control point coordinate error equation, but its weight should be determined according to the photon point coordinate accuracy, as shown in the following formula:

[0145] in, , the error in the photon point coordinates in the data of this embodiment is about 3 times the error in the control point coordinates, so the weight of the photon point coordinate observation value can be set to 1 / 9 of the weight of the control point coordinate observation value, that is:

[0146] Depend on The normal equation is established according to the principle of least square adjustment as follows to solve the compensation term coefficient and the object coordinate increment of the connection point:

[0147] The compensation coefficient correction number can be solved at the same time through the above formula And the unknown point coordinate correction number .

[0148] Step d: Determine whether the RPC parameters of the RFM model meet the preset accuracy conditions. If so, output the RPC parameters. If not, jump to step c.

[0149] In this embodiment, according to the second control point, the RPC parameters of the RFM model adjusted by affine transformation compensation are updated by using a terrain-independent algorithm, including: According to the second control points, grid points are extracted at equal intervals in the image plane, and the surface is divided into several layers according to the elevation range.

[0150] According to the image space coordinates of each grid point, the object space coordinates of each grid point on each surface layer are calculated using the RFM model.

[0151] A virtual control grid is constructed based on the object coordinates of each grid point on each surface layer.

[0152] Specifically, according to the left information of the second control point, the image plane is uniformly spaced to extract Grid points , and divide the surface elevation range into Layer, The elevation of the layer is For each grid point, the image coordinates , calculate the three-dimensional coordinates of the object point corresponding to each layer of the surface through the geometric relationship between the object and the image , thus constructing a The virtual control grid, such as Figure 4 As shown, Figure 4 Schematic diagram of the virtual control grid for calculating RPC parameters.

[0153] The specific steps of using the RFM model to calculate the object space coordinates of each grid point on each surface layer include: 1) Assign initial values ​​to the object plane coordinates of the point to be sought , .

[0154] 2) Settings Initial step size of coordinates .

[0155] 3) Calculate the object coordinate correction coefficient ,in The object coordinates are Substitution The image coordinates obtained are: = +d , = +d .

[0156] 4) Calculate the object coordinate correction number under the current step size: , get the corrected object coordinates = +∆ , = +∆ .

[0157] 5) Calculate the image coordinate deviation after object coordinate correction ,in Substitution The obtained image coordinates.

[0158] 6) Inspection Is it less than the tolerance (0.01 pixel), if so, output the object coordinate result , otherwise replace , and reduce the step size , repeat steps 3) to 6) until the accuracy requirements are met.

[0159] Based on the virtual control grid, the RPC parameters of the RFM model are updated using the least squares method.

[0160] It can be understood that this embodiment expands one adjustment into multiple adjustments and proposes an RFM iterative adjustment method with variable RPC. After performing an image square adjustment on RFM, a virtual control grid is constructed based on the adjustment results and the RPC parameters are recalculated. Since the new RPC parameters have eliminated the influence of some systematic errors, the positioning results are closer to the image point observation values. The next image square adjustment can be performed on this basis to eliminate the residual systematic errors. After several iterative calculations, the complex image distortion is gradually eliminated.

[0161] Example 3 This embodiment verifies the reliability of RFM iterative adjustment through simulation experiments based on the method proposed in Example 2, and after verifying that the adjustment method is effective and determining the optimal number of iterations, performs measurement adjustment on the actual survey area and evaluates the accuracy of the method of the present invention.

[0162] The simulation test of this embodiment uses RSM to generate 1:10000 scale dual-line array stereo camera simulated observation data, including RPC parameter initial values, object coordinates of target points and front and back image coordinates. The effectiveness of the adjustment method used is tested by comparing the target positioning accuracy before and after RFM adjustment. The specific simulation parameters are shown in Table 1. Among them, the simulated distortion value that conforms to the quadratic polynomial distribution is added to the image coordinates of the target points. , Corner elements , Main Point , main distance The corresponding errors ΔL, ΔA, ΔP, and Δf are introduced into the “true values” of the parameters, and the initial values ​​of the RPC parameters are calculated using the terrain-independent algorithm based on RSM.

[0163] Table 1 Optical satellite stereo image simulation parameters

[0164] When ΔL=100m, ΔA=100", ΔP=100 pixels, Δf=10mm, the image line edge is distorted The RFM adjustment positioning accuracy of different iteration times is shown in Table 2. They are the image-square positioning errors on the foresight and backsight images, respectively, in pixels. They are the object plane and elevation positioning errors respectively, in meters.

[0165] Table 2 Positioning results of RFM iterative adjustment simulation test

[0166] Since the front and back stereo intersection angles of the dual-line array camera are small, the errors of the internal and external orientation elements produce large object space elevation positioning errors. As shown in Table 2, when the number of iterations is 1, that is, the conventional RFM adjustment based on the affine transformation model is used, the front and back image space positioning accuracy is about 3 pixels. At this time, the error in elevation positioning is 8.5m, which is 5 times the error in plane, far from the positioning accuracy requirement of 1:10000 scale. After the second adjustment, the image space system error was greatly weakened, the image space positioning error was reduced by 94%, and the errors in plane and elevation positioning were reduced by 93%. The third adjustment still has a significant compensation effect on the system error. The error in image space positioning is further reduced by 22% to within 0.15 pixels, and the errors in object space plane and elevation positioning are reduced by 28% and 25%. At this time, the adjustment effect is close to the maximum limit. If the number of iterations is increased, the positioning accuracy of the image space and object space will basically not change. Taking into account the computational performance and cost, the number of RFM adjustment iterations can be set to 3.

[0167] This example uses the southwestern area of ​​Kunming, Yunnan as the actual measurement area for measurement adjustment, and uses the photon point cloud ATL03 and ATL08 data of ICESat-2 / ATLAS in this area on March 6, 2020, with a distance of 30km along the track as experimental data. Figure 5 As shown, Figure 5 The photon point cloud and the photogrammetry point cloud position relationship diagram of the real measurement area in this embodiment, wherein the red part is the photon point cloud, the terrain of this area varies, the highest altitude is 2655.29m, the lowest is 1413.57m, and the area is mainly covered with forest vegetation. The area of ​​the adjusted photogrammetry point cloud is about 48, and the number of point clouds is about 50 million. The photon point cloud ATL03 and ATL08 data of ICESat-2 / ATLAS used in this embodiment are processed and applied by the program after Gaussian coordinate conversion. The satellite image experimental data parameters are shown in Table 3.

[0168] Table 3 Satellite image experimental data parameters

[0169] This embodiment imports the corresponding ICESat-2 / ATLAS photon point cloud data and unadjusted photogrammetric point cloud data as processing files. Because this experiment involves verifying the accuracy and feasibility of the method, the adjusted photogrammetric point cloud will be imported as an evaluation standard, which is not required in actual applications.

[0170] After the file is correctly imported, click Start Debugging and wait for the result output. The first result outputted is, from top to bottom, the serial number of the strip, the total number of photon point clouds of the strip, and the number of control points selected after filtering the photon point clouds of the strip. Figure 6 shown.

[0171] like Figure 7 and Figure 8 As shown, Figure 7 This is a comparison of the elevation distribution of the photon point cloud and the photogrammetry point cloud. Figure 7 (a) is the elevation distribution map of the photon point cloud. Figure 7 (b) is the elevation distribution map of the photogrammetric point cloud. Figure 8 This is the registration data of photon point cloud and photogrammetry point cloud. From the image and data, we can see that the registration effect using simulated annealing algorithm and terrain undulation similarity as evaluation criteria is very good. Points with registration differences greater than 10m only account for 2.8% of the total points, and the maximum registration difference is 42.6m.

[0172] If a verification file is imported, the accuracy compared with the verification file will be output, that is, the offset between the verification file and the corresponding point cloud in the adjustment result, from top to bottom, before and after error correction The offset in three directions, such as Fig. 9 shown.

[0173] After completing the calibration (adjustment) process of the photogrammetric point cloud data, the system will generate the adjusted point cloud data. If the user needs, this data can be saved as a point cloud file. In addition, if the user imports a point cloud file for verification, the system will display two comparison windows: Window a: shows the point cloud in the validation file (in blue) and the photogrammetric point cloud adjusted using this method (in red). The purpose of this window is to visually compare the original validation point cloud and the calibrated point cloud to evaluate the effectiveness of the calibration.

[0174] Fig.10 This is the comparison view of the calibration results in window a, which focuses on the comparison between the calibrated point cloud (red) and the verification point cloud (blue) in a certain local area. Fig.12 This is a comparison view of the local mountain calibration results in window a, which is used to demonstrate the effect of the calibration process in complex terrain.

[0175] b Window: Displays the point cloud in the validation file (in blue) and the unadjusted photogrammetric point cloud (in red). This window is used to show the difference between the uncalibrated point cloud data and the validation point cloud.

[0176] Fig.11 This is the raw data comparison view of window b, showing the comparison between the uncalibrated point cloud (red) and the verification point cloud (blue) in a local area. Fig.13 This is a comparison view of the local original data of the mountain in window b, which further compares the representation of the uncalibrated point cloud data in complex terrain and the differences with the verification point cloud.

[0177] In this embodiment, after obtaining the adjusted photogrammetric point cloud coordinates: Calculate the average Euclidean distance between the preset verification point cloud and the adjusted photogrammetric point cloud. The expression is as follows:

[0178]

[0179]

[0180] in, is the average Euclidean distance in the x direction between the preset verification point cloud and the adjusted photogrammetric point cloud, is the average Euclidean distance in the y direction between the preset verification point cloud and the adjusted photogrammetric point cloud, is the average Euclidean distance in the z direction between the preset verification point cloud and the adjusted photogrammetric point cloud, The preset verification point cloud i Points x coordinate, The preset verification point cloud i Points y coordinate, The preset verification point cloud i Points z coordinate.

[0181] The accuracy of the photogrammetric point cloud after adjustment was evaluated based on the average Euclidean distance.

[0182] Table 4 Comparison results of the accuracy of the present invention and the existing method

[0183] From the results in Table 4, it can be seen that the point cloud files after adjustment by the present invention and the verification point cloud files are basically overlapped in both mountainous areas and forested hilly areas, while the point cloud files that have not been adjusted have poor overlap with the verification point cloud files. Compared with the existing methods, the present invention has doubled the accuracy in the horizontal direction and increased the accuracy tenfold in the vertical direction, which can effectively prove the feasibility, universality and high precision of the present invention.

[0184] Example 4 like Fig.14As shown, this embodiment proposes a photogrammetric calibration system based on satellite photon point cloud, which is applied to the photogrammetric point cloud adjustment processing method based on satellite photon point cloud as described in the above embodiment, including: a first acquisition module 100, a first screening module 200, a first calculation module 300, a second calculation module 400, a second screening module 500, a second acquisition module 600, an optimization module 700 and a third calculation module 800.

[0185] Among them, the first acquisition module 100 acquires photon point cloud data and photogrammetric point cloud data. The first screening module 200 divides the photon point cloud data into a two-dimensional grid, and screens out photon points that meet the preset conditions from the grid as the first control point. The first calculation module 300 calculates the similarity of terrain undulation between the first control point and the photogrammetric point cloud data. The second calculation module 400 calculates the offset between the first control point and the photogrammetric point cloud data using a simulated annealing algorithm based on the similarity of terrain undulation. The second screening module 500 performs registration processing on the first control point based on the offset, and screens out the second control point that meets the preset conditions from the first control point after registration. The second acquisition module 600 acquires the original RPC parameters of the photogrammetric point cloud data. The optimization module 700 performs iterative adjustment optimization on the original RPC parameters based on the second control point to obtain the RPC parameters that meet the accuracy conditions. The third calculation module 800 calculates the adjusted photogrammetric point cloud based on the RPC parameters that meet the accuracy conditions.

[0186] It should be noted that the aforementioned explanation of the embodiment of the photogrammetric point cloud adjustment processing method based on satellite photon point cloud is also applicable to the photogrammetric calibration system based on satellite photon point cloud of this embodiment, and will not be repeated here.

[0187] In the description of this specification, the description with reference to the terms "one embodiment", "some embodiments", "example", "specific example", or "some examples" etc. means that the specific features, structures, materials or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of the present application. In this specification, the schematic representations of the above terms are not necessarily directed to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described may be combined in any one or N embodiments or examples in a suitable manner. In addition, those skilled in the art may combine and combine the different embodiments or examples described in this specification and the features of the different embodiments or examples, without contradiction.

[0188] In addition, the terms "first" and "second" are used for descriptive purposes only and should not be understood as indicating or implying relative importance or implicitly indicating the number of the indicated technical features. Therefore, the features defined as "first" and "second" may explicitly or implicitly include at least one of the features. In the description of this application, the meaning of "N" is at least two, such as two, three, etc., unless otherwise clearly and specifically defined.

[0189] Any process or method description in a flowchart or otherwise described herein may be understood to represent a module, fragment or portion of code comprising one or more executable instructions for implementing the steps of a custom logical function or process, and the scope of the preferred embodiments of the present application includes alternative implementations in which functions may not be performed in the order shown or discussed, including performing functions in a substantially simultaneous manner or in reverse order depending on the functions involved, which should be understood by technicians in the technical field to which the embodiments of the present application belong.

[0190] It should be understood that the various parts of the present application can be implemented by hardware, software, firmware or a combination thereof. In the above-mentioned embodiment, the N steps or methods can be implemented by software or firmware stored in a memory and executed by a suitable instruction execution system. For example, if implemented by hardware, as in another embodiment, it can be implemented by any one of the following technologies known in the art or their combination: a discrete logic circuit having a logic gate circuit for implementing a logic function for a data signal, a dedicated integrated circuit having a suitable combination of logic gate circuits, a programmable gate array, a field programmable gate array, etc.

[0191] A person skilled in the art may understand that all or part of the steps in the above-mentioned embodiment method may be completed by instructing related hardware through a program, and the program may be stored in a computer-readable storage medium, which, when executed, includes one or a combination of the steps of the method embodiment.

[0192] Obviously, the above embodiments of the present invention are merely examples for clearly illustrating the present invention, and are not intended to limit the embodiments of the present invention. For those skilled in the art, other different forms of changes or modifications can be made based on the above description. It is not necessary and impossible to list all the embodiments here. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the claims of the present invention.

Claims

1. A photogrammetric point cloud adjustment processing method based on satellite photon point cloud, characterized in that: include: Acquire photon point cloud data and photogrammetry point cloud data; Dividing the photon point cloud data into a two-dimensional grid, and selecting photon points that meet preset conditions from the grid as first control points; Calculating the terrain relief similarity between the first control point and the photogrammetric point cloud data; According to the terrain undulation similarity, using a simulated annealing algorithm to calculate the offset between the first control point and the photogrammetric point cloud data; According to the offset, the first control point is registered, and second control points meeting preset conditions are selected from the registered first control points; Get the original RPC parameters of photogrammetric point cloud data; According to the second control point, the original RPC parameters are iteratively adjusted and optimized to obtain RPC parameters that meet the accuracy conditions; According to the RPC parameters that meet the accuracy conditions, the adjusted photogrammetric point cloud is calculated.

2. The photogrammetric point cloud adjustment processing method based on satellite photon point cloud according to claim 1, characterized in that: Dividing the photon point cloud data into a two-dimensional grid, and selecting photon points that meet preset conditions from the grid as control points, including: Divide the photon point cloud data into a two-dimensional grid with the y-axis as the horizontal coordinate and the z-axis as the vertical coordinate; Count the number of photon points in each grid , the number of photon points Reaching the threshold The photon points in the grid are retained and the corresponding matrix is ​​constructed , where the retained grid is in the matrix The corresponding position is 1, and the rest are 0. The specific expression is as follows: in, i Indicates the vertical position of the grid and is used to traverse or reference different rows in the matrix. j Indicates the horizontal position of the grid and is used to traverse or reference different columns in the matrix; According to the matrix Calculate the continuity score for each raster , whose expression is as follows: Group the rasters in each column and keep the rasters with the largest continuity score and a continuity score greater than the threshold in each raster group. Photon points in a grid As a control point, its expression is as follows: 。 3. The photogrammetric point cloud adjustment processing method based on satellite photon point cloud according to claim 1, characterized in that: Calculate the elevation difference between control points and photogrammetric point cloud data As terrain relief similarity, it includes: The photogrammetric points and control points The elevation is set to 0, the horizontal and vertical coordinates remain unchanged, and saved to the 2D photogrammetry point cloud carrier and 2D photon point cloud carrier The expression is as follows: ； For two-dimensional photon point cloud carrier For each control point in , find the 2D photogrammetric point cloud vector The photogrammetric point with the smallest Euclidean distance is taken as the corresponding point cloud , whose expression is as follows: in, Indicates the Euclidean distance between the calculated control point and the photogrammetric point; The absolute value of the elevation difference between the control point and its corresponding photogrammetric point is calculated as the elevation difference, and its expression is as follows: in, and are the elevations of the control point and its corresponding photogrammetric point, respectively.

4. The photogrammetric point cloud adjustment processing method based on satellite photon point cloud according to claim 3, characterized in that: According to the terrain undulation similarity, the offset between the control point and the photogrammetric point cloud data is calculated using a simulated annealing algorithm, including: Change the spatial position of the control point within a preset search area; Get the current temperature of the simulated annealing algorithm , whose expression is as follows: in, n represents the number of iterations, is the cooling rate; According to the current temperature, the control point is calculated x and y The offset in the direction is expressed as follows: in, for A random number; Update the coordinates of the control point according to the offset of the control point in the x and y directions; Calculate the elevation difference change of the control point based on the updated coordinates , whose expression is as follows: in, Represents the elevation difference between the control point and the photogrammetric point. n Indicates the number of iterations; Determine whether the elevation difference change meets the preset conditions. If not, iterate the above steps. If yes, output the current accumulated offsets in the x and y directions as the offsets of the optimal position of the control point on the horizontal plane according to the following formula: 。 5. The photogrammetric point cloud adjustment processing method based on satellite photon point cloud according to claim 4, characterized in that: After obtaining the offset of the optimal position of the control point on the horizontal plane, the method further includes: The calculation control point and the corresponding photogrammetric point are x , y and z Offset in direction , whose expression is as follows: in, represents the total number of photogrammetric point clouds, No. i control points x coordinate, For the i The control points correspond to the photogrammetric points x coordinate, No. i control points y coordinate, For the i The control points correspond to the photogrammetric points y coordinate, No. i control points z coordinate, For the i The control points correspond to the photogrammetric points z coordinate.

6. The photogrammetric point cloud adjustment processing method based on satellite photon point cloud according to claim 5, characterized in that: According to the second control point, the original RPC parameters are iteratively adjusted and optimized to obtain RPC parameters that meet the accuracy conditions, including: Step a: Construct the RFM model of the original RPC parameters; Step b: Performing adjustment processing of affine transformation compensation on the RFM model; Step c: according to the second control point, utilizing terrain-independent algorithm successively to carry out RPC parameter update to the RFM model processed by the adjustment of affine transformation compensation, and carrying out the adjustment of Fourier compensation; Step d: Determine whether the RPC parameters of the RFM model meet the preset accuracy conditions. If so, output the RPC parameters. If not, jump to step c.

7. The photogrammetric point cloud adjustment processing method based on satellite photon point cloud according to claim 6 is characterized in that: According to the second control point, the RPC parameters of the RFM model adjusted by affine transformation compensation are updated by using a terrain-independent algorithm, including: According to the second control point, grid points are extracted at equal intervals in the image plane, and the surface is divided into a plurality of layers according to the elevation range; According to the image coordinates of each grid point, the object coordinates of each grid point on each surface layer are calculated using the RFM model; Construct a virtual control grid based on the object coordinates of each grid point on each surface layer; Based on the virtual control grid, the RPC parameters of the RFM model are updated using the least squares method.

8. The photogrammetric point cloud adjustment processing method based on satellite photon point cloud according to claim 5, characterized in that: After obtaining the adjusted photogrammetric point cloud coordinates, the method further includes: Calculate the average Euclidean distance between the preset verification point cloud and the adjusted photogrammetric point cloud. The expression is as follows: in, is the average Euclidean distance in the x direction between the preset verification point cloud and the adjusted photogrammetric point cloud, is the average Euclidean distance in the y direction between the preset verification point cloud and the adjusted photogrammetric point cloud, is the average Euclidean distance in the z direction between the preset verification point cloud and the adjusted photogrammetric point cloud, The preset verification point cloud i Points x coordinate, The preset verification point cloud i Points y coordinate, The preset verification point cloud i Points z coordinate; The accuracy of the photogrammetric point cloud after adjustment was evaluated based on the average Euclidean distance.

9. The photogrammetric point cloud adjustment processing method based on satellite photon point cloud according to claim 4, characterized in that: Determining whether the elevation difference change amount meets a preset condition includes: Determine the current elevation difference change Is it less than the elevation difference change of the previous iteration? And the current temperature Is it less than the threshold, or , if so, then output the current accumulated offsets in the x and y directions as the offsets of the optimal position of the control point on the horizontal plane.

10. A photogrammetry point cloud adjustment processing system based on satellite photon point cloud, characterized in that: include: A first acquisition module acquires photon point cloud data and photogrammetry point cloud data; A first screening module divides the photon point cloud data into a two-dimensional grid, and screens out photon points that meet preset conditions from the grid as first control points; A first calculation module calculates the terrain relief similarity between the first control point and the photogrammetric point cloud data; A second calculation module calculates the offset between the first control point and the photogrammetric point cloud data using a simulated annealing algorithm according to the terrain undulation similarity; A second screening module performs registration processing on the first control points according to the offset, and screens out second control points that meet preset conditions from the registered first control points; The second acquisition module acquires the original RPC parameters of the photogrammetric point cloud data; An optimization module, according to the second control point, performs iterative adjustment optimization on the original RPC parameters to obtain RPC parameters that meet the accuracy condition; The third calculation module calculates the adjusted photogrammetric point cloud according to the RPC parameters that meet the accuracy conditions.

Images