A method for improving the accuracy of large area DSM using control points

Abstract

The application discloses a method for improving the precision of a large-area DSM by using control points, which comprises overall terrain correction, overall improvement of plane and elevation precision; dividing a target area into regular grids, counting control points in each grid, determining the elevation correction amount of each grid according to the relationship between the number of control points in each grid and the corresponding threshold and the longitude and latitude range of the control points; iteratively updating the elevation correction amount of the obtained grid with the effective elevation correction amount, obtaining the iteratively updated elevation correction amount of the corresponding grid; and based on the iteratively updated elevation correction amount of the corresponding grid, taking the center of each grid as a reference, bilinearly interpolating to calculate the elevation correction amount of each row and column of the grid, and updating the whole DSM. The method has the advantages that: on the basis of overall correction, the method combines the strategies of block processing and iterative correction, can improve the correction precision, and can ensure that the large-area DSM after correction does not have the edge traces caused by block division.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of photogrammetry and remote sensing, and particularly relates to a method for improving the precision of large-area DSM by using control points. BACKGROUND

[0002] Various error transmission factors in the production process of large-area DSM often cause local or global distortion, deformation, undulation and other forms of errors in DSM. These errors can be corrected and eliminated by using control points, laser points, reference terrain and other data. Due to the differences in the type, density and reliability of control points, it is difficult to achieve ideal correction effect by using a unified correction globally. SUMMARY

[0003] The present application aims to provide a method for improving the precision of large-area DSM by using control points, thereby solving the aforementioned problems in the prior art.

[0004] To achieve the above-mentioned purpose, the technical solution adopted by the present application is as follows:

[0005] A method for improving the precision of large-area DSM by using control points, comprising the following steps,

[0006] S1, overall terrain correction: using bilinear interpolation to obtain the elevation corresponding to the row and column coordinates of each control point on the DSM in the target area to obtain a second ground coordinate sequence; using the RANSAC method to obtain a new control point magnification coordinate sequence and an interpolated magnification coordinate sequence based on the ground three-dimensional coordinates in the control point coordinate file and the second ground coordinate sequence, and calculating the correlation coefficient between the corresponding elevation values of the two new coordinate sequences; modifying the upper left corner coordinate of the DSM based on the elevation value with the largest correlation coefficient, calculating the affine transformation coefficients between the elevation values of the two coordinate sequences based on the least square method, and replacing the original grid point elevation with the transformed grid point elevation calculated based on the affine transformation coefficients to update the elevation of all grid points of the entire DSM;

[0007] S2, correction parameter calculation: dividing the target area into regular grids, counting the control points inside each grid, and determining the elevation correction amount of each grid according to the relationship between the number of control points inside each grid and the latitude and longitude range of the control points and the corresponding threshold;

[0008] S3, correction parameter update: iteratively updating the elevation correction amount of the grid with effective elevation correction amount obtained in step S2 to obtain the iteratively updated elevation correction amount of the corresponding grid;

[0009] S4, DSM update: based on the elevation correction of the corresponding grid after iterative updating, taking each grid center as the reference, the elevation correction of each row and column of the grid is calculated by bilinear interpolation to update the whole DSM.

[0010] Preferably, in step S1, the raster DSM file to be corrected and the control point coordinate file are input data, the raster DSM file contains the number of rows R, the number of columns C, the row direction resolution ΔB, the column direction resolution ΔL and the upper left corner coordinate (L0, B0); wherein L0 is the longitude of the upper left corner grid point, and B0 is the latitude of the upper left corner grid point;

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

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

