Method for improving large-area DSM precision by using control points

Abstract

The invention discloses a method for improving large-area DSM precision through control points. The method comprises the steps of overall terrain correction and overall plane and elevation precision improvement. Dividing the target area into regular grids, counting control points in each grid, and determining the elevation correction amount of each grid according to the number of the control points in each grid and the relationship between the longitude and latitude ranges of the control points and corresponding thresholds; carrying out iterative updating on the obtained elevation correction quantity of the grid with the effective elevation correction quantity, and obtaining the elevation correction quantity of the corresponding grid after iterative updating; and on the basis of the elevation correction of the corresponding grid after iteration update, calculating the elevation correction of each row and column of the grid by bilinear interpolation by taking the center of each grid as a reference so as to update the whole DSM. The method has the advantages that on the basis of overall correction, the strategy of combining block processing and iterative correction is adopted, the correction precision can be improved, and meanwhile it can be guaranteed that no edge matching trace caused by block exists in the corrected large-area DSM.

Description

Technical Field

[0001] This invention relates to the field of photogrammetry and remote sensing technology, and in particular to a method for improving the accuracy of large-area DSM using control points. Background Technology

[0002] Various error propagation factors during large-area DSM production often lead to errors in the DSM, such as local or global distortion, deformation, and undulation. These errors can be eliminated through correction using control points, laser points, and reference terrain data. However, due to differences in the type, density, and reliability of control points, it is difficult to achieve ideal correction results by applying a uniform correction globally. Summary of the Invention

[0003] The purpose of this invention is to provide a method for improving the accuracy of large-area DSM by utilizing control points, thereby solving the aforementioned problems existing in the prior art.

[0004] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0005] A method for improving the accuracy of a large-area DSM using control points includes the following steps:

[0006] S1. Overall Terrain Correction: Bilinear interpolation is used to obtain the elevation corresponding to the row and column coordinates of each control point in the target area on the DSM, thus obtaining a second ground coordinate sequence. The control point magnified coordinate sequence and interpolated magnified coordinate sequence obtained based on the ground 3D coordinates in the control point coordinate file and the second ground coordinate sequence are respectively used to obtain a new control point magnified coordinate sequence and interpolated magnified coordinate sequence using the RANSAC method, and the correlation coefficient between the corresponding elevation values ​​of the two new coordinate sequences is calculated. The starting coordinate of the upper left corner of the DSM is modified based on the elevation value combination with the largest correlation coefficient. The affine transformation coefficient between the elevation values ​​of the corresponding two sets of coordinate sequences is calculated based on the least squares method, and the original grid point elevation is replaced with the transformed grid point elevation calculated based on the affine transformation coefficient to update the elevation of all grid points in the entire DSM.

[0007] S2. Correction parameter calculation: Divide the target area into regular grids, count the control points within each grid, and determine the elevation correction amount for each grid based on the number of control points within each grid and the relationship between the latitude and longitude range of the control points and the corresponding thresholds.

[0008] S3. Correction parameter update: Iteratively update the elevation correction of the grid with effective elevation correction obtained in step S2 to obtain the corresponding grid after iterative update.

[0009] S4. Full DSM Update: Based on the elevation correction of the corresponding grid after iterative update, the elevation correction of each row and column of the grid is calculated by bilinear interpolation with the center of each grid as the reference, so as to update the entire DSM.

[0010] Preferably, in step S1, the raster DSM file to be corrected and the control point coordinate file are used as input data. The raster DSM file includes the number of rows R, the number of columns C, the row direction resolution ΔB, the column direction resolution ΔL, and the coordinates of the top left corner (L0, B0); where L0 is the longitude of the top left corner raster point and B0 is the latitude of the top left corner raster point.

