A random sampling adjustment method for laser altimetry data and optical stereo images

Abstract

The present invention discloses a method for jointly processing laser altimetry data and optical stereo images by random sampling adjustment, and relates to the technical field of photogrammetry and remote sensing image processing. The laser altimetry points are preliminarily screened, the same-name points of the laser points between the stereo images are obtained, the laser elevation is ignored for uncontrolled regional network adjustment, a number of laser points are randomly selected as elevation control points for adjustment, the above steps are repeated, and a group of laser points with the largest number of elevation differences less than a threshold is retained as the maximum consistent set, and the maximum consistent set is used as the elevation control point for final adjustment. The present invention can automatically eliminate the influence of the uncertainty of the laser point elevation caused by complex terrain and objects, and the inconsistency between the laser point and the stereo image caused by the change of terrain and objects over time, and adaptively obtain the best coupling result of the laser altimetry data and the optical stereo image, so as to ensure the geometric accuracy and robustness of the final joint processing result.

Description

Technical Field

[0001] The invention belongs to the technical field of photogrammetry and remote sensing image processing, and in particular relates to a method for jointly processing laser altimetry data and optical stereo images by random sampling adjustment. Background Art

[0002] High-precision three-dimensional terrain data such as DEM and DSM play an important role in national production and construction. Satellite stereo mapping is one of the main ways to produce large-scale high-precision three-dimensional terrain data. With the steady progress of three-dimensional real scene construction at the national level, the demand for producing high-precision three-dimensional terrain data using high-resolution satellite stereo images is becoming more and more vigorous. At present, the ground resolution of domestic optical stereo satellites has been continuously improved and has reached the sub-meter level, but the accuracy of optical satellite image stereo mapping, especially the elevation accuracy, has not been improved to the same extent. The main reason is that it is affected by the satellite attitude and orbit measurement accuracy, the manufacturing level of various sensor hardware, and atmospheric conditions. Satellite laser altimetry can obtain high-precision surface elevation data, which is an important supplement to make up for the insufficient elevation accuracy of optical remote sensing satellite stereo measurement. Therefore, improving the geometric positioning accuracy of remote sensing images through satellite laser altimetry has become a development trend of satellite photogrammetry. The laser altimeter carried by the "Resource-3" 02 / 03 satellite has increased the accuracy of optical stereo image uncontrolled elevation measurement from 5m to about 3m, meeting the 1:50000 scale mapping accuracy requirements. The Gaofen-7 satellite is equipped with a two-line array sub-meter resolution stereo camera and a two-beam laser altimeter. It uses laser altimetry data to assist optical imaging for stereo mapping, which is used for the production of high-precision DSM / DEM at a scale of 1:10,000.

[0003] Laser altimetry is used to improve the geometric positioning accuracy of high-resolution stereo images, mainly by using laser elevation points as elevation control for the regional network adjustment of stereo images. The traditional regional network adjustment method only requires a small amount of control data, but the absolute accuracy of the control data itself is extremely high. However, it is difficult to guarantee the absolute elevation accuracy of laser altimetry data due to the influence of terrain and objects. Usually, the accuracy is good only on flat hard surfaces, resulting in uneven distribution of laser points that meet the accuracy conditions. It is difficult to obtain laser points that meet the accuracy conditions in many areas. At the same time, it is difficult to distinguish between laser points with good and poor accuracy by relying solely on laser altimetry data itself, which greatly limits the applicability of laser altimetry data as control data in the process of satellite stereo mapping. Therefore, it is urgent to explore a more adaptable joint processing method of laser altimetry data and optical stereo images.

[0004] At present, satellite laser altimetry technology is developing in the direction of multi-beam high repetition rate, which can provide a sufficient number of laser points. The ICESat-2 satellite launched by NASA in 2018 uses single-photon detection technology, with a laser repetition rate of 10KHz and a ground interval of only 0.7m between laser points. The Land Ecosystem Carbon Monitoring Satellite (CM-1) launched by my country in 2022 uses a multi-beam linear detection system with a laser repetition rate of 40Hz. Therefore, it is necessary to include as many laser points as possible in the satellite stereo image regional network adjustment, give full play to the advantages of laser altimetry technology, and expand the scope of application of laser altimetry data and optical remote sensing images. To this end, a joint processing method of laser altimetry data and optical stereo images with random sampling adjustment is proposed to solve the bottleneck problem of laser altimetry data-assisted high-resolution optical stereo image automated joint mapping. Summary of the invention

