Photogrammetric Point Cloud Adjustment Processing Method and System Based on Satellite Photon Point Cloud

Abstract

The present invention relates to the technical field of photogrammetry, and discloses a photogrammetric point cloud adjustment processing method and system based on satellite photon point clouds. The method includes obtaining photon point cloud data and photogrammetric point cloud data; dividing the photon point cloud data into grids of two-dimensional planes and screening out first control points; calculating the terrain undulation similarity between the first control points and the photogrammetric point cloud data, and using the simulated annealing algorithm to calculate the offset between the first control points and the photogrammetric point cloud data; according to the offset, performing registration processing on the first control points, and screening out second control points that meet the preset conditions from the registered first control points; obtaining and iteratively adjusting and optimizing the original RPC parameters of the photogrammetric point cloud data to obtain RPC parameters that meet the accuracy conditions; according to the RPC parameters that meet the accuracy conditions, calculating the adjusted photogrammetric point cloud. The present invention has the advantages of high universality and improves the accuracy of the photogrammetric point cloud.

Description

Technical Field

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

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

[0003] Currently, some scholars have attempted to jointly process satellite image data and photon radar data to establish the corresponding relationship 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 laser points as corresponding image points. Other scholars have filtered the photon point cloud data and selected points with qualified accuracy as control points for adjustment. However, these methods have some problems: the former has inaccurate laser point object-image correspondence due to inevitable errors in the positioning model, especially with larger deviations in mountainous forest areas with rugged terrain. The latter has low universality for mountainous forest areas with rugged terrain due to the overly single registration method, 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 accuracy existing in the prior art, the present invention proposes the following technical solutions:

[0005] In a first aspect, the present invention proposes a method for adjusting and processing photogrammetric point clouds based on satellite photon point clouds, including:

[0006] Obtain photon point cloud data and photogrammetric point cloud data.

[0007] Divide the photon point cloud data into grids in a two-dimensional plane, and screen out photons that meet the preset conditions from the grids as the first control points.

[0008] Calculate the terrain undulation similarity between the first control points and the photogrammetric point cloud data.

[0009] According to the terrain undulation similarity, use the simulated annealing algorithm to calculate the offset between the first control points and the photogrammetric point cloud data.

[0010] According to the offset, perform registration processing on the first control points, and screen out the second control points that meet the preset conditions from the registered first control points.

[0011] Obtain the original RPC parameters of the photogrammetric point cloud data.

[0012] According to the second control points, perform iterative adjustment and optimization on the original RPC parameters to obtain RPC parameters that meet the accuracy conditions.

[0013] Calculate the adjusted photogrammetric point cloud based on the RPC parameters that meet the accuracy conditions.

[0014] As a preferred technical solution, divide the photon point cloud data into grids of a two-dimensional plane, and screen out photon points that meet the preset conditions from the grids, including:

[0015] Taking the y-axis as the abscissa and the z-axis as the ordinate, divide the photon point cloud data into grids of a two-dimensional plane.

[0016] Count the number of photon points in each grid , and the number of photon points reaching the threshold retain the photon points in the grid and construct the corresponding matrix , where the retained grids are 1 at the corresponding positions in the matrix , and the rest are 0. The specific expression is as follows:

[0017]

[0018]

[0019] Among them, i represents the position of the grid in the vertical direction, used to traverse or reference different rows in the matrix, j represents the position of the grid in the horizontal direction, used to traverse or reference different columns in the matrix.

[0020] Calculate the continuity score of each grid according to the matrix , and its expression is as follows:

[0021]

[0022] Set the grids in each column as a group, and retain the photon points in the grid with the largest continuity score and a continuity score greater than the threshold in each grid group as the control points, and its expression is as follows:

[0023]

[0024] .

[0025] As a preferred technical solution, calculate the elevation difference between the control points and the photogrammetric point cloud data As the terrain undulation similarity, it includes:

[0026] Respectively set the elevations of the photogrammetric points and the control points to 0, keep the abscissa and ordinate unchanged, and save them to the two-dimensional photogrammetric point cloud carrier and the two-dimensional photon point cloud carrier The expression is as follows:

[0027] ;

[0028]

[0029] For each control point in the two-dimensional photon point cloud carrier , find the photogrammetric point with the smallest Euclidean distance in the two-dimensional photogrammetric point cloud carrier as the corresponding point cloud , and its expression is as follows:

[0030]

[0031]

[0032] Among them, represents the calculation of the Euclidean distance between the control point and the photogrammetric point.

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

[0034]

[0035] Among them, and are the elevations of the control point and its corresponding photogrammetric point respectively.