[0011] The control point coordinate file contains N ground three-dimensional coordinates (Li', Bi', Hi'); where Li', Bi', and Hi' are the longitude, latitude, and elevation of the i'th control point, respectively; i' = 0, 1, ..., N-1.

[0012] Preferably, step S1 specifically includes the following:

[0013] S11. Set the search window to M*N and the search step size to ΔS. Calculate the row and column coordinates (ri', ci') of each control point on the DSM based on its latitude and longitude. Obtain the corresponding elevation H1i based on the bilinear interpolation of the row and column coordinates, and obtain the second ground coordinate sequence (Li', Bi', H1i'). Multiply the latitude and longitude coordinates of the coordinate sequences (Li', Bi', Hi') and (Li', Bi', H1i') by 100000 to form two new object coordinate sequences, namely the control point magnified coordinate sequence and the interpolated magnified coordinate sequence. Filter the gross errors of the control point magnified coordinate sequence and the interpolated magnified coordinate sequence using the RANSAC method with an error threshold δ as the limit, forming a new control point magnified coordinate sequence and an interpolated magnified coordinate sequence. Calculate the correlation coefficient θ between the elevation values ​​of the new control point magnified coordinate sequence and the interpolated magnified coordinate sequence. Transform the control point latitude and longitude coordinates (Li', Bi') into (L1i', B1i').

[0014] L1i'=Li'+(mM / 2)*ΔL

[0015] B1i'=Bi'+(nN / 2)*ΔB;

[0016] S12. For each control point coordinate sequence corresponding to the m and n combination, calculate the correlation coefficient θm,n using the method in step S11; find the m and n combination (mmax, nmax) corresponding to the maximum correlation coefficient, and modify the starting coordinate of the upper left corner of the DSM to (L10, B10).

[0017] L10 = L0 + (mmax - M / 2) * ΔL

[0018] B10 = B0 + (nmax - N / 2) * ΔB

[0019] Where m and n are the elevation values ​​corresponding to the new control point magnified coordinate sequence and the interpolated magnified coordinate sequence, respectively; m = 0, 1, ..., M-1; n = 0, 1, ..., N-1;

[0020] Using interpolated latitude and longitude as variables, the affine transformation coefficients (a, b, c, d) between two sets of elevation values ​​corresponding to (mmax, nmax) are calculated based on the least squares method. The affine transformation coefficients conform to the following elevation transformation formula.

[0021] Hi'=H1i'+a+b*Li'+c*Bi'+d*Li'*Bi';

[0022] Except for invalid grid points marked with -999, the elevation of each grid point in the DSM is calculated based on the row and column numbers and the latitude and longitude of the top left corner after translation. The transformed grid point elevation is calculated according to the elevation transformation formula, and the transformed grid point elevation replaces the original grid point elevation, thereby updating the elevation of all grid points in the entire DSM.

[0023] Preferably, step S2 specifically includes the following:

[0024] S21. Set the initial block intervals in the longitude and latitude directions of the target area to Lt0 and Bt0, and calculate the actual number of blocks Ct and Rt in the longitude and latitude directions;

[0025] Ct=(L0+(C-1)*ΔL) / Lt0

[0026] Rt=(B0+(R-1)*ΔB) / Bt0;

[0027] Calculate the actual block intervals Lt and Bt in the longitude and latitude directions based on Ct and Rt;

[0028] Lt=(L0+(C-1)*ΔL) / Ct

[0029] Bt=(B0+(R-1)*ΔB) / Rt

[0030] Wherein, the number of rows in each block is Rs = Bt / Rt+1, and the number of columns is Cs = Lt / Ct+1;

[0031] S22. Let the overlap of two adjacent blocks in the longitude and latitude directions be Lc and Bc, respectively. Then the longitude range of the block in the m1th row and n1st column is:

[0032] Minimum longitude Lmin_m1n1=L0+n1*Lt-Lc

[0033] Maximum longitude Lmax_m1n1=L0+(n1+1)*Lt+Lc

[0034] The latitude range of the block in row m1 and column n1 is:

[0035] Minimum latitude Bmin_m1n1=B0+m1*Bt-Bc

[0036] Maximum latitude Bmax_m1n1=B0+(m1+1)*Bt+Bc

[0037] Where m1 = 0, ..., Rt-1; n1 = 0, ..., Ct-1;

[0038] S23. Based on the longitude and latitude range of the block in row m1 and column n1, count the control points whose coordinates are within the range, and count the longitude range dL_m1n1 and latitude range dB_m1n1 of the control points within the range; based on the relationship between the number of control points, the longitude range, the latitude range and the corresponding threshold, determine the elevation correction of the block in row m1 and column n1 using the appropriate method.

[0039] Preferably, step S23 specifically involves:

[0040] If the number of control points is less than a threshold Nmin, or dL_m1n1 is less than a certain threshold, or dB_m1n1 is less than a certain threshold, then only the bilinear interpolation method in step S1 is used to obtain the interpolation high program sequence and the high program sequence of control points within the range. The corresponding elevations of the two high program sequences are subtracted to obtain an elevation difference sequence. All values ​​in the elevation difference sequence are sorted, and the median is taken as the elevation correction amount Hjzh_m1n1 of the m1-th row and n1-th column block. If the number of control points is less than a threshold Nmin / 2, then the elevation correction amount of the m1-th row and n1-th column block is recorded as an invalid value.

[0041] If the number of control points exceeds a threshold Nmin and both dL_m1n1 and dB_m1n1 exceed certain thresholds, then the row and column coordinates (ri', ci') of each control point within this range are calculated on the DSM based on its longitude and latitude. The corresponding elevation H1i' is obtained through bilinear interpolation of the row and column coordinates, forming a second ground coordinate sequence (Li', Bi', H1i'). An error threshold δ is set, and the longitude and latitude coordinates of the coordinate sequences (Li', Bi', Hi') and (Li', Bi', H1i') are then compared. Multiply by 100000 to form two new object coordinate sequences: a control point magnified coordinate sequence and an interpolated magnified coordinate sequence. Filter out gross errors in the control point magnified coordinate sequence and the interpolated magnified coordinate sequence using RANSAC method with δ as the limit difference, forming new control point magnified coordinate sequences and interpolated coordinate sequences. Subtract the corresponding elevations from the control point magnified coordinate sequence and the interpolated coordinate sequence to obtain an elevation difference sequence. Sort all values ​​in the elevation difference sequence and take the median as the elevation correction Hjzh_m1n1 for the m1-th row and n1-th column block.

[0042] Preferably, step S3 specifically includes the following:

[0043] S31. Let the blocks with effective elevation correction obtained in step S2 be called correction blocks, and the number of blocks be Njzh; record the elevation correction and the corresponding row and column number of each correction block to form the record information of the correction block.

[0044] S32. For the m2th correction block, set a two-dimensional connection array with a length and width of Ct*Rt+Njzh. The values ​​in the connection array are initialized to 0, indicating no connection; where m2=0,…,Njzh-1;

[0045] For the i-th row of the concatenated array, if i < Ct * Rt, then its adjacent row block number is (m1i_j, n1i_j):

[0046] m1i_0=i / Ct-1, n1i_0=i%Ct

[0047] m1i_1=i / Ct, n1i_1=i%Ct-1

[0048] m1i_2=i / Ct, n1i_2=i%Ct+1

[0049] m1i_3=i / Ct+1, n1i_3=i%Ct

[0050] Where i starts from 0, and j = 0, 1, 2, 3;

[0051] If m1i_j is not less than 0 and not greater than Rt-1, and n1i_j is not less than 0 and not greater than Ct-1, then set the value of the i-th row and m1i_j*Ct+n1i_j column of the connection array to 1, indicating that there is a connection; at the same time, search all correction information. If the m1 in the correction information record is equal to m1i_j and n1 is equal to n1i_j, then set the value of the i-th row and Ct*Rt+m2 column of the connection array to 1.

[0052] S33. Set up an array ifjzh of length Ct*Rt+Njzh to store whether the block DSM corresponding to the m3th array element in row m3%Ct needs to be corrected, where 0 means no correction is needed and 1 means correction is needed; where m3 = 0, ..., Ct*Rt+Njzh-1;

[0053] S34. Set up two elevation correction arrays lstHjzh of length Ct*Rt+Njzh to store the elevation correction of the previous iteration. The n3rd array element represents the previous elevation correction of the block DSM in row n3 / Ct and column n3%Ct. All elements of the two arrays are initialized to 0. Wherein, n3 = 0, ..., Ct*Rt+Njzh-1.

[0054] Preferably, the update strategy for each iteration is as follows:

[0055] A1. Set up two elevation correction arrays Hjzh with length Ct*Rt+Njzh. The n3rd element represents the elevation correction of block DSM in row n3 / Ct and column n3%Ct. The first Ct*Rt elements of both arrays are initialized to Hjzh_m1n1 calculated in step S2, where m1 = n3 / Ct and n1 = n3%Ct. All elements after the Ct*Rt element are initialized to 0.

[0056] A2. Set two correction count arrays jlnum with length Ct*Rt+Njzh, and initialize all elements of both arrays to 0;

[0057] A3. If the m4th row and n4th column of the concatenated array are both 1; where m4 = 0, ..., Ct*Rt-1 and n4 = m4+1, ..., Ct*Rt-1; then calculate the elevation correction:

[0058] Hjzh[m4]=Hjzh[m4]+0.5*(lstHjzh[n4]-lstHjzh[m4])

[0059] Hjzh[n4]=Hjzh[n4]+0.5*(lstHjzh[m4]-lstHjzh[n4])

[0060] jlnum[m4]=jlnum[m4]+1

[0061] jlnum[n4] = jlnum[n4] + 1

[0062] if jzh[m4] == 1

[0063] if jzh[n4] == 1

[0064] After calculating all correction amounts in a loop, save the average correction amount to the lstHjzh array:

[0065] lstjzh[m4] = Hjzh[m4] / jlnum[m4];

[0066] A4. If the above iterative update step is executed 100 times to end the iteration, the elevation correction amount of the block in the m1-th row and n1-th column is lstjzh[m1 * Ct + n1].

[0067] Preferably, step S4 specifically includes the following content

[0068] S41. Update each column corresponding to each row block; where the elevation correction amount of the j-th column of the block in the m1-th row and n1-th column is Hjzh_m1_n1_j

[0069] (1) When n1 is equal to 0:

[0070] If j < Cs / 2, then Hjzh_m1_n1_j = 0; if j ≥ Cs / 2, then Hjzh_m1_n1_j = (1 - (j - Cs / 2) / Cs) * lstjzh[m1 * Ct + n1] + (j - Cs / 2) / Cs * lstjzh[m1 * Ct + n1 + 1];

[0071] (2) When n1 = Ct - 1:

[0072] If j ≥ Cs / 2, then Hjzh_m1_n1_j = 0; if j < Cs / 2, then Hjzh_m1_n1_j = (1 - j / Cs) * lstjzh[m1 * Ct + n1] + j / Cs * lstjzh[m1 * Ct + n1 - 1];

[0073] (3) When n1 is greater than 0 and n1 < Ct - 1:

[0074] If j ≥ Cs / 2, then Hjzh_m1_n1_j = ((1 - (j - Cs / 2) / Cs) * lstjzh[m1 * Ct + n1] + (j - Cs / 2) / Cs * lstjzh[m1 * Ct + n1 + 1]); if j < Cs / 2, then Hjzh_m1_n1_j = (1 - j / Cs) * lstjzh[m1 * Ct + n1] + j / Cs * lstjzh[m1 * Ct + n1 - 1];

[0075] S42. Update each row corresponding to each column block; the elevation correction amount of the i-th row of the block in the m1-th row and n1-th column is Hjzh_m1_n1_i;

[0076] (1) When m1 is equal to 0: '

[0077] If i < Rs / 2, then Hjzh_m1_n1_i = 0; if i ≥ Rs / 2, then Hjzh_m1_n1_i = (1 - (i - Rs / 2) / Rs) * lstjzh[m1 * Ct + n1] + (i - Rs / 2) / Rs * lstjzh[(m1 + 1) * Ct + n1];

[0078] (2) When m1 = Rt - 1:

[0079] If i ≥ Rs / 2, then Hjzh_m1_n1_i = 0; if i < Rs / 2, then Hjzh_m1_n1_i = (1 - i / Rs) * lstjzh[m1 * Ct + n1] + i / Rs * lstjzh[(m1 - 1) * Ct + n1];

[0080] (3) When m1 is greater than 0 and m1 < Rt - 1:

[0081] If i ≥ Rs / 2, then Hjzh_m1_n1_i = ((1 - (i - Rs / 2) / Rs) * lstjzh[m1 * Ct + n1] + (i - Rs / 2) / Rs * lstjzh[(m1 + 1) * Ct + n1]; if i < Rs / 2, then Hjzh_m1_n1_i = (1 - i / Rs) * lstjzh[m1 * Ct + n1] + i / Rs * lstjzh[(m1 - 1) * Ct + n1].

[0082] The beneficial effects of this invention are: 1. The method of this invention performs overall terrain correction on the DSM using all control points, improving the accuracy of horizontal and vertical alignment. 2. The method of this invention divides the large-area DSM into blocks, and each block is individually corrected according to the distribution of control points, obtaining correction information with higher reliability for some blocks. 3. The method of this invention uses iterative correction to propagate the correction information globally through adjacent blocks, achieving overall correction of all blocks. 4. The method of this invention obtains the correction information for each grid point of each block by interpolating the inverse distances in both the row and column directions. Attached Figure Description

[0083] Figure 1 This is a flowchart of the method in an embodiment of the present invention. Detailed Implementation

[0084] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0085] like Figure 1 As shown in this embodiment, a method for improving the accuracy of a large-area DSM using control points is provided. This method, based on overall correction, employs a strategy combining block processing and iterative correction, which can improve correction accuracy while ensuring that the corrected large-area DSM does not have edge marks caused by block division. In addition to correcting the stitched large-area DSM, this method can also correct large-area multi-frame single-scene DSMs before stitching, treating the single-scene DSM as a block DSM according to this invention. The method specifically includes the following parts:

[0086] I. Overall Topographic Correction

[0087] Bilinear interpolation is used to obtain the elevation corresponding to the row and column coordinates of each control point in the target area on the DSM, thus obtaining a second ground coordinate sequence. The control point magnified coordinate sequence and interpolated magnified coordinate sequence obtained based on the ground 3D coordinates in the control point coordinate file and the second ground coordinate sequence are respectively used to obtain a new control point magnified coordinate sequence and interpolated magnified coordinate sequence using the RANSAC method, and the correlation coefficient between the corresponding elevation values ​​of the two new coordinate sequences is calculated. The starting coordinate of the upper left corner of the DSM is modified based on the elevation value combination with the largest correlation coefficient. The affine transformation coefficient between the elevation values ​​of the corresponding two sets of coordinate sequences is calculated based on the least squares method, and the original grid point elevation is replaced with the transformed grid point elevation calculated based on the affine transformation coefficient to update the elevation of all grid points in the entire DSM.

[0088] The input data consists of a raster DSM file to be calibrated and a control point coordinate file. The raster DSM file contains the number of rows R, the number of columns C, the row resolution ΔB, the column resolution ΔL, and the coordinates from the top-left corner (L0, B0), where L0 is the longitude of the top-left raster point and B0 is the latitude of the top-left raster point. The control point coordinate file contains N ground 3D coordinates (Li', Bi', Hi') (Li', Bi', and Hi' are the longitude, latitude, and elevation of the i'th point, respectively, i' = 0, ..., N-1). The search window size is set to M*N, where M and N are odd numbers, and the search step size is ΔS in degrees.

[0089] 1.1 Calculate the row and column coordinates (ri', ci') of each control point on the DSM based on its longitude and latitude. Obtain the corresponding elevation H1i' using bilinear interpolation of the row and column coordinates. This forms a second ground coordinate sequence (Li', Bi', H1i') (i' = 0, ..., N-1). Set an error threshold δ. Multiply the longitude and latitude coordinates of the coordinate sequences (Li', Bi', Hi') and (Li', Bi', H1i') by 100000 to form two new object coordinate sequences (referred to as the control point magnified coordinate sequence and the interpolated magnified coordinate sequence, respectively). Filter out outliers in the control point magnified coordinate sequence and the interpolated magnified coordinate sequence using the RANSAC method with δ as the limit of error, forming new control point magnified coordinate sequences and interpolated magnified coordinate sequences. Calculate the correlation coefficient θ between the elevation values ​​of the two coordinate sequences.

[0090] 1.2. Transform the control point's latitude and longitude coordinates (Li', Bi') into (L1i', B1i'):

[0091] L1i'=Li'+(mM / 2)*ΔL

[0092] B1i'=Bi'+(nN / 2)*ΔB;

[0093] For each control point coordinate sequence corresponding to a combination of m and n, calculate a correlation coefficient θ using the method described above. m，n (m and n are the elevation values ​​corresponding to the new control point magnified coordinate sequence and the interpolated magnified coordinate sequence, respectively; m = 0, 1, ..., M-1; n = 0, 1, ..., N-1).

[0094] Find the m and n combination corresponding to the maximum correlation coefficient, let's say it's (mmax, nmax), and modify the starting coordinates of the top left corner of the DSM to (L10, B10):

[0095] L10 = L0 + (mmax - M / 2) * ΔL

[0096] B10 = B0 + (nmax - N / 2) * ΔB

[0097] Using interpolated latitude and longitude as variables, the affine transformation coefficients (a, b, c, d) between the elevation values ​​of the two corresponding elevation sequence (mmax, nmax) are calculated based on the least squares method. The affine transformation coefficients conform to the following elevation transformation formula:

[0098] Hi'=H1i'+a+b*Li'+c*Bi'+d*Li'*Bi'

[0099] Except for invalid grid points marked with -999, the elevation of each grid point in the DSM is calculated based on the row and column numbers and the latitude and longitude of the top left corner after translation. The transformed grid point elevation is calculated according to the elevation transformation formula, and the transformed grid point elevation replaces the original grid point elevation, thereby updating the elevation of all grid points in the entire DSM.

[0100] II. Calculation of Correction Parameters

[0101] The target area is divided into regular grids, and the control points within each grid are counted. Based on the number of control points within each grid and the relationship between the latitude and longitude range of the control points and the corresponding thresholds, the elevation correction amount for each grid is determined.

[0102] 2.1. Set the initial block intervals in the longitude and latitude directions to Lt0 and Bt0, respectively, and calculate the actual number of blocks Ct and Rt in the longitude and latitude directions:

[0103] Ct=(L0+(C-1)*ΔL) / Lt0

[0104] Rt=(B0+(R-1)*ΔB) / Bt0

[0105] Calculate the actual block intervals Lt and Bt in the longitude and latitude directions based on Ct and Rt:

[0106] Lt=(L0+(C-1)*ΔL) / Ct

[0107] Bt=(B0+(R-1)*ΔB) / Rt

[0108] In this case, the number of rows in each block is Rs = Bt / Rt+1, and the number of columns is Cs = Lt / Ct+1.

[0109] 2.2. Let the overlap of two adjacent blocks in the longitude and latitude directions be Lc and Bc, respectively. Then the longitude range of the block in the m1th row and n1th column (m1 = 0, ..., Rt-1; n1 = 0, ..., Ct-1) is:

[0110] Minimum longitude Lmin_m1n1=L0+n1*Lt-Lc

[0111] Maximum longitude Lmax_m1n1=L0+(n1+1)*Lt+Lc

[0112] The latitude range of the block in row m1 and column n1 is:

[0113] Minimum latitude Bmin_m1n1=B0+m1*Bt-Bc

[0114] Maximum latitude Bmax_m1n1 = B0 + (m1 + 1) * Bt + Bc;

[0115] 2.3. Based on the longitude and latitude range of the m1th row block and the n1th column block, count the control points whose coordinates are within the range, and count the longitude range dL_m1n1 and latitude range dB_m1n1 of the control points within the range.

[0116] (1) If the number of control points is less than a threshold Nmin, or dL_m1n1 is less than a certain threshold, or dB_m1n1 is less than a certain threshold, then only the bilinear interpolation method in step S1 is used to obtain the interpolation high program sequence and the high program sequence of control points within the range. Then, the corresponding elevations of the two high program sequences are subtracted to obtain an elevation difference sequence. All values ​​in the elevation difference sequence are sorted, and the median is taken as the elevation correction amount Hjzh_m1n1 of the m1-th row and n1-th column block. If the number of control points is less than a threshold Nmin / 2, then the elevation correction amount of the m1-th row and n1-th column block is recorded as an invalid value.

[0117] (2) If the number of control points is greater than a threshold Nmin and both dL_m1n1 and dB_m1n1 are greater than a certain threshold, then calculate the row and column coordinates (ri', ci') of each control point on the DSM based on its longitude and latitude. Obtain the corresponding elevation H1i' based on the bilinear interpolation of the row and column coordinates to form a second ground coordinate sequence (Li', Bi', H1i'). Set an error threshold δ and then convert the longitude and latitude of the coordinate sequences (Li', Bi', Hi') and (Li', Bi', H1i') into a single coordinate sequence. Multiply the degree coordinates by 100000 to form two new object coordinate sequences: the control point magnified coordinate sequence and the interpolated magnified coordinate sequence. Filter out gross errors in the control point magnified coordinate sequence and the interpolated magnified coordinate sequence using the RANSAC method with δ as the limit difference, forming new control point magnified coordinate sequences and interpolated coordinate sequences. Subtract the corresponding elevations from the control point magnified coordinate sequence and the interpolated coordinate sequence to obtain an elevation difference sequence. Sort all values ​​in the elevation difference sequence and take the median as the elevation correction Hjzh_m1n1 for the m1-th row and n1-th column block.