[0005] The present invention aims to provide a random sampling adjustment laser altimetry data and optical stereo image joint processing method, through the random sampling and adjustment technology, to achieve accurate joint processing of laser altimetry data and optical stereo image, to improve the elevation accuracy.

[0006] In order to solve the above technical problems, the present invention is achieved through the following technical solutions:

[0007] The present invention is a method for jointly processing laser altimetry data and optical stereo images by random sampling adjustment, comprising the following steps:

[0008] S1, preliminary selection of laser height measurement points: preliminary screening of laser height measurement points, and elimination of laser points with obviously uncertain elevations;

[0009] S2, obtaining the same-name points: calculating the same-name points of the laser point between each stereo image according to the laser point coordinates and the geometric positioning model of the stereo image;

[0010] S3, uncontrolled regional network adjustment: without considering the laser elevation, uncontrolled regional network adjustment is performed to remove the wrongly matched laser points;

[0011] S4, random sampling adjustment: randomly select several laser points as elevation control points for adjustment, use the adjustment results to calculate the elevations of other laser points, and count the laser points whose elevation difference is less than the threshold as the consistent set;

[0012] S5, repeatedly screening the maximum consistent set: repeat S4 several times, each time retaining a group of laser points with the largest number of elevation differences less than the threshold as the current maximum consistent set, until the final maximum consistent set with the largest number is found;

[0013] S6, final adjustment: Use the final maximum consistent set as the elevation control point for adjustment to obtain the final adjustment result.

[0014] As a preferred technical solution of the present invention, S2 comprises the following steps:

[0015] S21, calculating the image point coordinates of the laser point on the stereo image according to the laser point coordinates and the geometric positioning model of the stereo image, and obtaining the ground object image of the area near the laser point;

[0016] S22, for each laser point, taking a certain ground object image as a reference, performing feature matching with other ground object images to obtain the geometric transformation relationship between the ground object images;

[0017] S23, calculating the image point coordinates of the laser point in other ground object images according to the geometric transformation relationship;

[0018] S24, calculating the image point coordinates of the laser point on the original image;

[0019] S25, matching a certain number of connection points according to requirements.

[0020] As a preferred technical solution of the present invention, the geometric positioning model used in S21 is an RPC model, which is abbreviated as:

[0021] As a preferred technical solution of the present invention, S3 comprises the following steps:

[0022] S31, using the original RPC parameters, forward intersection calculates the object coordinates of all laser points and connection points;

[0023] S32, combining the desired object coordinates and the corresponding image point coordinates, constructing an error equation point by point according to the connection points;

[0024] S33, all error equations are combined and the RPC compensation parameters and object coordinates are solved by the least square method;

[0025] S34, taking into account the RPC compensation parameters, calculating the image square residuals of each laser point and the connection point, and eliminating the image points whose residuals are greater than the threshold.

[0026] As a preferred technical solution of the present invention, S4 includes:

[0027] S41, obtaining a number of mutually unequal random numbers, and selecting some laser points as elevation control points according to the random numbers;

[0028] S42, construct the error equation of the elevation control point;

[0029] S43, constructing error equations for other laser points and connection points;

[0030] S44, combine all error equations to solve RPC compensation parameters and object coordinates;

[0031] S45, calculating the absolute difference of elevations of other laser points;

[0032] S46, counting the laser points whose absolute elevation difference is less than a threshold value as the current consistent set, and recording the number of laser points in the point set.

[0033] As a preferred technical solution of the present invention, S5 includes:

[0034] S51, repeat step S4, compare the number of points in the current consistent set with the maximum consistent set, and update the maximum consistent set;

[0035] S52, when the number of repetitions reaches the limit I max When , step S5 ends.

[0036] As a preferred technical solution of the present invention, the limit value I in S52 max The calculation is shown as

[0037]

[0038] In the formula, w i is the probability that the absolute difference of the laser point elevation is less than the threshold after the i-th iteration, that is, p is the preset confidence level, which is generally between 0.95 and 0.98, and K is the number of randomly selected laser elevation control points.

[0039] As a preferred technical solution of the present invention, S6 includes:

[0040] S61, construct the elevation control point error equation based on the maximum consistent set of laser points;

[0041] S62, constructing a connection point error equation based on other laser points and the connection point;

[0042] S63, combining the error equations in the above two steps to solve the final RPC compensation parameters and object space coordinates as the regional network adjustment results.

[0043] The present invention has the following beneficial effects:

[0044] The main advantage is that it can automatically eliminate the influence of uncertainty in laser point elevation caused by complex terrain and objects and the inconsistency between laser points and stereo images caused by changes in terrain and objects over time, adaptively obtain the best coupling result of laser altimetry data and optical stereo images, and ensure the geometric accuracy and robustness of the final joint processing results.