[0036] As a preferred technical solution, according to the terrain undulation similarity, use the simulated annealing algorithm to calculate the offset between the control points and the photogrammetric point cloud data, including:

[0037] Change the spatial position of the control point within the preset search area.

[0038] Obtain the current temperature of the simulated annealing algorithm , and its expression is as follows:

[0039]

[0040] Among them, nrepresents the number of iterations, is the cooling rate.

[0041] According to the current temperature, calculate the offsets of the control points in the x and y directions, and its expression is as follows:

[0042]

[0043]

[0044] where, is a random number.

[0045] Update the coordinates of the control points according to the offsets of the control points in the x and y directions.

[0046] Calculate the change in the elevation difference of the control points according to the updated coordinates , and its expression is as follows:

[0047]

[0048] where, represents the elevation difference between the control point and the photogrammetric point, n represents the number of iterations.

[0049] Judge whether the change in the elevation difference satisfies a preset condition. If not, iterate and execute the above steps. If so, according to the following formula, output the cumulative offsets in the x and y directions as the offsets of the optimal position of the control point on the horizontal plane:

[0050]

[0051] .

[0052] As a preferred technical solution, after obtaining the offsets of the optimal position of the control point on the horizontal plane, the method further includes:

[0053] Calculate the offsets between the control point and the corresponding photogrammetric point in the x , y and z directions , and its expression is as follows:

[0054]

[0055]

[0056]

[0057] where, Indicates the total number of photogrammetric point clouds, The i coordinates of the x nth control point, are the i coordinates of the photogrammetric points corresponding to the x nth control point, The i coordinates of the y nth control point, are the i coordinates of the photogrammetric points corresponding to the y nth control point, The i coordinates of the z nth control point, are the i coordinates of the photogrammetric points corresponding to the z nth control point.

[0058] As a preferred technical solution, according to the second control point, iterative adjustment optimization is performed on the original RPC parameters to obtain RPC parameters that meet the accuracy conditions, including:

[0059] Step a: Construct an RFM model of the original RPC parameters.

[0060] Step b: Perform adjustment processing with affine transformation compensation on the RFM model.

[0061] Step c: According to the second control point, sequentially use terrain-independent algorithms to update the RPC parameters of the RFM model that has undergone adjustment processing with affine transformation compensation, and perform adjustment processing with Fourier compensation.

[0062] 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 for execution.

[0063] As a preferred technical solution, according to the second control point, use terrain-independent algorithms to update the RPC parameters of the RFM model that has undergone adjustment processing with affine transformation compensation, including:

[0064] According to the second control point, equally spaced grid points are extracted in the image plane, and the ground surface is divided into several layers according to the elevation range.

[0065] According to the image coordinates of each grid point, use the RFM model to calculate the object coordinates of each grid point on each layer of the ground surface.

[0066] According to the object coordinates of each grid point on each layer of the ground surface, construct a virtual control grid.

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

[0068] As a preferred technical solution, after obtaining the adjusted photogrammetric point cloud coordinates, the method further includes:

[0069] Calculating the average Euclidean distance between a preset verification point cloud and the adjusted photogrammetric point cloud, and its expression is as follows:

[0070]

[0071]

[0072]

[0073] Wherein, 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, is the i th point in the preset verification point cloud x coordinate, is the i th point in the preset verification point cloud y coordinate, is the i th point in the preset verification point cloud z coordinate.

[0074] According to the average Euclidean distance, evaluate the accuracy of the adjusted photogrammetric point cloud.

[0075] As a preferred technical solution, determining whether the change amount of the elevation difference meets a preset condition includes:

[0076] Determining whether the current change amount of the elevation difference is less than the change amount of the elevation difference in the previous iteration and whether the current temperature is less than the threshold value, or , if so, output the cumulative offsets in the x and y directions as the offsets of the optimal position of the control point on the horizontal plane.

[0077] In a second aspect, the present invention also proposes a photogrammetric calibration system based on satellite photon point clouds, which is applied to the photogrammetric point cloud adjustment processing method based on satellite photon point clouds described in any one of the solutions in the first aspect, and includes:

[0078] The first acquisition module acquires photon point cloud data and photogrammetric point cloud data.

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

[0080] The first calculation module calculates the terrain undulation similarity between the first control points and the photogrammetric point cloud data.

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

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

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

[0084] The optimization module performs iterative adjustment and optimization on the original RPC parameters according to the second control points to obtain RPC parameters that meet the accuracy conditions.

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

[0086] The beneficial effects of the present invention at least include:

[0087] By dividing the photon point cloud data into grids on a two-dimensional plane and screening out photons that meet the preset conditions as control points, the present invention effectively utilizes the high-precision characteristics of the ICESat-2 satellite photon point cloud data, providing a reliable reference point for subsequent adjustment processing. Using the terrain undulation similarity as an index and adopting the simulated annealing algorithm to calculate the offset between the control points and the photogrammetric point cloud to be processed, compared with the traditional grid search method, the registration accuracy and universality are greatly improved, especially suitable for mountainous areas with complex terrain. By acquiring the original RPC parameters of the photogrammetric point cloud data and using the selected high-quality second control points to perform iterative adjustment and optimization on these parameters until the accuracy requirements are met. Finally, using the optimized RPC parameters to calculate the adjusted photogrammetric point cloud not only improves the reliability and accuracy of the control points, ensures the high accuracy of the RPC parameters, but also significantly improves the overall accuracy of the final photogrammetric point cloud. Description of the Drawings

[0088] Figure 1 It is a schematic flowchart of the photogrammetric point cloud adjustment processing method based on satellite photon point cloud provided in Embodiment 1.

[0089] Figure 2 Flow chart of the simulated annealing algorithm provided in Embodiment 2.

[0090] Figure 3 Flow chart of the iterative adjustment and optimization of RPC parameters provided in Embodiment 2.

[0091] Figure 4 Schematic diagram of the virtual control grid for calculating RPC parameters provided in Embodiment 2.

[0092] Figure 5 Schematic diagram of the positional relationship between the photon point cloud and the photogrammetric point cloud in the real measurement area provided in Embodiment 3.

[0093] Figure 6 Schematic diagram of the number of control points selected from the photon point cloud provided in Embodiment 3.

[0094] Figure 7 (a) Elevation distribution map of the photon point cloud provided in Embodiment 3, Figure 7 (b) Elevation distribution map of the photogrammetric point cloud provided in Embodiment 3.

[0095] Figure 8 Registration data map of the photon point cloud and the photogrammetric point cloud provided in Embodiment 3.

[0096] Figure 9 Schematic diagram of the accuracy evaluation result provided in Embodiment 3.

[0097] Figure 10 Comparison view of the calibration result of window a provided in Embodiment 3.

[0098] Figure 11 Comparison view of the original data of window b provided in Embodiment 3.

[0099] Figure 12 Comparison view of the local calibration result of mountain area of window a provided in Embodiment 3.

[0100] Figure 13 Comparison view of the local original data of mountain area of window b provided in Embodiment 3.

[0101] Figure 14 Architecture diagram of the photogrammetric calibration system based on satellite photon point cloud provided in Embodiment 4. Specific implementation mode

[0102] The embodiments of the present invention will be described below 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 content disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments. Various details in this specification can also be modified or changed 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 explaining the present invention, rather than limiting the protection scope of the present invention.

[0103] It should be noted that the diagrams provided in the following embodiments only illustrate the basic concept of the present invention in a schematic manner. Therefore, only the components related to the present invention are shown in the diagrams, rather than being drawn according to the number, shape, and size of the components in actual implementation. The type, quantity, and ratio of each component in actual implementation can be arbitrarily changed, and the component layout type may also be more complex.

[0104] In the following description, a large number of details are explored 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.

[0105] Embodiment 1

[0106] This embodiment proposes a photogrammetric point cloud adjustment processing method based on satellite photon point clouds, as Figure 1 shown. Figure 1 is a schematic flowchart of a photogrammetric point cloud adjustment processing method based on satellite photon point clouds provided in this embodiment. The method includes the following steps:

[0107] S1: Obtain photon point cloud data and photogrammetric point cloud data.

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

[0109] S3: Calculate the terrain undulation similarity between the first control points and the photogrammetric point cloud data.

[0110] S4: According to the terrain undulation similarity, use the simulated annealing algorithm to calculate the offset between the first control points and the photogrammetric point cloud data.

[0111] S5: According to the offset, perform registration processing on the first control points, and select the second control points that meet the preset conditions from the registered first control points.

[0112] S6: Obtain the original RPC parameters of the photogrammetric point cloud data.

[0113] S7: According to the second control points, perform iterative adjustment and optimization on the original RPC parameters to obtain RPC parameters that meet the accuracy conditions.

[0114] S8: Calculate the adjusted photogrammetric point cloud according to the RPC parameters that meet the accuracy conditions.

