Sliding Windows Method Based on Terrain Self-Similarity for Higher DEM Resolution in Flood Simulating Modeling

: A digital elevation model (DEM) is a quantitative representation of terrain and an important tool for Earth science and hydrological applications. A high-resolution DEM provides accurate basic Geodata and plays a crucial role in related scientiﬁc research and practical applications. However, in reality, high-resolution DEMs are often difﬁcult to obtain. Due to the self-similarity present within terrains, we proposed a method using the original DEM itself as a sample to expand the DEM using sliding windows method (SWM) and generate a higher resolution DEM. The main processes of SWM include downsampling the original DEM and constructing mapping sets, searching for the optimal matching, window replacement. Then, we repeat these processes with the small-scale expansion factor. In this paper, the grid resolution of the Taitou Basin was expanded from 30 to 10 m. Overall, the superresolution reconstruction results showed that the method could achieve better outcomes than other commonly used techniques and exhibited a slight deviation (root mean square error ( RMSE ) = 3.38) from the realistic DEM. The generated high-resolution DEM prove to be signiﬁcant in the application of ﬂood simulation modeling.


Introduction
DEM is a quantitative representation of the Earth's surface, which provides basic information about the terrain relief [1]. DEM and its derived attributes (slope, aspect, drainage area and network, curvature, topographic index, etc.) are important parameters for assessment of any process using terrain analysis and prerequisite in different applications such as flood simulation and management, landform analysis, terrain visualization and mapping [2][3][4][5][6]. A high-resolution DEM provides accurate basic Geodata and plays a crucial role in related scientific research and practical applications [7].
DEM is generated using different techniques such as LiDAR technology [8,9], the photogrammetric method using stereo data [10,11], aerial photographs [12,13] or interferometry [14]. The acquisition of quality DEM data over a large area is a challenging task because of its complicated generation process. In addition, there are many available opensource DEMs, including the Shuttle Radar Topography Mission (SRTM, 1 for the USA and 3 for other areas) [15], the Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM, 30 m) [16] and the Global 30 Arc-Second Elevation (GTOPO, 30 which is ∼1000 m) model [17], that still cannot meet the finer resolution requirements of some applications [18].
Currently, interpolation techniques, which include nearest-neighbor interpolation (NNI), bilinear interpolation (BI), inverse distance weighting (IDW), ordinary kriging (OK), and simulated, then rainfall thresholds can be estimated and the recurrence of flash floods can be better avoided.
We use the 30 m precision DEM data as the input source data. At the same time, 10 m precision DEM data are used to verify and support the effectiveness of the methods. The implemented 10 m precision DEM data are unpublished but important research data provided by the Hebei Meteorological Disaster Prevention Center and were collected after the severe mountain flood disaster occurred in the Taitou Basin.

Experiment Flowchart
As shown in Figure 2, the main processes of the experiment include the following steps. We use the 30 m precision DEM data as the input source data. At the same time, 10 m precision DEM data are used to verify and support the effectiveness of the methods. The implemented 10 m precision DEM data are unpublished but important research data provided by the Hebei Meteorological Disaster Prevention Center and were collected after the severe mountain flood disaster occurred in the Taitou Basin.

Experiment Flowchart
As shown in Figure 2, the main processes of the experiment include the following steps.

Mapping Set Construction
For an input sequence, several values are taken at intervals such that the new sequence is the downsampling of the original sequence. At this time, we compress DEM I into a lower resolution DEM I λ according to parameter λ using MATLAB software and DEM I has a size of M × N. After λ downsampling operations are performed, a low-resolution DEM I λ of size M/λ × N/λ will be obtained.
Then, we calculate the slopes of DEMs I and I λ to obtain the corresponding digital slope models S and S λ , respectively.  mapping sets, searching for the optimized matching and expansion. During the expansion process, each small-sized window is replaced by its best match window's corresponding larger window according to its coordinates, thereby constructing a higher-resolution DEM. Then, we repeat these processes with the small-scale expansion factor until the dem resolution reach the requirement.