[0118] III. Correction Parameter Update

[0119] The elevation correction values ​​of the obtained grids with effective elevation correction values ​​are iteratively updated to obtain the corresponding grids after iterative updates.

[0120] 3.1 Assume that the number of blocks for obtaining effective elevation corrections in the second step is Njzh, hereinafter referred to as correction blocks. For each record, the elevation correction and the corresponding row and column number of each correction block are recorded. For example, the record information of the m2th (m2=0,…,Njzh-1) correction block is Inf_m2(Hjzh_m1n1,m1,n1). Set up a two-dimensional connection array (hereinafter referred to as connection array) with a length and width of Ct*Rt+Njzh. The values ​​in the array are initialized to 0, indicating no connection.

[0121] 3.2 For the i-th row of the concatenated array (i starts from 0), if i < Ct*Rt, then its adjacent row block number is (m1i_j, n1i_j), j = 0, 1, 2, 3:

[0122] m1i_0=i / Ct-1, n1i_0=i%Ct

[0123] m1i_1=i / Ct, n1i_1=i%Ct-1

[0124] m1i_2=i / Ct, n1i_2=i%Ct+1

[0125] m1i_3=i / Ct+1, n1i_3=i%Ct

[0126] If m1i_j is not less than 0 and not greater than Rt-1, and n1i_j is not less than 0 and not greater than Ct-1, then set the value of the i-th row and m1i_j*Ct+n1i_j column of the connection array to 1, indicating that there is a connection. Simultaneously, retrieve all correction information Inf_m2. If m1 equals m1i_j and n1 equals n1i_j in the information record, then set the value of the i-th row and Ct*Rt+m2 column of the connection array to 1.