[0115] It can be understood that by dividing the photon point cloud data into two-dimensional plane grids and screening photons that meet the preset conditions as control points, the present invention effectively utilizes the high-precision characteristics of the ICESat-2 satellite photon point cloud data, providing a reliable reference point for subsequent adjustment processing. Using the terrain undulation similarity as an index and adopting the simulated annealing algorithm to calculate the offset between the control points and the photogrammetric points to be processed, compared with the traditional grid search method, the registration accuracy and universality are greatly improved, especially 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 perform iterative adjustment and optimization on these parameters until the accuracy requirements are met. Finally, using the optimized RPC parameters to calculate the adjusted photogrammetric point cloud not only improves the reliability and accuracy of the control points, ensures the high precision of the RPC parameters, but also significantly improves the overall accuracy of the final photogrammetric point cloud.

[0116] The present invention is applicable to various terrain conditions, especially has significant advantages in complex terrain areas lacking control points or without control points, and can be widely applied to fields such as mountain natural resource survey, geological disaster monitoring, and infrastructure construction. The entire processing flow has a high degree of automation, reduces manual intervention, improves processing efficiency, and at the same time reduces the influence 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 high universality, providing an effective technical means for solving photogrammetry problems in complex terrains such as mountainous areas.

[0117] Embodiment 2

[0118] This embodiment makes improvements on the basis of the photogrammetric point cloud adjustment processing method based on satellite photon point cloud proposed in Embodiment 1.

[0119] In this embodiment, dividing the photon point cloud data into grids of two-dimensional plane and screening photons that meet the preset conditions from the grids as control points includes:

[0120] Taking the y-axis as the abscissa and the z-axis as the ordinate, divide the two-dimensional plane of the photon point cloud data into a length of a height of The grid is set according to the minimum and maximum values of the y coordinate and the minimum and maximum values of the z coordinate of the photon point cloud data, and the side length of the grid is set.

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

[0122]

[0123]

[0124] Among them, i represents the position of the grid in the vertical direction and is used to traverse or reference different rows in the matrix. j represents the position of the grid in the horizontal direction and is used to traverse or reference different columns in the matrix.

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

[0126]

[0127] Set the grids in each column as a group, and retain the photon points in the grid with the largest continuity score and a continuity score greater than the threshold as the control points. The expression is as follows:

[0128]

[0129] .

[0130] In this embodiment, the sum of the absolute values of the elevation differences of the nearest points in the plane coordinates is used to evaluate the global terrain similarity. The smaller the sum of the absolute values of the elevation differences, the more similar the terrain undulations are. The smaller.

[0131] The elevation difference between the control point and the photogrammetry point The calculation steps include:

[0132] Respectively, the photogrammetry points and the control points The elevations are all set to 0, the abscissa and ordinate remain unchanged, and they are saved to the two-dimensional photogrammetric point cloud carrier and the two-dimensional photon point cloud carrier Among them, the expression is as follows:

[0133] ;

[0134]

[0135] For each control point in the two-dimensional photon point cloud carrier find the photogrammetric point with the smallest Euclidean distance in the two-dimensional photogrammetric point cloud carrier as the corresponding point cloud The expression is as follows:

[0136]

[0137]

[0138] Among them, represents calculating the Euclidean distance between the control point and the photogrammetric point.

[0139] Calculate the absolute value of the elevation difference between the control point and its corresponding photogrammetric point as the elevation difference value. The expression is as follows:

[0140]

[0141] Among them, and are the elevations of the control point and its corresponding photogrammetric point respectively.

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

[0143]

[0144] Among them, represents the total number of photon point clouds.

[0145] 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 deficiencies of the grid search method in terms of accuracy and the excessive dependence on the normal vector and long calculation time of the ICP registration method. The basic logic of the simulated annealing algorithm is to change the spatial position of the control point to obtain evaluations of different terrain undulation similarities, and finally record the change amount of the control point at the position where the terrain undulation is the most similar, that is, the smallest and , as Figure 2 shownFigure 2 Flow chart of the simulated annealing algorithm provided in this embodiment:

[0146] According to the terrain undulation similarity, use the simulated annealing algorithm to calculate the offset between the control points and the photogrammetric point cloud data, including:

[0147] Change the spatial position of the control points within a preset search area.

[0148] Obtain the current temperature of the simulated annealing algorithm , and its expression is as follows:

[0149]

[0150] Where n represents the number of iterations, is the cooling rate.

[0151] According to the current temperature, calculate the offsets of the control points in the x and y directions, and its expression is as follows:

[0152]

[0153]

[0154] Where is a random number.

[0155] Update the coordinates of the control points according to the offsets of the control points in the x and y directions.

[0156] Calculate the change in elevation difference of the control points according to the updated coordinates , and its expression is as follows:

[0157]

[0158] Where represents the elevation difference between the control point and the photogrammetric point, n represents the number of iterations.