Mapping Set Construction
For an input sequence, several values are taken at intervals such that the new sequence is the downsampling of the original sequence. At this time, we compress DEM I into a lower resolution DEM according to parameter λ using MATLAB software and DEM I has a size of × . After λ downsampling operations are performed, a low-resolution DEM of size / × / will be obtained. Then, we calculate the slopes of DEMs I and to obtain the corresponding digital slope models S and , respectively.
The acquisition of the digital slope model requires calculating the overall terrain slope based on the DEM elevation data. The slope of the surface is the angle between the normal direction of the surface tangent and the Z axis; therefore, the slope ( , ) can be calculated as where Slope ( , ) represents the rate of elevation change of the terrain surface. ( , ) represents the rate of elevation change in the north-south direction and ( , ) represents the rate of elevation change in the east-west direction. The common slope extraction algorithms used to calculate ( , ) and ( , ) include the simple difference, three-order unweighted difference and third-order inverse Figure 2. Flowchart of the sliding windows method (SWM). The main processes of the experiment include constructing mapping sets, searching for the optimized matching and expansion. During the expansion process, each small-sized window is replaced by its best match window's corresponding larger window according to its coordinates, thereby constructing a higher-resolution DEM. Then, we repeat these processes with the small-scale expansion factor until the dem resolution reach the requirement.
The acquisition of the digital slope model requires calculating the overall terrain slope based on the DEM elevation data. The slope of the surface is the angle between the normal direction of the surface tangent and the Z axis; therefore, the slope S(i, j) can be calculated as where Slope S(i, j) represents the rate of elevation change of the terrain surface. f x (i, j) represents the rate of elevation change in the north-south direction and f y (i, j) represents the rate of elevation change in the east-west direction. The common slope extraction algorithms used to calculate f x (i, j) and f y (i, j) include the simple difference, three-order unweighted difference and third-order inverse distance squared weighted difference. The third-order inverse distance weight difference is used in our algorithm because of its higher calculation accuracies [44,45]. f x (i, j) and f y (i, j) can be, respectively, expressed as Remote Sens. 2021, 13, 3604 where g is the resolution of the grid and z(i, j) is the center value of the grid in a 3 × 3 window. Digital slope model calculation instructions are shown in Figure 3. f x (i, j) and f y (i, j) are calculated by collecting the elevations of the eight points around z(i, j) in a 3 × 3 gride (Equations (2) and (3)) and then we can obtain slope S(i, j) (Equation (1)). By traversing each grid as a center point, we complete the slope calculation of the entire DEM based on the MATLAB platform and, finally, obtain the digital slope model corresponding to the DEM.
Remote Sens. 2021, 13, x FOR PEER REVIEW 5 of 15 distance squared weighted difference. The third-order inverse distance weight difference is used in our algorithm because of its higher calculation accuracies [44,45].
where g is the resolution of the grid and ( , ) is the center value of the grid in a 3 × 3 window. Digital slope model calculation instructions are shown in Figure 3. ( , ) and ( , ) are calculated by collecting the elevations of the eight points around ( , ) in a 3 × 3 gride (Equations (2) and (3)) and then we can obtain slope ( , ) (Equation (1)). By traversing each grid as a center point, we complete the slope calculation of the entire DEM based on the MATLAB platform and, finally, obtain the digital slope model corresponding to the DEM. Processes and descriptions of sliding-window expansion. Take λ = 2 as an example. In the first step, we compress DEM I into a lower resolution DEM according to parameter 2. Due to downsampling, for each 2 × 2 sliding window in DEM and digital slope model , there is a corresponding 3 × 3 window in original DEM I and digital slope model S according to its coordinates. In the second step, for any window of the normalized set , we can search each window of set to obtain the optimal matching results. The matching search step includes height data matching and slope . Processes and descriptions of sliding-window expansion. Take λ = 2 as an example. In the first step, we compress DEM I into a lower resolution DEM I λ according to parameter 2. Due to downsampling, for each 2×2 sliding window in DEM I λ and digital slope model S λ , there is a corresponding 3 × 3 window in original DEM I and digital slope model S according to its coordinates. In the second step, for any window of the normalized set IW 1 , we can search each window of set IW λ to obtain the optimal matching results. The matching search step includes height data matching and slope data matching. In the third step, the 2 × 2 windows in set IW 1 can be replaced by the 3 × 3 windows in set IW 0 . The red and green frames represent different windows, represents known elevation data and is the extended data obtained from window replacement. The black solid symbols and G represent the overlapping data present in two windows and the values in the overlapped part are averaged. The windows should recover the normalization after replacement. After the expansion of every window, a high-resolution DEM T can be reconstructed.
After downsampling the original DEM I to I λ and calculating the slope of them, we slide the windows on the low-resolution DEM I λ and the digital slope model S λ , divide them into fixed-size windows and then put the height matrix data and slope matrix data contained in the windows into the to-be-matched elevation set IW λ and to-be-matched slope set SW λ , respectively.
Since our purpose is to use the original DEM itself as a sample to expand the DEM, we slide the windows on the original DEM I and the digital slope model S, divide them into fixed-size windows and then put the height matrix data and slope matrix data contained in the windows into the sample elevation set IW 0 and sample slope set SW 0 , respectively. In addition, the original DEM I needs to be expanded, so we put the height matrix data and slope matrix data contained in the windows from DEM I into the to-be-extended elevation set IW 1 and to-be-extended slope set SW 1 , respectively. It can be seen that the to-beextended elevation set IW 1 is equal to the sample elevation set IW 0 . Due to downsampling, for each n × n sliding window in DEM I λ , there is a corresponding (λ(n − 1) + 1) × (λ(n − 1) + 1) window in original DEM I according to its coordinates. Therefore, there are also mapping relationships f between the data contained in the windows from IW λ and the data contained in the windows from original IW 0 , which can be expressed as