[0127] 3.3. Set up an array ifjzh of length Ct*Rt+Njzh to store whether the m3th array element (m3=0,…,Ct*Rt+Njzh-1) corresponds to the m3 / Ct row block and m3%Ct column block DSM need to be corrected. 0 means no correction is needed and 1 means correction is needed.

[0128] Set up two elevation correction arrays lstHjzh of length Ct*Rt+Njzh to store the elevation correction from the previous iteration. The n3rd element (n3 = 0, ..., Ct*Rt+Njzh-1) represents the previous elevation correction for the n3 / Ct row block and n3%Ct column block DSM. All elements in both arrays are initialized to 0.

[0129] In this embodiment, the update strategy for each iteration is as follows:

[0130] A1. Set up two elevation correction arrays Hjzh of length Ct*Rt+Njzh. The n3rd element (n3=0,…,Ct*Rt+Njzh-1) represents the elevation correction of the n3 / Ct row block and n3%Ct column block DSM. The first Ct*Rt elements of both arrays are initialized to Hjzh_m1n1 calculated in step S2, where m1=n3 / Ct and n1=n3%Ct. All elements after the Ct*Rt element are initialized to 0.

[0131] A2. Set two correction count arrays jlnum with length Ct*Rt+Njzh, and initialize all elements of both arrays to 0.