[0159] Judge whether the current change in elevation difference is less than the change in elevation difference of the previous iteration and whether the current temperature is less than the threshold, or , if so, then according to the following formula, output the cumulative offsets in the x and y directions at present as the offsets of the optimal position of the control point on the horizontal plane:

[0160]

[0161] .

[0162] In this embodiment, after obtaining the offset of the optimal position of the control point on the horizontal plane, the method further includes:

[0163] The calculation control point and the corresponding photogrammetric point are x , y and z Offset in direction , whose expression is as follows:

[0164]

[0165]

[0166]

[0167] 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.

[0168] 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:

[0169] Step a: Build the RFM model of the original RPC parameters.

[0170] As an exemplary illustration, the RFM (rational function model) model is further extended on the basis of the traditional Direct Linear Transformation (DLT) model. The transformation from the geodetic coordinates in the object space to the image point coordinates in 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:

[0171]

[0172] Among them, is the normalized image point coordinate in the image plane coordinate system, is the normalized geodetic coordinate in the object space coordinate system. The purpose of the normalization parameters is to ensure that the numerical values of the image coordinates and the object coordinates are both in the range of -1 to 1, avoiding calculation errors caused by differences in the coordinate orders of magnitude. Specifically:

[0173]

[0174]

[0175]

[0176]

[0177]

[0178] Among them, represents the image point coordinate in the image plane coordinate system. represents the object coordinate of the ground point, and can adopt geocentric rectangular coordinates, geodetic coordinates or map projection coordinates. is the normalization parameter of the image coordinates. is the normalization parameter of the object coordinates. 、 、 、 all represent polynomials. Here, a cubic polynomial is introduced. Taking as an example, the form is as shown in the following formula:

[0179]

[0180] Among them, is the polynomial coefficient, that is, the RPC parameter, 、 The constant terms of take 1.

[0181] Step b: Perform adjustment processing for affine transformation compensation on the RFM model.

[0182] Step c: According to the second control point, sequentially update the RPC parameters of the RFM model that has undergone adjustment processing with affine transformation compensation using a terrain-independent algorithm, and perform adjustment processing with Fourier compensation.

[0183] As an exemplary illustration, to correct the positioning error of the RFM model, in this embodiment, without changing the original RPC parameters, a systematic error compensation term for the image coordinates is added, and the compensation term coefficient is used as an adjustment parameter. The observation equation is shown as follows:

[0184]

[0185] Where is the calculated value of the image point coordinates obtained using the original RPC, is the error compensation term of, which can generally be expressed as a mathematical function of the image point observation value . The closer the error compensation term used in the adjustment conforms to the true error distribution, the higher the adjustment accuracy.

[0186] In this embodiment, adjustment is performed through an affine transformation compensation term and a Fourier compensation term. Specifically:

[0187] A first-order polynomial, i.e., an affine transformation formula, can effectively eliminate the systematic error in the image plane. The form of the affine transformation compensation term in this embodiment is shown as follows:

[0188]

[0189] The Fourier series can theoretically approximate any function that meets the continuity condition within a two-dimensional interval. Research shows that the binary Fourier polynomial has a good compensation effect on complex image plane distortions. The form of the Fourier compensation term in this embodiment is shown as follows:

[0190]

[0191] Where, and , and are Fourier polynomial coefficients. and are the cosine component and sine component respectively, that is, shown as follows:

[0192]

[0193]

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

[0195]

[0196] Furthermore, if the error model the error compensation terms in respectively contain compensation term coefficients, considering the unknown object coordinates of the homologous connection points extracted from the stereo images, for linearization, the error equations for the RFM adjustment are as follows:

[0197]

[0198] Writing it in the form of an error equation matrix:

[0199]

[0200] Among them, is the residual vector of the image point coordinate observations. is the coefficient matrix, that is, the partial derivative matrix of the unknowns. is the compensation term coefficient increment vector. is the increment vector of the object space coordinates of the connection points. is the constant term, where is the image point coordinate observation value, is the image point coordinate calculated by substituting the approximate values of the unknowns into Equation . is the weight matrix, which is set as the identity matrix when the accuracies of all observation points are equal. For ground control points, the control point coordinate error equations can be listed. Specifically

[0201]

[0202] Among them, it can be set that , the weight of the control point coordinate observation value

[0203]

[0204] For photon elevation points, the form of their coordinate error equations is the same as that of the general control point coordinate error equations, but their weight values should be determined according to the accuracy of the photon point coordinates, as shown in the following formula:

[0205]

[0206] Among them, , in the data of this embodiment, the mean error of the photon point coordinates is about 3 times that of the control point coordinates. Therefore, 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:

[0207]

[0208] From According to the least squares adjustment principle, establish the normal equation as follows to solve the compensation term coefficients and the object space coordinate increments of the connection points:

[0209]

[0210] Through the above formula, the correction number of the compensation term coefficient and the correction number of the unknown point coordinates can be solved simultaneously.

[0211] 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 for execution.

[0212] In this embodiment, according to the second control points, use the terrain-independent algorithm to update the RPC parameters of the RFM model that has undergone affine transformation compensation and adjustment, including:

[0213] According to the second control points, equally spaced grid points are extracted in the image plane, and the ground surface is divided into several layers according to the elevation range.

[0214] According to the image coordinates of each grid point, use the RFM model to calculate the object space coordinates of each grid point on each layer of the ground surface.

[0215] According to the object space coordinates of each grid point on each layer of the ground surface, construct a virtual control grid.

[0216] Specifically, according to the left information of the second control points, equally spaced grid points are extracted in the image plane. At the same time, the ground surface elevation range is divided into layers. The elevation of the th layer is . For the image coordinates of each grid point, calculate its three-dimensional object space coordinates of the corresponding object points on each layer of the ground surface through the object-image geometric relationship , thereby constructing a virtual control grid with a size of , as shown in Figure 4 . Figure 4 is a schematic diagram of the virtual control grid for calculating RPC parameters.

[0217] The specific steps for calculating the object space coordinates of each grid point on each layer of the ground surface using the RFM model include:

[0218] 1) Assign initial values to the object-space plane coordinates of the point to be solved. , .

[0219] 2) Set the initial step size of the coordinates .

[0220] 3) Calculate the object-space coordinate correction factor , where are the image-space coordinates obtained by substituting the object-space coordinates into Equation , = + d , = + d .

[0221] 4) Calculate the object-space coordinate correction at the current step size: , to obtain the corrected object-space coordinates = + ∆ , = + ∆ .

[0222] 5) Calculate the deviation of the image-space coordinates after correcting the object-space coordinates , where is the image-space coordinate obtained by substituting into Equation .

[0223] 6) Check whether is less than the tolerance (0.01 pixel). If so, output the object-space coordinate result , otherwise use to replace , and reduce the step size , and repeat steps 3) - 6) until the accuracy requirement is met.

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

[0225] It can be understood that in this embodiment, a single adjustment is extended to multiple adjustments, and an iterative adjustment method for RFM with variable RPC is proposed. After performing a single image-space adjustment on the RFM, a virtual control grid is constructed based on the adjustment result and the RPC parameters are recalculated. Since the new RPC parameters have excluded the influence of some systematic errors, the positioning result is closer to the image point observation value, and the next image-space adjustment can be performed on this basis to eliminate the residual systematic errors. After several iterative calculations in this way, the complex image distortion is gradually eliminated.

[0226] Example 3

[0227] Based on the method proposed in Example 2, this example uses simulation experiments to verify the reliability of the RFM iterative adjustment, and after verifying the effectiveness of the adjustment method and determining the optimal number of iterations, it performs survey adjustment on the real survey area and evaluates the accuracy of the method of the present invention.

[0228] The simulation experiment of this example uses RSM to generate 1:10000 scale stereo camera simulation observation data for a pushbroom camera, including the initial values of RPC parameters, the object coordinates of the target points, and the forward and backward image coordinates. By comparing the target positioning accuracy before and after the RFM adjustment, the effectiveness of the adjustment method used is verified. The specific simulation parameters are shown in Table 1. Among them, simulated distortion values that conform to the quadratic polynomial distribution are added to the image coordinates of the target points. For the line elements , angular elements , principal points , principal distances and other parameter "true values", corresponding errors ΔL, ΔA, ΔP, Δf are introduced respectively, and the initial values of RPC parameters are calculated according to the terrain-independent algorithm based on RSM.

[0229] Table 1 Optical satellite stereo image simulation parameters

[0230]

[0231] When ΔL = 100m, ΔA = 100", ΔP = 100 pixels, Δf = 10mm, and the distortion at the image row edge is pixels, the RFM adjustment positioning accuracy for different numbers of iterations is shown in Table 2, where are respectively the image positioning mean errors on the forward and backward images, with the unit of pixel. are respectively the object plane and elevation positioning mean errors, with the unit of meter.

[0232] Table 2 RFM iterative adjustment simulation test positioning results

[0233]