Search for the Optimal Matching
Each matrix in sets should be normalized to extract the rules of elevation or slope changes. Only after this step, can the similarities on elevation and slope between the matrix from set IW 1 and the matrix from set IW λ to be calculated. We consider a n × n (n < N λ and n < M λ ) sliding window that contains an n-dimensional matrix A. The normalization process needs to search the maximum and minimum values in matrix A, where the maximum value is denoted as Max and the minimum value is denoted as Min. The normalization process is expressed as where E is n × n identity matrix. Through the normalization operation, the original data in the matrix can be transformed into the data better reflecting the change rules of the whole matrix. Then, we can obtain normalized sets IW 0 , SW 0 , IW 1 , SW 1 , IW λ and SW λ by normalizing the matrices in sets IW 0 , SW 0 , IW 1 , SW 1 , IW λ and SW λ . In addition, the set IW λ and the set IW 0 also have the mapping relationships f .
Let two normalized elevation matrices are A and B and the normalized slope matrices that correspond to the matrices A and B are S A and S B . We can calculate the sum of the squared elevation difference and the sum of the squared slope difference to obtain the overall sum of the squared difference. The smaller the overall sum of the squared difference is, the higher the similarities and correlation are and the higher the degree of matching is. The correlation Correlation A B between matrices A and B can be, respectively expressed as where D is the difference matrix of two matrices A and B . In addition, S D is the difference matrix of two matrices S A and S B .The elements of matrices D and S D are d (i, j) and s d (i, j), respectively. ω is the slope weight, which is adjusted empirically. The matching process is shown in Figure 3. During this process, the elevation and slope values are all needed, which ensures that both the elevation rules and slope rules in the two windows are as consistent as possible. For any matrix of the normalized to-beexpanded set IW 1 , we can search each matrix of set IW λ to obtain the optimal matching results with the highest correlation. Therefore, we can obtain a mapping relationship g between matrices from set IW 1 and set IW λ .

Window Replacement
Due to the optimal matching relationship g between matrices from set IW 1 and set IW λ and the mapping relationships f between matrices from set IW λ and set IW 0 , there are also mapping relationships between matrices from set IW 1 and set IW 0 . The n × n windows in set IW 1 can be replaced by the (λ(n − 1) + 1) × (λ(n − 1) + 1) windows in set IW 0 since they are larger and contain more elevation information. The windows in set IW 0 should recover the normalization after replacement.
Meanwhile, there are duplicate data between sliding windows created during replacement. In addition, if there are several columns (rows) of overlap between two initial to-be-expanded sliding windows, there will also be the same number of columns (rows) that overlap between two expanded windows and the values in the overlapped part are averaged. Therefore, the obtained image will have better continuity and smoothness.
As shown in Figure 3, due to the overlap of sliding windows, the black solid symbols and G represent the overlapping data present in two windows. Since G is already-known elevation data, in the third replacement operation, the to-be-extended data at need to be averaged on both sides.
After the expansion of every window, a high-resolution DEM T can be reconstructed from a low-resolution DEM I.
The overall process description of the above three steps is shown in Figure 3.