[0045] Of course, any product implementing the present invention does not necessarily need to achieve all of the advantages described above at the same time. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings required for describing the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other accompanying drawings can be obtained based on these accompanying drawings without paying creative work.

[0047] Figure 1 A flow chart according to the present invention;

[0048] Figure 2 It is a schematic diagram of preliminary screening of laser points during the specific implementation of the present invention;

[0049] Figure 3 It is a schematic diagram of the laser point and the DSM range of the stereoscopic image in the specific implementation process of the present invention;

[0050] Figure 4 A schematic diagram of matching points of laser points on different stereo images during the specific implementation of the present invention;

[0051] Figure 5 It is a schematic diagram of laser elevation control points and laser point consistency sets of a certain random sampling adjustment in a specific implementation process of the present invention;

[0052] Figure 6 It is a schematic diagram of the maximum consistent set of laser elevation control points and laser points during the specific implementation of the present invention. DETAILED DESCRIPTION

[0053] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention. Specific embodiment:

[0055] See also Figure 1-6 As shown, a specific application of this embodiment is: comprising the following steps:

[0056] Step S1, preliminarily screening the laser height measurement points and eliminating the laser points with uncertain heights;

[0057] Affected by complex terrain and objects, the elevation obtained by satellite laser altimetry is uncertain, and it is impossible to accurately distinguish whether the measured elevation is the ground elevation value. Therefore, it is necessary to preliminarily eliminate the laser points with uncertain elevation. For the laser points obtained by the linear laser altimetry technology, the laser points with uncertain elevation are eliminated according to the number of waveform Gaussian decompositions and waveform characteristic parameters. Figure 1As shown in the figure, laser points with a single full-waveform laser echo peak and a Gaussian distribution usually have better elevation accuracy, and the relevant parameters are screened according to the laser altimetry data product. For laser points obtained by single-photon detection laser altimetry technology, laser points with uncertain elevation are eliminated according to terrain feature parameters. Figure 2 As shown in Figure 1, the laser altimetry data of the single-photon detection system can provide a large number of surface laser points, which can obtain information such as surface slope and roughness, thereby extracting laser points with better accuracy. The distribution of the selected laser points and the stereo image DSM is shown in Figure 1. Figure 3 shown.

[0058] Step S2, obtaining the same-name points of laser points between the three-dimensional images, and matching additional connection points according to the distribution of laser points;

[0059] Step S21, according to the laser point coordinates and the geometric positioning model of the stereo image, calculate the coordinates of the image points on each stereo image to obtain the ground object image in the area near the laser point. The geometric positioning model of the image usually adopts the RPC model, and its calculation formula is as follows:

[0060]

[0061] Where, NuML, DenL, NumS, DenS are third-order polynomials, and the coefficients of the polynomials are the polynomial coefficients of the known stereo image RPC model. P, L, H are the regularized object coordinates, l, s are the regularized image coordinates, and the calculation formula is as follows:

[0062]

[0063] Where, LAT_OFF, LAT_SCALE, LON_OFF, LON_SCALE, HEIGHT_OFF, HEIGHT_SCALE, SAMP_OFF, SAMP_SCALE, LINE_OFF, LINE_SCALE are all known regularization parameters of the stereo image RPC model. The RPC model can be simply written as:

[0064]

[0065] Step S22, for each laser point, taking a certain ground object image as a reference, feature matching is performed with other ground object images to obtain the geometric transformation relationship between the image and other ground object images;

[0066] Step S23, calculating the coordinates of the laser point in other ground object images according to the geometric transformation relationship;

[0067] Step S24, calculate the image point coordinates of the laser point on the corresponding original image, such as Figure 4 As shown;

[0068] In particular, when the number of laser points used is insufficient or unevenly distributed, a certain number of connection points can be matched on each image and included in the subsequent regional block adjustment calculation to ensure the robustness of the regional network.

[0069] Step S3, ignoring the laser elevation, performing uncontrolled regional network adjustment, and eliminating the laser points that are mismatched;

[0070] Block adjustment uses additional parameters to compensate for positioning errors between images. The compensation model can use first-order or second-order polynomials, where the first-order polynomial compensation model is as follows:

[0071]

[0072] In the formula, a0, a1, a2, b0, b1, and b2 are compensation parameters and the quantities to be solved for the block adjustment. The above block adjustment model is a nonlinear model and needs to be linearized. S and F L Performing a Taylor expansion and taking the first term yields:

[0073]

[0074] The above error equation can be simplified as:

[0075] V1=A1X-B1 (6)

[0076] Where V1 is v x 、v y The residual matrix composed of The error equation coefficient matrix is ​​composed of B1, B2 is the constant matrix composed of the constant terms of the error equation, and P is the weight matrix of the observations, which is generally a unit matrix.

[0077] The regional block adjustment is to solve the RPC model compensation parameters of the image and the object coordinates of the image points. The least squares method is used to solve equation (6):

[0078]

[0079] In particular, when solving the error equation, the forward intersection-backward intersection alternating iteration method can be used to avoid the inversion calculation of a very large matrix.

[0080] Step S31, using the original RPC parameters, forward intersection calculates the object coordinates of all laser points and connection points;

[0081] Step S32, combining the desired object coordinates and the corresponding image point coordinates, and formulating an error equation according to equation (5);

[0082] Step S33, all error equations are combined, and the RPC compensation parameters of the image and the object coordinates of the laser point and the connection point are solved by the least square method according to formula (7);

[0083] Step S34, taking into account the RPC compensation parameters, the image square residuals of each laser point and the connection point are calculated, and the image points with residuals greater than the threshold are eliminated.

[0084] Step S4 randomly selects several laser points and uses them as elevation control points for adjustment. The adjustment result is used to calculate the elevation of other laser points. The laser point set with elevation difference less than the threshold value is counted as the consistent set of laser points and stereo images.

[0085] Step S41, using a random number generation algorithm to obtain K mutually unequal random numbers {s1, s2, s3, ..., s K};

[0086] Step S42: for each random number k∈{s1,s2,s3,…,s K}, take the kth laser point as the elevation control point, and the linearization error equation is as follows:

[0087]

[0088] The above error equation can be simplified as:

[0089] V2=A2X-B2 (9)

[0090] Step S43, for other laser points and connection points, the error equation is still listed according to formula (5);

[0091] Step S44, combine all error equations and use the least square method to solve the RPC compensation parameters of the image and the object coordinates of the laser point and the connection point:

[0092] X=(A T PA) -1 A T PB (10)

[0093] In the formula,

[0094] Step S45, using the RPC of the image and its compensation parameters, calculate the elevation H' of other laser points, and calculate the absolute difference with the actual elevation H of the laser point:

[0095] ΔH=|HH′| (11)

[0096] Step S46, the laser points whose elevation difference is less than the threshold h are counted as the consistent set of laser points and stereo images, which is denoted as C. The consistent set of laser elevation control points and laser points is as follows: Figure 5 shown.

[0097] Step S5, repeat step S4 several times, and retain the set of laser points with the largest number of consistent sets of laser points and stereo images as the maximum consistent set C max .

[0098] Step S51, repeat step S4, compare the laser point obtained by the current repetition number i with the stereo image consistent set C i With the current maximum consistent set C max The number of i >C max , then update C max , that is, C max =C i , and update the repetition limit I max :

[0099] The calculation is shown as

[0100]

[0101] In the formula, w i is the probability that the absolute difference of the laser point elevation is less than the threshold after the i-th iteration, that is, p is the preset confidence level, which is generally between 0.95 and 0.98, and K is the number of randomly selected laser elevation control points.

[0102] If C i ≤C max , then keep C max , and repeat step S4.

[0103] Step S52: when the current number of repetitions is greater than the limit, that is, i>I max , and the repetitive steps are ended. After several iterations, the maximum consistent set of laser elevation control points and laser points is as follows Figure 6 shown.

[0104] Step S6, using the obtained maximum consistent set C max The laser points in are used as elevation control points for adjustment to obtain the final adjustment result.

[0105] Step S61: the maximum consistent set C max The error equation for the laser point in is listed according to formula (8);

[0106] Step S62, for other laser points and connection points, error equations are listed according to formula (5);

[0107] Step S63, all error equations are combined and solved by least square method according to equation (10) as the final result of regional block adjustment.

[0108] The beneficial effect of the present invention is that a random sampling adjustment laser altimetry data and optical stereo image joint processing method is realized. The main advantage of this method is that it can automatically eliminate the influence of the uncertainty of the laser point elevation caused by complex terrain and objects and the inconsistency between the laser point and the stereo image caused by the change of terrain and objects over time, adaptively obtain the best coupling result of the laser altimetry data and the optical stereo image, and ensure the geometric accuracy and robustness of the final joint processing result.