[0234] Due to the small stereo intersection angle of the forward and backward views of the dual-line array camera, large object-space elevation positioning errors are generated by the errors of the interior and exterior orientation elements. It can be seen from Table 2 that when the number of iterations is 1, that is, when using the conventional RFM adjustment based on the affine transformation model, the positioning accuracy of the forward and backward views in the image space is about 3 pixels. At this time, the root mean square error of elevation positioning reaches 8.5 m, exceeding 5 times the root mean square error of the plane, far from meeting the positioning accuracy requirements of the 1:10,000 scale. After the second adjustment, the systematic error in the image space is greatly reduced. The root mean square error of image-space positioning is reduced by 94%, and the root mean square errors of plane and elevation positioning decrease by 93% accordingly. The third adjustment still has an obvious compensation effect on the systematic error. The root mean square error of image-space positioning is further reduced by 22% to within 0.15 pixels, and the root mean square errors of object-space plane and elevation positioning decrease by 28% and 25%. At this time, the adjustment effect is close to the maximum. Increasing the number of iterations further, the positioning accuracies in the image space and object space basically no longer change. Considering the computational performance and cost comprehensively, the number of iterations of RFM adjustment can be set to 3.

[0235] In this embodiment, the southwestern region of Kunming, Yunnan is used as the actual measurement area for adjustment, and the photon point cloud data of ICESat-2 / ATLAS, namely ATL03 and ATL08 on March 6, 2020 in this area, with an along-track distance of 30 km is used as the experimental data. As Figure 5 shown, Figure 5 is a schematic diagram of the position relationship between the photon point cloud and the photogrammetric point cloud in the actual measurement area of this embodiment. The red part is the photon point cloud. The terrain in this area varies greatly, with the highest altitude being 2655.29 m and the lowest being 1413.57 m. The area is mainly covered by forest vegetation. The area of the photogrammetric point cloud for adjustment is about 48, and the number of point clouds is about fifty million. The photon point cloud data of ICESat-2 / ATLAS, namely ATL03 and ATL08 used in this embodiment, are processed and applied after Gaussian coordinate transformation. The parameters of the satellite image experimental data are shown in Table 3.

[0236] Table 3 Satellite Image Experimental Data Parameters

[0237]

[0238] In this embodiment, the corresponding photon point cloud data of ICESat-2 / ATLAS and the unadjusted photogrammetric point cloud data are imported as processing files. Since this experiment involves verifying the accuracy and feasibility of the method, the adjusted photogrammetric point cloud is imported as the evaluation standard, which is not required in actual applications.

[0239] After the files are imported correctly, click Start Debugging and wait for the result output. The first result output from top to bottom is: the strip number, the total number of photon point clouds in this strip, and the number of photon point clouds selected as control points after filtering in this strip. As Figure 6as shown

[0240] As Figure 7 and Figure 8 shown Figure 7 is a comparison chart of the elevation distribution of the photon point cloud and the elevation distribution of the photogrammetric point cloud, where Figure 7 (a) is the elevation distribution chart of the photon point cloud, Figure 7 (b) is the elevation distribution chart of the photogrammetric point cloud, Figure 8 is the registration data chart of the photon point cloud and the photogrammetric point cloud. It can be seen from the image and data that the registration effect using the simulated annealing algorithm and the terrain undulation similarity as the evaluation criterion is very good. The points with a registration difference greater than 10m only account for 2.8% of the total number of points, and the maximum registration difference is 42.6m.

[0241] If a verification file is imported, the accuracy obtained by comparing with the verification file will be output, that is, the offset of the corresponding point cloud between the verification file and the adjusted result. From top to bottom are the offsets in the three directions before and after error correction as Figure 9 shown

[0242] 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, these 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:

[0243] Window a: Displays the point cloud in the verification file (in blue) and the photogrammetric point cloud after adjustment using this method (in red). The purpose of this window is to visually compare the original verification point cloud and the calibrated point cloud to evaluate the calibration effect.

[0244] Figure 10 is the comparison view of the calibration results in window a, highlighting the comparison of the calibrated point cloud (red) and the verification point cloud (blue) in a certain local area. Figure 12 is the comparison view of the local calibration results of the mountain area in window a, used to show the effect of the calibration process in complex terrain.

[0245] Window b: Displays the point cloud in the verification 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 verification point cloud.

[0246] Figure 11 is the comparison view of the original data in window b, showing the comparison of the uncalibrated point cloud (red) and the verification point cloud (blue) in a certain local area. Figure 13It is a comparison view of the original data of the mountain area of window b, further comparing the representation form of the uncalibrated point cloud data in complex terrain and the difference from the verification point cloud.

[0247] In this embodiment, after obtaining the adjusted photogrammetric point cloud coordinates:

[0248] Calculate the average Euclidean distance between the preset verification point cloud and the adjusted photogrammetric point cloud, and its expression is as follows:

[0249]

[0250]

[0251]

[0252] Among them, 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, is the i th point in the preset verification point cloud x coordinate, is the i th point in the preset verification point cloud y coordinate, is the i th point in the preset verification point cloud z coordinate.