Small-Scale Expansion Factor
Since the expansion ratio may be large or the amount of sample data may be too small, multilevel segmented upsampling may be adopted to reduce the error in the expansion [46]. By using a small-scale expansion factor to repeat the expansion process step by step, a DEM with a specified resolution can finally be obtained.

Quantitative Evaluation
The mean error (ME) and root mean square error (RMSE) are often used as indicators to determine the reconstruction accuracy. The smaller the absolute values of ME and RMSE are, the better the reconstruction quality.
In this study, the estimated height (ẑ i ) derived from the SWM and the selected interpolation technique were compared at each point to the observed height (z i ) using ME and RMSE.
In addition, for the different interpolation methods, we calculate the PEP1.5 (percentage of error points with an elevation error greater than 1.5%) for the whole image and we mark the error points with red dots on the generated image to ensure that we can intuitively see the number and distribution of the error points generated by each method.

Visual Evaluation
In this study, the interpolated result's degree of smoothing to the realistic DEM is the main indicator for judging visual quality [29]. The smoothness is measured considering the topographic roughness, which is a secondary terrain parameter derived from the DEM that is used in the geosciences and environmental studies (Equation (13), [47,48]). A large roughness value often indicates a minor smoothness. To conform to the overall visual perception ability of humans, we used the 3 × 3 windows' average roughness instead of the roughness obtained at the pixel level for visual evaluation. The roughness difference between the realistic DEM and the generated results is calculated by Equation (12).
whereŝ l i is the estimated slope, n represents the window's pixels, "realistic" represents the realistic DEM and "estimated" represents the estimated DEM.

Simulated Flooding Event Evaluation
The goodness of fit index (FIT A ) is used to evaluate the influence of DEM resolution on the calculated flood extent [49]: where FIT A is a goodness-of-fit index and Faobs and Famod are the observed and modeled flood extents, respectively. A large FIT A value often indicates a better simulation effect.

Parameters of Sliding Windows
In this paper, we took the Taitou Basin as the study area and gradually expanded the DEM resolution from 30 to 10 m in MATLAB. The parameter λ was 2 and the size of the sliding windows was 2 × 2. Each sliding window corresponded to a 3 × 3 sample window of the original DEM. Meanwhile, we tried different values and found that the best result was obtained when the slope weight w was 0.1. This mapping method expands the original DEM by two times. After repeating this operation twice, we obtain a 7.5 m precision DEM with four times better accuracy compared to the original DEM, to which we can then apply downsampling to obtain a 10 m precision DEM.

Image Generation
In Figure 4, compared with the 30 m precision original contour image, the 10 m precision contour image after expansion has more textural information. Meanwhile, the 10 m precision image generated by the SWM and the realistic contour image largely have consistent contour curves, indicating that the SWM performs well and produces a terrain that largely fits the realistic situation.

Accuracy for Altitude Estimation
We selected interpolation techniques with a local neighborhood or geostatistical approach since they are commonly used in geomorphological research. The techniques include NNI, BI, IDW and OK. In addition, we evaluated the performances of the different techniques by comparing the different results generated by these methods.
Remote Sens. 2021, 13, x FOR PEER REVIEW 9 of 15 precision image generated by the SWM and the realistic contour image largely have consistent contour curves, indicating that the SWM performs well and produces a terrain that largely fits the realistic situation.

Accuracy for Altitude Estimation
We selected interpolation techniques with a local neighborhood or geostatistical approach since they are commonly used in geomorphological research. The techniques include NNI, BI, IDW and OK. In addition, we evaluated the performances of the different techniques by comparing the different results generated by these methods.