[0013] S11, set the search window to M*N, and the search step to ΔS; calculate the row and column coordinates (ri', ci') of each control point on the DSM according to its longitude and latitude, and obtain the height H1i corresponding to the row and column coordinates by bilinear interpolation, and obtain the second ground coordinate sequence (Li', Bi', H1i'); multiply the longitude and latitude coordinates of the coordinate sequences (Li', Bi', Hi') and (Li', Bi', H1i') by 100000 to form two new object coordinates sequences, which are the control point enlarged coordinate sequence and the interpolation enlarged coordinate sequence respectively; filter the gross error points based on the RANSAC method with the error threshold δ as the limit difference, form the new control point enlarged coordinate sequence and the interpolation enlarged coordinate sequence, and calculate the correlation coefficient θ between the height values of the new control point enlarged coordinate sequence and the interpolation enlarged coordinate sequence; perform coordinate translation on the control point longitude and latitude coordinates (Li', Bi') to transform them into (L1i', B1i');

[0014] L1i'=Li'+(m-M / 2)*ΔL

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

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

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

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

[0019] wherein m, n are elevation values corresponding to new control point magnified coordinate sequence and interpolation magnified coordinate sequence respectively; m = 0, 1, …, M-1; n = 0, 1, …, N-1;

[0020] Taking the interpolation longitude and latitude as variables, the affine transformation coefficients (a, b, c, d) between the two groups of elevation sequences corresponding to (mmax, nmax) are calculated based on the least square method, and 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 the invalid grid points marked with -999, the elevation of each grid point of the DSM is calculated according to the row number and column number and the longitude and latitude of the upper left corner after translation, the longitude and latitude of the current grid point are calculated, the elevation of the transformed grid point is calculated according to the elevation transformation formula, the original grid point elevation is replaced by the transformed grid point elevation, and all the grid point elevations of the entire DSM are updated.

[0023] Preferably, step S2 specifically comprises the following contents,

[0024] S21, setting the initial block interval of longitude and latitude direction of the target area as Lt0 and Bt0, and calculating the actual block number Ct and Rt of longitude and latitude direction;

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

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

[0027] According to Ct and Rt, the actual block interval Lt and Bt of longitude and latitude direction are calculated;

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

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

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

[0031] S22, setting the overlap of adjacent two blocks in longitude and latitude direction as Lc and Bc, then the longitude range of the m1th row and n1th column block is,

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

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

[0034] The latitude range of the block in the m1th row and the n1th column is,

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

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

[0037] wherein m1=0,…,Rt-1; n1=0,…,Ct-1;

[0038] S23, according to the longitude and latitude range of the block in the m1th row and the n1th column, statistics control points located in the range, and statistics the longitude range dL_m1n1 and the latitude range dB_m1n1 of the control points in the range; according to the size relationship between the control point number, the longitude range and the latitude range and the corresponding threshold, a corresponding method is used to determine the height correction amount of the block in the m1th row and the n1th column.

[0039] Preferably, step S23 is specifically,

[0040] If the control point number 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 through the method of step S1 to obtain the interpolation height sequence and the height sequence of the control points in the range by bilinear interpolation, subtract the corresponding height of the two height sequences to obtain a height difference sequence, sort all values of the height difference sequence, and take the median as the height correction amount Hjzh_m1n1 of the block in the m1th row and the n1th column; if the control point number is less than a threshold Nmin / 2, then record the height correction amount of the block in the m1th row and the n1th column as an invalid value;

[0041] 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 the row and column coordinates (r'i', c'i') of the point on the DSM are calculated according to the longitude and latitude of each control point in the range, the elevation H1i' corresponding to the row and column coordinates is obtained by bilinear interpolation according to the row and column coordinates, a second ground coordinate sequence (Li', Bi', H1i') is formed, an error threshold δ is set, the longitude and latitude coordinates of the coordinate sequences (Li', Bi', Hi') and (Li', Bi', H1i') are multiplied by 100000 to form two new object coordinates sequences, which are the control point enlarged coordinate sequence and the interpolation enlarged coordinate sequence, the control point enlarged coordinate sequence and the interpolation enlarged coordinate sequence are filtered for gross error points based on the RANSAC method with δ as the limit difference to form new control point enlarged coordinate sequences and interpolation enlarged coordinate sequences, the elevations corresponding to the control point enlarged coordinate sequence and the interpolation enlarged coordinate sequence are subtracted to obtain an elevation difference sequence, all values of the elevation difference sequence are sorted, and the median is taken as the elevation correction amount Hjzh_m1n1 of the m1th row and n1th column block.

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