[0132] A3. If the m4th row (m4 = 0, ..., Ct*Rt-1) and the n4th column (n4 = m4+1, ..., Ct*Rt-1) of the concatenated array are both 1, then calculate the correction:

[0133] Hjzh[m4]=Hjzh[m4]+0.5*(lstHjzh[n4]-lstHjzh[m4]);

[0134] Hjzh[n4]=Hjzh[n4]+0.5*(lstHjzh[m4]-lstHjzh[n4]);

[0135] jlnum[m4] = jlnum[m4] + 1;

[0136] jlnum[n4] = jlnum[n4] + 1;

[0137] ifjzh[m4] = 1;

[0138] ifjzh[n4] = 1;

[0139] After iterating through all corrections, save the average correction to the lstHjzh array:

[0140] lstjzh[m4]=Hjzh[m4] / jlnum[m4];

[0141] A4. The above iterative update steps are executed 100 times to end the iteration. Then the elevation correction of the block in row m1 and column n1 (m1=0,…,Rt-1;n1=0,…,Ct-1) is lstjzh[m1*Ct+n1].

[0142] IV. Complete DSM Update

[0143] Based on the elevation correction amount of the corresponding grid after iterative update, taking the center of each grid as the reference, bilinear interpolation is used to calculate the elevation correction amount of each row and column of the grid to update the entire DSM.

[0144] 4.1. Update each column corresponding to each row block; among them, the elevation correction amount of the j-th column in the block of the m1-th row and n1-th column is Hjzh_m1_n1_j;