Error Indicators
The mean error (ME) and the root mean square error (RMSE) between estimations and observations of altitude at the study sites are presented in Table 1. Compared with other methods, the SWM has the smallest RMSE and ME, which indicates that the SWM better represents the realistic DEM than the other tested interpolation methods. Thus, the accuracy and reconstruction quality of the SWM are better than those of other commonly used algorithms when reconstructing a high-resolution DEM.
The PEP1.5 s in the whole images generated by different methods are shown in Table 2. Figure 5 and Table 2 illustrate that IDW and NNI have significantly larger PEP1.5 values than the other interpolation techniques. NNI uses simple assignment operations, which is prone to large errors. IDW is easily affected by extreme values and performs well in uniform samples. However, it is not suitable for mountainous areas. Therefore, the performance of IDW is not so good in this experiment. In addition, OK has achieved better results based on the function model. However, it still has more error points with elevation errors greater than 1.5% of the boundary of the study area than BI and the SWM. BI and the SWM seem to have more error points greater than 1.5% on the edge and steep areas than flat areas, but we can see that the error points with elevation errors greater than 1.5% of the SWM account for the lowest proportion because SWM consider both the elevation

Error Indicators
The mean error (ME) and the root mean square error (RMSE) between estimations and observations of altitude at the study sites are presented in Table 1. Compared with other methods, the SWM has the smallest RMSE and ME, which indicates that the SWM better represents the realistic DEM than the other tested interpolation methods. Thus, the accuracy and reconstruction quality of the SWM are better than those of other commonly used algorithms when reconstructing a high-resolution DEM.
The PEP1.5 s in the whole images generated by different methods are shown in Table 2.  Figure 5 and Table 2 illustrate that IDW and NNI have significantly larger PEP1.5 values than the other interpolation techniques. NNI uses simple assignment operations, which is prone to large errors. IDW is easily affected by extreme values and performs well in uniform samples. However, it is not suitable for mountainous areas. Therefore, the performance of IDW is not so good in this experiment. In addition, OK has achieved better results based on the function model. However, it still has more error points with elevation errors greater than 1.5% of the boundary of the study area than BI and the SWM. BI and the SWM seem to have more error points greater than 1.5% on the edge and steep areas than flat areas, but we can see that the error points with elevation errors greater than 1.5% of the SWM account for the lowest proportion because SWM consider both the elevation rules and slope rules. Therefore, the error of SWM is smaller than other interpolation techniques. In general, the SWM is more suitable for complex terrain and performs better than other interpolation techniques. other interpolation techniques.

Visual Comparison
We extract the regions from the reconstructed DEM grayscale image and contour map and then compare them with the corresponding part of the realistic 10 m DEM image.
The DRs covering the entire images generated by different methods are shown in Table 3.

Visual Comparison
We extract the regions from the reconstructed DEM grayscale image and contour map and then compare them with the corresponding part of the realistic 10 m DEM image.
The DRs covering the entire images generated by different methods are shown in Table 3. As shown in Figure 6, the high-resolution DEM data reconstructed by the traditional interpolation algorithm easily lose a considerable amount of terrain detail information, while the image generated by the SWM provides more realistic and detailed information. From a visual point of view, the expanded DEM image generated by the SWM is closer to the true 10 m DEM elevation data. In addition, as we can see, generally, the generated results with less smoothness appear to be more similar to the realistic DEM and DR values that close to zeros reveal this well, which we consider obtain better visual quality. Table 3 shows that the SWM result has a smaller absolute DR value, while the proposed methods have larger absolute DR values, which indicates that the SWM generates more realistic textures and produces better visual quality than the proposed methods.
As shown in Figure 6, the high-resolution DEM data reconstructed by the traditional interpolation algorithm easily lose a considerable amount of terrain detail information, while the image generated by the SWM provides more realistic and detailed information. From a visual point of view, the expanded DEM image generated by the SWM is closer to the true 10 m DEM elevation data. In addition, as we can see, generally, the generated results with less smoothness appear to be more similar to the realistic DEM and DR values that close to zeros reveal this well, which we consider obtain better visual quality. Table 3 shows that the SWM result has a smaller absolute DR value, while the proposed methods have larger absolute DR values, which indicates that the SWM generates more realistic textures and produces better visual quality than the proposed methods.

Results of the Madian Basin
As we can see that SWM works well in the Taitou Basin, we apply SWM on another terrain case to test whether this model is still useful. The original input data are based on the 30 m resolution DEM of the Madian Basin. We gradually expanded the DEM resolution from 30 to 10 m using NNI, BI, IDW and SWM techniques.
As shown in Table 4, compared with other methods, the SWM has the smallest RMSE, ME, PEP1.5 and DR, which indicates that the SWM represents the realistic DEM better than the other tested interpolation methods. Thus, it is concluded that SWM has an advantage in quantitative accuracy and generalization capability when reconstructing a high-resolution DEM.