[0109] In the description of this specification, the description with reference to the terms "one embodiment", "example", "specific example", 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 invention. In this specification, the schematic representation of the above terms does not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described can be combined in any one or more embodiments or examples in a suitable manner.

[0110] The preferred embodiments of the present invention disclosed above are only used to help illustrate the present invention. The preferred embodiments do not describe all the details in detail, nor do they limit the invention to the specific implementation methods described. Obviously, many modifications and changes can be made according to the content of this specification. This specification selects and specifically describes these embodiments in order to better explain the principles and practical applications of the present invention, so that those skilled in the art can understand and use the present invention well. The present invention is limited only by the claims and their full scope and equivalents.

Claims

1. A method for jointly processing laser altimetry data and optical stereo images by random sampling adjustment, characterized in that: The following steps are involved: S1, preliminary selection of laser height measurement points: preliminary screening of laser height measurement points, and elimination of laser points with obviously uncertain elevations; S2, obtaining the same-name points: calculating the same-name points of the laser point between each stereo image according to the laser point coordinates and the geometric positioning model of the stereo image; S3, uncontrolled regional network adjustment: without considering the laser elevation, uncontrolled regional network adjustment is performed to remove the wrongly matched laser points; S4, random sampling adjustment: randomly select several laser points as elevation control points for adjustment, use the adjustment results to calculate the elevations of other laser points, and count the laser points whose elevation difference is less than the threshold as the consistent set; S5, repeatedly screening the maximum consistent set: repeat S4 several times, each time retaining a group of laser points with the largest number of elevation differences less than the threshold as the current maximum consistent set, until the final maximum consistent set with the largest number is found; The S5 includes: S51, repeat step S4, compare the number of points in the current consistent set with the maximum consistent set, and update the maximum consistent set; S52, when the number of repetitions reaches the limit I max When , end step S5; The limit value I in S52 max The calculation is shown as In the formula, w i is the probability that the absolute difference of the laser point elevation is less than the threshold after the i-th iteration, that is, p is the preset confidence level, which ranges from 0.95 to 0.98, and K is the number of randomly selected laser height control points; S6, final adjustment: Use the final maximum consistent set as the elevation control point for adjustment to obtain the final adjustment result.

2. The laser altimetry data and optical stereo image joint processing method of a random sampling adjustment according to claim 1, characterized in that: The S2 comprises the following steps: S21, calculating the image point coordinates of the laser point on the stereo image according to the laser point coordinates and the geometric positioning model of the stereo image, and obtaining the ground object image of the area near the laser point; S22, for each laser point, taking a certain ground object image as a reference, performing feature matching with other ground object images to obtain the geometric transformation relationship between the ground object images; S23, calculating the image point coordinates of the laser point in other ground object images according to the geometric transformation relationship; S24, calculating the image point coordinates of the laser point on the original image; S25, matching a certain number of connection points according to requirements.

3. The laser altimetry data and optical stereo image joint processing method of a random sampling adjustment according to claim 2, characterized in that: The geometric positioning model used in S21 is the RPC model, which is abbreviated as:

4. The laser altimetry data and optical stereo image joint processing method of a random sampling adjustment according to claim 3, characterized in that: The S3 comprises the following steps: S31, using the original RPC parameters, forward intersection calculates the object coordinates of all laser points and connection points; S32, combining the desired object coordinates and the corresponding image point coordinates, constructing an error equation point by point according to the connection points; S33, all error equations are combined and the RPC compensation parameters and object coordinates are solved by the least square method; S34, taking into account the RPC compensation parameters, calculating the image square residuals of each laser point and the connection point, and eliminating the image points whose residuals are greater than the threshold.

5. The laser altimetry data and optical stereo image joint processing method of a random sampling adjustment according to claim 2, characterized in that: The S4 includes: S41, obtaining a number of mutually unequal random numbers, and selecting some laser points as elevation control points according to the random numbers; S42, constructing the error equation of the elevation control point; S43, constructing error equations for other laser points and connection points; S44, combine all error equations to solve RPC compensation parameters and object coordinates; S45, calculating the absolute difference of elevations of other laser points; S46, counting the laser points whose absolute elevation difference is less than a threshold value as the current consistent set, and recording the number of laser points in the point set.

6. The laser altimetry data and optical stereo image joint processing method of a random sampling adjustment according to claim 2, characterized in that: The S6 comprises: S61, construct the error equation of the elevation control point based on the maximum consistent set of laser points; S62, constructing a connection point error equation based on other laser points and the connection point; S63, combining the error equations in the above two steps to solve the final RPC compensation parameters and object space coordinates as the regional network adjustment results.

Images