[0145] (1) When n1 is equal to 0:

[0146] If j < Cs / 2, then Hjzh_m1_n1_j = 0; if j ≥ Cs / 2, then Hjzh_m1_n1_j = (1 - (j - Cs / 2) / Cs) * lstjzh[m1 * Ct + n1] + (j - Cs / 2) / Cs * lstjzh[m1 * Ct + n1 + 1];

[0147] (2) When n1 = Ct - 1:

[0148] If j ≥ Cs / 2, then Hjzh_m1_n1_j = 0; if j < Cs / 2, then Hjzh_m1_n1_j = (1 - j / Cs) * lstjzh[m1 * Ct + n1] + j / Cs * lstjzh[m1 * Ct + n1 - 1];

[0149] (3) When n1 is greater than 0 and n1 < Ct - 1:

[0150] If j ≥ Cs / 2, then Hjzh_m1_n1_j = ((1 - (j - Cs / 2) / Cs) * lstjzh[m1 * Ct + n1] + (j - Cs / 2) / Cs * lstjzh[m1 * Ct + n1 + 1]; if j < Cs / 2, then Hjzh_m1_n1_j = (1 - j / Cs) * lstjzh[m1 * Ct + n1] + j / Cs * lstjzh[m1 * Ct + n1 - 1];

[0151] 4.2. Update each row corresponding to each column block; the elevation correction amount of the i-th row in the block of the m1-th row and n1-th column is Hjzh_m1_n1_i;

[0152] (1) When m1 is equal to 0:

[0153] If i < Rs / 2, then Hjzh_m1_n1_i = 0; if i ≥ Rs / 2, then Hjzh_m1_n1_i = (1 - (i - Rs / 2) / Rs) * lstjzh[m1 * Ct + n1] + (i - Rs / 2) / Rs * lstjzh[(m1 + 1) * Ct + n1];​​​

[0155] If i ≥ Rs / 2, then Hjzh_m1_n1_i = 0; if i < Rs / 2, then Hjzh_m1_n1_i = (1-i / Rs)*lstjzh[m1*Ct+n1] + i / Rs*lstjzh[(m1-1)*Ct+n1];

[0156] (3) When m1 is greater than 0 and m1 < Rt-1:

[0157] If i≥Rs / 2, then Hjzh_m1_n1_i=((1-(i-Rs / 2) / Rs)*lstjzh[m1*Ct+n1]+(i-Rs / 2) / Rs*lstjzh[(m1+1) *Ct+n1]; if i<Rs / 2, then Hjzh_m1_n1_i=(1-i / Rs)*lstjzh[m1*Ct+n1]+i / Rs*lstjzh[(m1-1)*Ct+n1].

[0158] By adopting the above-disclosed technical solution of this invention, the following beneficial effects are obtained:

[0159] This invention provides a method for improving the accuracy of a large-area topographic surveying and mapping (DSM) using control points. The method performs overall terrain correction on the DSM using all control points, improving both horizontal and vertical accuracy. The method divides the large-area DSM into blocks, and each block is individually corrected based on the distribution of control points, obtaining correction information with higher reliability from some blocks. The method iteratively propagates the correction information globally through adjacent blocks, achieving overall correction for all blocks. Finally, the method uses inverse distance interpolation in both row and column directions to obtain the correction information for each grid point in each block.

[0160] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method for improving the accuracy of a large-area DSM using control points, characterized in that: Includes the following steps, S1. Overall Terrain Correction: Bilinear interpolation is used to obtain the elevation corresponding to the row and column coordinates of each control point in the target area on the DSM, thus obtaining a second ground coordinate sequence. The control point magnified coordinate sequence and interpolated magnified coordinate sequence obtained based on the ground 3D coordinates in the control point coordinate file and the second ground coordinate sequence are respectively used to obtain a new control point magnified coordinate sequence and interpolated magnified coordinate sequence using the RANSAC method, and the correlation coefficient between the corresponding elevation values ​​of the two new coordinate sequences is calculated. The starting coordinate of the upper left corner of the DSM is modified based on the elevation value combination with the largest correlation coefficient. The affine transformation coefficient between the elevation values ​​of the corresponding two sets of coordinate sequences is calculated based on the least squares method, and the original grid point elevation is replaced with the transformed grid point elevation calculated based on the affine transformation coefficient to update the elevation of all grid points in the entire DSM. S2. Correction parameter calculation: Divide the target area into regular grids, count the control points within each grid, and determine the elevation correction amount for each grid based on the number of control points within each grid and the relationship between the latitude and longitude range of the control points and the corresponding thresholds. S3. Correction parameter update: Iteratively update the elevation correction of the grid with effective elevation correction obtained in step S2 to obtain the corresponding grid after iterative update. S4. Full DSM Update: Based on the elevation correction of the corresponding grid after iterative update, the elevation correction of each row and column of the grid is calculated by bilinear interpolation with the center of each grid as the reference, so as to update the entire DSM.

2. The method for improving the accuracy of a large-area DSM using control points according to claim 1, characterized in that: In step S1, the raster DSM file to be corrected and the control point coordinate file are used as input data. The raster DSM file includes the number of rows R, the number of columns C, the row direction resolution ΔB, the column direction resolution ΔL, and the coordinates of the top left corner (L0, B0); where L0 is the longitude of the top left corner raster point and B0 is the latitude of the top left corner raster point. The control point coordinate file contains N ground three-dimensional coordinates (Li', Bi', Hi'); where Li', Bi', and Hi' are the longitude, latitude, and elevation of the i'th control point, respectively; i' = 0, 1, ..., N-1.

3. The method for improving the accuracy of a large-area DSM using control points according to claim 2, characterized in that: Step S1 specifically includes the following: S11. Set the search window to M*N and the search step size to ΔS. Calculate the row and column coordinates (ri', ci') of each control point on the DSM based on its latitude and longitude. Obtain the corresponding elevation H1i based on the bilinear interpolation of the row and column coordinates, and obtain the second ground coordinate sequence (Li', Bi', H1i'). Multiply the latitude and longitude coordinates of the coordinate sequences (Li', Bi', Hi') and (Li', Bi', H1i') by 100000 to form two new object coordinate sequences, namely the control point magnified coordinate sequence and the interpolated magnified coordinate sequence. Filter the gross errors of the control point magnified coordinate sequence and the interpolated magnified coordinate sequence using the RANSAC method with an error threshold δ as the limit, forming a new control point magnified coordinate sequence and an interpolated magnified coordinate sequence. Calculate the correlation coefficient θ between the elevation values ​​of the new control point magnified coordinate sequence and the interpolated magnified coordinate sequence. Transform the control point latitude and longitude coordinates (Li', Bi') into (L1i', B1i'). L1i'=Li'+(mM / 2)*ΔL B1i'=Bi'+(nN / 2)*ΔB; S12. For each control point coordinate sequence corresponding to the m and n combination, calculate the correlation coefficient θm,n using the method in step S11; find the m and n combination (mmax, nmax) corresponding to the maximum correlation coefficient, and modify the starting coordinate of the upper left corner of the DSM to (L10, B10). L10 = L0 + (mmax - M / 2) * ΔL B10 = B0 + (nmax - N / 2) * ΔB Where m and n are the elevation values ​​corresponding to the new control point magnified coordinate sequence and the interpolated magnified coordinate sequence, respectively; m = 0, 1, ..., M-1; n = 0, 1, ..., N-1; Using interpolated latitude and longitude as variables, the affine transformation coefficients (a, b, c, d) between two sets of elevation values ​​corresponding to (mmax, nmax) are calculated based on the least squares method. The affine transformation coefficients conform to the following elevation transformation formula. Hi'=H1i'+a+b*Li'+c*Bi'+d*Li'*Bi'; Except for invalid grid points marked with -999, the elevation of each grid point in the DSM is calculated based on the row and column numbers and the latitude and longitude of the top left corner after translation. The transformed grid point elevation is calculated according to the elevation transformation formula, and the transformed grid point elevation replaces the original grid point elevation, thereby updating the elevation of all grid points in the entire DSM.

4. The method for improving the accuracy of a large-area DSM using control points according to claim 3, characterized in that: Step S2 specifically includes the following: S21. Set the initial block intervals in the longitude and latitude directions of the target area to Lt0 and Bt0, and calculate the actual number of blocks Ct and Rt in the longitude and latitude directions; Ct=(L0+(C-1)*ΔL) / Lt0 Rt=(B0+(R-1)*ΔB) / Bt0; Calculate the actual block intervals Lt and Bt in the longitude and latitude directions based on Ct and Rt; Lt=(L0+(C-1)*ΔL) / Ct Bt=(B0+(R-1)*ΔB) / Rt Wherein, the number of rows in each block is Rs = Bt / Rt+1, and the number of columns is Cs = Lt / Ct+1; S22. Let the overlap of two adjacent blocks in the longitude and latitude directions be Lc and Bc; then the longitude range of the block in the m1th row and n1st column is... Minimum longitude Lmin_m1n1=L0+n1*Lt-Lc Maximum longitude Lmax_m1n1=L0+(n1+1)*Lt+Lc The latitude range of the block in row m1 and column n1 is: Minimum latitude Bmin_m1n1=B0+m1*Bt-Bc Maximum latitude Bmax_m1n1=B0+(m1+1)*Bt+Bc Where m1 = 0, ..., Rt-1; n1 = 0, ..., Ct-1; S23. Based on the longitude and latitude range of the block in row m1 and column n1, count the control points whose coordinates are within the range, and count the longitude range dL_m1n1 and latitude range dB_m1n1 of the control points within the range; based on the relationship between the number of control points, the longitude range, the latitude range and the corresponding threshold, determine the elevation correction of the block in row m1 and column n1 using the appropriate method.

5. The method for improving the accuracy of a large-area DSM using control points according to claim 4, characterized in that: Step S23 specifically involves, If the number of control points is less than a threshold Nmin, or dL_m1n1 is less than a certain threshold, or dB_m1n1 is less than a certain threshold, then only the bilinear interpolation method in step S1 is used to obtain the interpolation high program sequence and the high program sequence of control points within the range. The corresponding elevations of the two high program sequences are subtracted to obtain an elevation difference sequence. All values ​​in the elevation difference sequence are sorted, and the median is taken as the elevation correction amount Hjzh_m1n1 of the m1-th row and n1-th column block. If the number of control points is less than a threshold Nmin / 2, then the elevation correction amount of the m1-th row and n1-th column block is recorded as an invalid value. If the number of control points exceeds a threshold Nmin and both dL_m1n1 and dB_m1n1 exceed certain thresholds, then the row and column coordinates (ri', ci') of each control point within this range are calculated on the DSM based on its longitude and latitude. The corresponding elevation H1i' is obtained through bilinear interpolation of the row and column coordinates, forming a second ground coordinate sequence (Li', Bi', H1i'). An error threshold δ is set, and the longitude and latitude coordinates of the coordinate sequences (Li', Bi', Hi') and (Li', Bi', H1i') are then compared. Multiply by 100000 to form two new object coordinate sequences: a control point magnified coordinate sequence and an interpolated magnified coordinate sequence. Filter out gross errors in the control point magnified coordinate sequence and the interpolated magnified coordinate sequence using RANSAC method with δ as the limit difference, forming new control point magnified coordinate sequences and interpolated coordinate sequences. Subtract the corresponding elevations from the control point magnified coordinate sequence and the interpolated coordinate sequence to obtain an elevation difference sequence. Sort all values ​​in the elevation difference sequence and take the median as the elevation correction Hjzh_m1n1 for the m1-th row and n1-th column block.

6. The method for improving the accuracy of a large-area DSM using control points according to claim 5, characterized in that: Step S3 specifically includes the following: S31. Let the blocks with effective elevation correction obtained in step S2 be called correction blocks, and the number of blocks be Njzh; record the elevation correction and the corresponding row and column number of each correction block to form the record information of the correction block. S32. For the m2th correction block, set a two-dimensional connection array with a length and width of Ct*Rt+Njzh. The values ​​in the connection array are initialized to 0, indicating no connection; where m2=0,…,Njzh-1; For the i-th row of the concatenated array, if i < Ct * Rt, then its adjacent row block number is (m1i_j, n1i_j): m1i_0=i / Ct-1, n1i_0=i%Ct m1i_1=i / Ct, n1i_1=i%Ct-1 m1i_2=i / Ct, n1i_2=i%Ct+1 m1i_3=i / Ct+1, n1i_3=i%Ct Where i starts from 0, and j = 0, 1, 2, 3; If m1i_j is not less than 0 and not greater than Rt-1, and n1i_j is not less than 0 and not greater than Ct-1, then set the value of the i-th row and m1i_j*Ct+n1i_j column of the connection array to 1, indicating that there is a connection; at the same time, search all correction information. If the m1 in the correction information record is equal to m1i_j and n1 is equal to n1i_j, then set the value of the i-th row and Ct*Rt+m2 column of the connection array to 1. S33. Set up an array ifjzh of length Ct*Rt+Njzh to store whether the block DSM corresponding to the m3th array element in row m3%Ct needs to be corrected, where 0 means no correction is needed and 1 means correction is needed; where m3 = 0, ..., Ct*Rt+Njzh-1; S34. Set up two elevation correction arrays lstHjzh of length Ct*Rt+Njzh to store the elevation correction of the previous iteration. The n3rd array element represents the previous elevation correction of the block DSM in row n3 / Ct and column n3%Ct. All elements of the two arrays are initialized to 0. Wherein, n3 = 0, ..., Ct*Rt+Njzh-1.

7. The method for improving the accuracy of a large-area DSM using control points according to claim 6, characterized in that: The update strategy for each iteration is as follows: A1. Set up two elevation correction arrays Hjzh with length Ct*Rt+Njzh. The n3rd element represents the elevation correction of block DSM in row n3 / Ct and column n3%Ct. The first Ct*Rt elements of both arrays are initialized to Hjzh_m1n1 calculated in step S2, where m1 = n3 / Ct and n1 = n3%Ct. All elements after the Ct*Rt element are initialized to 0. A2. Set two arrays jlnum of correction counts with length Ct*Rt+Njzh, and initialize all elements of both arrays to 0; A3. If the m4th row and n4th column of the concatenated array are both 1; where m4 = 0, ..., Ct*Rt-1 and n4 = m4+1, ..., Ct*Rt-1; then calculate the elevation correction: Hjzh[m4]=Hjzh[m4]+0.5*(lstHjzh[n4]-lstHjzh[m4]) Hjzh[n4] = Hjzh[n4] + 0.5 * (lstHjzh[m4] - lstHjzh[n4]) jlnum[m4] = jlnum[m4] + 1 jlnum[n4] = jlnum[n4] + 1 if jzh[m4] == 1 if jzh[n4] == 1 After calculating all correction amounts in a loop, save the average correction amount to the lstHjzh array: lstjzh[m4] = Hjzh[m4] / jlnum[m4] A4. If the above iterative update steps are executed 100 times to end the iteration, the elevation correction amount of the block in the m1-th row and n1-th column is lstjzh[m1 * Ct + n1].

8. The method for improving the accuracy of a large-area DSM using control points according to claim 7, characterized in that: Step S4 specifically includes the following content, S41. Update each column corresponding to each row block; among them, the elevation correction amount of the j-th column of the block in the m1-th row and n1-th column is Hjzh_m1_n1_j; (1) When n1 is equal to 0: If j < Cs / 2, then Hjzh_m1_n1_j = 0; if j ≥ Cs / 2, then Hjzh_m1_n1_j = (1 - (j - Cs / 2) / Cs) * lstjzh[m1 * Ct + n1] + (j - Cs / 2) / Cs * lstjzh[m1 * Ct + n1 + 1]; (2) When n1 = Ct - 1: If j ≥ Cs / 2, then Hjzh_m1_n1_j = 0; if j < Cs / 2, then Hjzh_m1_n1_j = (1 - j / Cs) * lstjzh[m1 * Ct + n1] + j / Cs * lstjzh[m1 * Ct + n1 - 1]; (3) When n1 is greater than 0 and n1 < Ct - 1: If j ≥ Cs / 2, then Hjzh_m1_n1_j = ((1 - (j - Cs / 2) / Cs) * lstjzh[m1 * Ct + n1] + (j - Cs / 2) / Cs * lstjzh[m1 * Ct + n1 + 1]; if j < Cs / 2, then Hjzh_m1_n1_j = (1 - j / Cs) * lstjzh[m1 * Ct + n1] + j / Cs * lstjzh[m1 * Ct + n1 - 1]; S42. Update each row corresponding to each column block; the elevation correction amount of the i-th row of the block in the m1-th row and n1-th column is Hjzh_m1_n1_i; (1) When m1 is equal to 0: If i < Rs / 2, then Hjzh_m1_n1_i = 0; if i ≥ Rs / 2, then Hjzh_m1_n1_i = (1 - (i - Rs / 2) / Rs) * lstjzh[m1 * Ct + n1] + (i - Rs / 2) / Rs * lstjzh[(m1 + 1) * Ct + n1]; (2) When m1 = Rt - 1: If i ≥ Rs / 2, then Hjzh_m1_n1_i = 0; if i < Rs / 2, then Hjzh_m1_n1_i = (1 - i / Rs) * lstjzh[m1 * Ct + n1] + i / Rs * lstjzh[(m1 - 1) * Ct + n1]; (3) When m1 is greater than 0 and m1 < Rt-1: If i≥Rs / 2, then Hjzh_m1_n1_i=((1-(i-Rs / 2) / Rs)*lstjzh[m1*Ct+n1]+(i-Rs / 2) / Rs*lstjzh[(m1+1) *Ct+n1]; if i<Rs / 2, then Hjzh_m1_n1_i=(1-i / Rs)*lstjzh[m1*Ct+n1]+i / Rs*lstjzh[(m1-1)*Ct+n1].

Images