Application of High-Resolution DEMs in Flood Modeling
We simulated the flash flood disaster that occurred in the Taitou Basin on 19 July 2016 by using a DEM with different resolutions. Since there were no flow data directly observed in the basin, the HEC-HMS model was first used to simulate the process of rainfall runoff and then the hydrograph of the July 2016 flash flood produced by the HEC-HMS model was used as the input data of the FLO-2D model to simulate the inundation of this event.
In Figure 7, it can be seen intuitively that the flood extent and the difference compared to the actual flood extent decrease with increasing DEM resolutions. As shown in Table 5, the flood extent decreases gradually with increasing water depth under the same DEM resolution. With the increase in DEM resolution, the flood extent changes slightly in areas with water depths between 1 m and 3 m, decreases considerably in areas with water depths < 1 m and increases considerably in areas with water depths ≥ 3 m. Therefore, the flood extent of areas with extreme water depths (<1 m and ≥3 M) exhibits a notable response to changes in the DEM resolution. In addition, when the resolution of the DEM changes from 30 m to 10 m, the total flood range of the basin is reduced from 0.35 km 2 to 0.25 km 2 and FIT A also increases from 0.56 to 0.74. The experimental results indicate that a higher resolution DEM can improve the accuracy of the flood extent simulation.

Conclusions and Recommendations
This paper proposed a method based on sliding windows to improve the resolution of DEMs. Traditional interpolation methods tend to lose large amounts of detailed terrain information when they expand DEMs. Based on the geographical self-similarity rule, we used the DEM itself as a sample to expand the original DEM through sliding windows and generate a higher resolution DEM. When searching for the best matching window, both elevation and slope information from two windows were considered and the matching rule was adjusted for different terrains to better explore the data trends and obtain the best results.
After that step, we used the Taitou Basin as the study area to evaluate the accuracy of the SWM and chose several common interpolation techniques used to generate DEMs (NNI, BI, IDW and OK) to evaluate the performance by comparing their different results. The final experimental results showed that the high-resolution DEM data reconstructed by the traditional interpolation algorithm easily lose a considerable amount of terrain detail information. OK and BI yielded better estimates flat areas but did not perform well on edges or in steep areas. NNI exhibited a large deviation and was not suitable for practical

Conclusions and Recommendations
This paper proposed a method based on sliding windows to improve the resolution of DEMs. Traditional interpolation methods tend to lose large amounts of detailed terrain information when they expand DEMs. Based on the geographical self-similarity rule, we used the DEM itself as a sample to expand the original DEM through sliding windows and generate a higher resolution DEM. When searching for the best matching window, both elevation and slope information from two windows were considered and the matching rule was adjusted for different terrains to better explore the data trends and obtain the best results.
After that step, we used the Taitou Basin as the study area to evaluate the accuracy of the SWM and chose several common interpolation techniques used to generate DEMs (NNI, BI, IDW and OK) to evaluate the performance by comparing their different results. The final experimental results showed that the high-resolution DEM data reconstructed by the traditional interpolation algorithm easily lose a considerable amount of terrain detail information. OK and BI yielded better estimates flat areas but did not perform well on edges or in steep areas. NNI exhibited a large deviation and was not suitable for practical applications. The SWM had similar defects, but it was still superior to other common algorithms in both subjective visual effects and reconstruction indicators and the image generated by the SWM has more detailed information. The sliding window is a good choice for geographic information systems (GIS) specialists to generate higher resolution DEMs.
Based on a generated high-resolution DEM, a digital basin can be generated and an effective hydrological model can be applied to the DEM grid to calculate the process of rainfall runoff and inundation. According to the simulation of the flood extent and water depth distribution, once the warning limit is reached, flood warnings can be sent.
However, the accuracy of the DEMs generated by the SWM can be improved furthermore. In practical applications, high-resolution DEMs of local areas can be obtained in many ways. Therefore, in our future work, we will consider adding partial high-resolution DEMs as samples to extract the terrain rule and use adaptive methods to update better window matching rules so that the original low-resolution DEMs can be more accurate.