[0043] S31, set the block with valid elevation correction amount obtained in step S2, that is, the correction block, the number of blocks is Njzh; record the elevation correction amount corresponding to each correction block and the corresponding row and column sequence number 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 a width of Ct*Rt+Njzh, the values in the connection array are initialized to 0, indicating no connection; wherein, m2=0,…,Njzh-1;

[0045] For the i-th row of the connection array, if i

[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] Wherein, i is counted from 0, j=0,1,2,3;

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

[0052] S33, set an array ifjzh with a length of Ct*Rt+Njzh, which stores whether the m3-th block DSM corresponding to the m3 / Ct-th row and the m3%Ct-th column needs to be corrected, 0 represents no need, and 1 represents need; wherein, m3=0,…,Ct*Rt+Njzh-1;

[0053] S34, set two height correction amount arrays lstHjzh with a length of Ct*Rt+Njzh, which store the height correction amount of the last iteration, and the n3-th array element represents the height correction amount of the n3 / Ct-th row and the n3%Ct-th column block DSM; all elements of the two arrays are initialized to 0; wherein, n3=0,…,Ct*Rt+Njzh-1.

[0054] Preferably, the update strategy of each iteration is,

[0055] A1, set two height correction amount arrays Hjzh with a length of Ct*Rt+Njzh, and the n3-th element represents the height correction amount of the n3 / Ct-th row and the n3%Ct-th column block DSM; the first Ct*Rt elements of the two arrays are initialized to Hjzh_m1n1 calculated in step S2, where m1=n3 / Ct and n1=n3%Ct, and all elements after the Ct*Rt-th element are initialized to 0;

[0056] A2, set two correction number arrays jlnum with a length of Ct*Rt+Njzh, and all elements of the two arrays are initialized to 0;

[0057] A3, if the m4-th row and the n4-th column of the connection array are 1; wherein, m4=0,…,Ct*Rt-1, n4=m4+1,…,Ct*Rt-1; calculate the height correction amount:

[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 the loop calculation of all correction amounts, save the average correction amount to the lstHjzh array:

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

[0066] A4, the iteration update step is executed 100 times to end the iteration, and the height correction amount of the m1th row and the n1th column block is lstjzh[m1*Ct+n1].

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

[0068] S41, update each column corresponding to each row block; wherein the jth column height correction amount of the m1th row and the n1th column block 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 i-th row height correction of the m1-th row and n1-th column block 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 = R1-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 < R1-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 the present application are: 1. The method of the present application corrects the overall terrain of DSM through all control points, and improves the plane and elevation accuracy. 2. The method of the present application divides the large area DSM into blocks, and each single block is corrected according to the distribution of control points to obtain correction information with higher reliability. 3. The method of the present application realizes the overall correction of all blocks through global propagation of correction information through adjacent blocks by iterative correction. 4. The method of the present application obtains the correction information of each grid point of each block through two times of distance reciprocal interpolation in the row block direction and the column block direction. BRIEF DESCRIPTION OF DRAWINGS

[0083] Figure 1 is a flowchart of the method in the embodiment of the present application. DETAILED DESCRIPTION

[0084] In order to make the purpose, technical scheme and advantages of the present application more clear, the present application is further described in detail below with reference to the drawings. It should be understood that the specific embodiments described herein are only used to explain the present application and do not limit the present application.

[0085] As shown in Figure 1 , in the embodiment, a method for improving the accuracy of large area DSM by using control points is provided. The method adopts the strategy of combining block processing and iterative correction on the basis of overall correction, which can improve the accuracy of correction and ensure that there is no edge trace caused by block in the large area DSM after correction. In addition to correcting the large area DSM after splicing, the same method can also be used to correct the large area multi-scene DSM before splicing, that is, the single-scene DSM is treated as the block DSM of the present application. The method specifically includes the following contents:

[0086] I. Overall terrain correction

[0087] The elevation corresponding to the row and column coordinates of each control point in the target area on the DSM is obtained by using bilinear interpolation to obtain a second ground coordinate sequence. The control point magnification coordinate sequence and the interpolation magnification coordinate sequence obtained based on the ground three-dimensional coordinates in the control point coordinate file and the second ground coordinate sequence are used to obtain new control point magnification coordinate sequence and interpolation magnification coordinate sequence by using the RANSAC method, and the correlation coefficient between the corresponding elevation values of the two new coordinate sequences is calculated. The left upper corner origin coordinates of the DSM are modified based on the elevation value with the maximum correlation coefficient, the affine transformation coefficients between the elevation values of the two coordinate sequences are calculated based on the least square method, and the grid point elevation calculated according to the affine transformation coefficients is used to replace the original grid point elevation to update the elevation of all grid points of the whole DSM.

[0088] The grid DSM file and the control point coordinate file to be corrected are input data. The grid DSM file contains row number R, column number C, row direction resolution ΔB, column direction resolution ΔL, top-left corner calculation coordinates (L0, B0), wherein L0 is the longitude of the top-left corner grid point, and B0 is the latitude of the top-left corner grid point. The control point coordinate file contains N ground three-dimensional coordinates (Li', Bi', Hi') (Li', Bi', Hi' are the longitude, latitude and height of the ith point respectively, i' = 0,..., N-1). The search window size is set to M*N, M and N are both odd numbers, and the search step is ΔS, the unit is degree.

[0089] 1.1, according to the longitude and latitude of each control point, the row and column coordinates (ri', ci') of the point on the DSM are calculated, and the height H1i' corresponding to the row and column coordinates is obtained by bilinear interpolation according to the row and column coordinates. Thus, a second ground coordinate sequence (Li', Bi', H1i') (i' = 0,..., N-1) is formed. 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 coordinates sequences (referred to as control point enlarged coordinate sequence and interpolation enlarged coordinate sequence respectively). The control point enlarged coordinate sequence and the interpolation enlarged coordinate sequence are filtered based on the RANSAC method with δ as the limit error to form new control point enlarged coordinate sequence and interpolation enlarged coordinate sequence. The correlation coefficient θ between the height values of the two coordinate sequences is calculated.

[0090] 1.2, the control point longitude and latitude coordinates (Li', Bi') are translated and transformed into (L1i', B1i'):

[0091] L1i' = Li' + (m-M / 2)*ΔL

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

[0093] For each control point coordinate sequence corresponding to the m, n combination, the above method is used to calculate a correlation coefficient θ m，n (m, n are the height values corresponding to the new control point enlarged coordinate sequence and the interpolation enlarged coordinate sequence respectively; m = 0, 1,..., M-1; n = 0, 1,..., N-1).

[0094] Find the m, n combination corresponding to the maximum correlation coefficient, assume it is (mmax, nmax), modify the top-left corner calculation coordinates of the DSM as (L10, B10):

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

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

[0097] With the interpolation longitude and latitude as variables, the affine transformation coefficients (a, b, c, d) between the two groups of high program elevation values corresponding to (mmax, nmax) are calculated based on the least square method, and the affine transformation coefficients conform to the following elevation transformation formula:

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

[0099] In addition to the invalid grid points identified by -999, the elevation of each grid point of the DSM is calculated according to the row number and column number and the longitude and latitude of the upper left corner after translation, the transformed grid point elevation is calculated according to the elevation transformation formula, and the original grid point elevation is replaced by the transformed grid point elevation, thereby updating all grid point elevations of the entire DSM.

[0100] II. Correction parameter calculation

[0101] The target area is divided into regular grids, and the control points inside each grid are counted. According to the relationship between the number of control points inside each grid and the latitude and longitude range of the control points and the corresponding threshold, the elevation correction amount of each grid is determined.

[0102] 2.1, Set the initial block interval of longitude and latitude direction as Lt0 and Bt0, and calculate the actual block number Ct and Rt of longitude and latitude direction:

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

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

[0105] According to Ct and Rt, the actual block interval Lt and Bt of longitude and latitude direction is calculated:

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

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

[0108] Wherein, the number of rows of each block Rs = Bt / Rt + 1, and the number of columns Cs = Lt / Ct + 1.

[0109] 2.2, Set the overlap of adjacent two blocks in longitude and latitude direction as Lc and Bc, then the longitude range of the m1th row and n1th column block (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 the m1th row and the n1th column is,

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

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

[0115] 2.3. According to the longitude and latitude ranges of the block in the m1th row and the block in the n1th column, the control points located in the range are counted, and the longitude range dL_m1n1 and the latitude range dB_m1n1 of the control points in the range are counted.

[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 after obtaining the interpolation height sequence and the height sequence of the control points in the range by the method of step S1, the elevation difference sequence is obtained by subtracting the elevations corresponding to the two height sequences, the values of the elevation difference sequence are sorted, and the median is taken as the height correction Hjzh_m1n1 of the block in the m1th row and the n1th column. If the number of control points is less than a threshold Nmin / 2, the height correction of the block in the m1th row and the n1th column 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 according to the longitude and latitude of each control point in the range, the row and column coordinates (ri', ci') of the point on the DSM are calculated, the elevation H1i' corresponding to the row and column coordinates is obtained by bilinear interpolation, a second ground coordinate sequence (Li', Bi', H1i') is formed, an error threshold δ is set, the longitude and latitude coordinates of the coordinate sequences (Li', Bi', Hi') and (Li', Bi', H1i') are multiplied by 100000 to form two new object coordinates sequences, which are the control point enlarged coordinate sequence and the interpolation enlarged coordinate sequence, the control point enlarged coordinate sequence and the interpolation enlarged coordinate sequence are filtered for gross error points based on the RANSAC method with δ as the limit difference to form new control point enlarged coordinate sequences and interpolation enlarged coordinate sequences, and the elevations corresponding to the control point enlarged coordinate sequences and the interpolation enlarged coordinate sequences are subtracted to obtain an elevation difference sequence, the values of the elevation difference sequence are sorted, and the median is taken as the height correction Hjzh_m1n1 of the block in the m1th row and the n1th column.

[0118] III. Correction parameter update

[0119] Iterative updating is performed on the obtained grid of effective height correction values to obtain the height correction values of the grid after iterative updating.

[0120] 3.1, assuming that the number of blocks of effective height correction values obtained in the second step is Njzh, hereinafter referred to as correction blocks, the height correction value corresponding to each correction block and the corresponding row and column number are recorded for each correction block, and the record information of the m2th (m2=0,…, Njzh-1) correction block is Inf_m2(Hjzh_m1n1, m1, n1). A two-dimensional connection array (hereinafter referred to as a connection array) with a length of Ct*Rt+Njzh and a width of Ct*Rt+Njzh is set, and the values in the array are initialized to 0, indicating no connection.

[0121] 3.2, for the i-th row (i is counted from 0) of the connection array, if i

[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, n1i_j is not less than 0 and not greater than Ct-1, the value of the i-th row and the m1i_j*Ct+n1i_j-th column of the connection array is set to 1, indicating that there is a connection. At the same time, all correction information Inf_m2 is searched, if the m1 of the information record is equal to m1i_j and n1 is equal to n1i_j, the value of the i-th row and the Ct*Rt+m2-th column of the connection array is set to 1.

[0127] 3.3, set an array ifjzh with a length of Ct*Rt+Njzh, store whether the DSM corresponding to the m3th array element (m3=0,…, Ct*Rt+Njzh-1) needs to be corrected, 0 represents no need, 1 represents need.

[0128] Set two height correction value arrays lstHjzh with a length of Ct*Rt+Njzh, store the height correction values of the last iteration, the n3th element (n3=0,…, Ct*Rt+Njzh-1) represents the last height correction value of the DSM of the n3 / Ct row block and the n3%Ct column block. All elements of the two arrays are initialized to 0.

[0129] In this embodiment, the update strategy of each iteration is:

[0130] A1, set two arrays Hjzh with length of Ct*Rt+Njzh, the nth3 element (n3=0,…,Ct*Rt+Njzh-1) represents the height correction of the n3 / Ct row block and n3%Ct column block DSM. The first Ct*Rt elements of the two arrays are initialized as Hjzh_m1n1 calculated in step S2, where m1=n3 / Ct and n1=n3%Ct, and all elements after the Ct*Rt elements are initialized as 0.

[0131] A2, set two arrays jlnum with length of Ct*Rt+Njzh, all elements of the two arrays are initialized as 0.

[0132] A3, if the m4th row (m4=0,…,Ct*Rt-1) and n4th column (n4=m4+1,…,Ct*Rt-1) of the connection array is 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 calculating all the corrections, save the average correction to the lstHjzh array:

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

[0141] A4, the above iteration update step is executed 100 times to end the iteration. The height correction of the m1th row and n1th column block (m1=0,…,Rt-1; n1=0,…,Ct-1) is lstjzh[m1*Ct+n1].

[0142] Four, update the entire DSM

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

[0144] 4.1, update each column corresponding to each row block; wherein the elevation correction amount of the jth column of the m1th row and the n1th column block is Hjzh_m1_n1_j;

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

[0146] If j is less than Cs / 2, Hjzh_m1_n1_j=0; if j is greater than or equal to Cs / 2, 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 is greater than or equal to Cs / 2, Hjzh_m1_n1_j=0; if j is less than Cs / 2, 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 is less than Ct-1:

[0150] If j is greater than or equal to Cs / 2, Hjzh_m1_n1_j=((1-(j-Cs / 2) / Cs)*lstjzh[m1*Ct+n1]+(j-Cs / 2) / Cs*lstjzh[m1*Ct+n1+1]; if j is less than Cs / 2, 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 ith row of the m1th row and the n1th column block is Hjzh_m1_n1_i;

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

[0153] If i is less than Rs / 2, Hjzh_m1_n1_i=0; if i is greater than or equal to Rs / 2, Hjzh_m1_n1_i=(1-(i-Rs / 2) / Rs)*lstjzh[m1*Ct+n1]+(i-Rs / 2) / Rs*lstjzh[(m1+1)*Ct+n1];

[0154] (2) when m1=Rt-1:

[0155] if i >= Rs / 2, then Hjzh_m1_n1_i = (1-i / Rs)*lstjzh[m1*Ct+n1]+i / 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];

[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)*lstjzh[m1*Ct+n1]+i / 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 technical solutions disclosed by the application, the following beneficial effects are obtained:

[0159] The application provides a method for improving the precision of a large-area DSM by using control points, which improves the plane and height precision by correcting the overall terrain of the DSM by all control points. The method divides the large-area DSM into blocks, and separately corrects each block according to the distribution of the control points to obtain correction information with higher reliability for some blocks. The method realizes the overall correction of all blocks by iteratively correcting and globally propagating the correction information through adjacent blocks, and obtains the correction information of each grid point of each block by twice distance reciprocal interpolation in the row block direction and the column block direction.

[0160] The above only describes the preferred embodiments of the application, and it should be noted that, for those skilled in the art, without departing from the principles of the application, some improvements and refinements can be made, and these improvements and refinements should also be considered as the protection scope of the application.

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 correction count arrays jlnum 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 at 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 at 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 at 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