[0253] Evaluate the accuracy of the adjusted photogrammetric point cloud according to the average Euclidean distance.

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

[0255]

[0256] It can be seen from the results in Table 4 that the point cloud file after adjustment by the present invention and the verification point cloud file are basically coincident in both mountainous areas and mountain forest and hilly areas, while the coincidence effect of the unadjusted point cloud file and the verification point cloud file is poor. The present invention improves the accuracy by a factor of two in the horizontal direction and by a factor of ten in the vertical direction compared with the existing method, which can effectively prove the feasibility, universality and high accuracy of the present invention.

[0257] Example 4

[0258] As Figure 14As shown in the figure, this embodiment proposes a photogrammetric calibration system based on satellite photon point clouds, which is applied to the photogrammetric point cloud adjustment processing method based on satellite photon point clouds as described in the above embodiment, and includes: 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.

[0259] 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 grids of a two-dimensional plane, and screens out photons that meet the preset conditions from the grids as the first control points. The first calculation module 300 calculates the terrain undulation similarity between the first control points and the photogrammetric point cloud data. The second calculation module 400 calculates the offset between the first control points and the photogrammetric point cloud data using the simulated annealing algorithm according to the terrain undulation similarity. The second screening module 500 performs registration processing on the first control points according to the offset, and screens out second control points that meet the preset conditions from the registered first control points. The second acquisition module 600 acquires the original RPC parameters of the photogrammetric point cloud data. The optimization module 700 performs iterative adjustment and optimization on the original RPC parameters according to the second control points to obtain RPC parameters that meet the accuracy conditions. The third calculation module 800 calculates the adjusted photogrammetric point cloud according to the RPC parameters that meet the accuracy conditions.

[0260] It should be noted that the foregoing explanation of the embodiment of the photogrammetric point cloud adjustment processing method based on satellite photon point clouds also applies to the photogrammetric calibration system based on satellite photon point clouds in this embodiment, and will not be elaborated here.

[0261] In the description of this specification, the descriptions with reference to terms such as "one embodiment", "some embodiments", "example", "specific example", or "some examples" mean that the specific features, structures, materials, or characteristics described in connection with the embodiment or example are included in at least one embodiment or example of the present application. In this specification, the schematic expressions of the above terms are not necessarily directed to the same embodiment or example. Moreover, the specific features, structures, materials, or characteristics described can be combined in a suitable manner in any one or N embodiments or examples. In addition, without conflict, those skilled in the art can combine and combine the different embodiments or examples described in this specification and the features of different embodiments or examples.

[0262] In addition, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the quantity of the indicated technical features. Thus, features defined with "first" and "second" may explicitly or implicitly include at least one such feature. In the description of the present application, the meaning of "N" is at least two, such as two, three, etc., unless otherwise specifically defined.

[0263] Any process or method description shown in a flowchart or described otherwise herein may be understood to represent a module, segment, or portion of code including one or more executable instructions for implementing a customized logical function or process. The scope of the preferred embodiments of the present application includes additional implementations where functions may be executed in a substantially simultaneous manner or in an order opposite to that shown or discussed, according to the functions involved, which should be understood by those skilled in the art to which the embodiments of the present application pertain.

[0264] It should be understood that each part of the present application may be implemented by hardware, software, firmware, or a combination thereof. In the above embodiments, the N steps or methods may 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, any one or a combination of the following techniques well known in the art may be used: discrete logic circuits having logic gate circuits for implementing logical functions on data signals, application specific integrated circuits having appropriate combinational logic gate circuits, programmable gate arrays, field programmable gate arrays, etc.

[0265] Those of ordinary skill in the art of the present technology can understand that all or part of the steps carried by the methods of the above embodiments can be completed by instructing relevant hardware through a program. The program can be stored in a computer-readable storage medium. When the program is executed, it includes one or a combination of the steps of the method embodiments.

[0266] 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 of ordinary skill 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 enumerate all the embodiments here. Any modifications, equivalent replacements, and improvements made within the spirit and principles of the present invention shall be included within 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, 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, retain the photon points in the grid where the number of photon points reaches the threshold and construct the corresponding matrix; Calculate the continuity score for each grid based on the matrix; The grids of each column are set as a group, and in each grid group, the photon point in the grid with the largest continuity score and a continuity score greater than a threshold is retained as the first control point; 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: 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 photogrammetry 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 photogrammetric point cloud adjustment processing system based on satellite photon point cloud, applied to the photogrammetric point cloud adjustment processing method based on satellite photon point cloud as claimed in any one of claims 1 to 9, 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