Modeling Change of Topographic Spatial Structures with DEM Resolution Using Semi-Variogram Analysis and Filter Bank

Abstract: In this paper, the way topographic spatial information changes with resolution was investigated using semi-variograms and an Independent Structures Model (ISM) to identify the mechanisms involved in changes of topographic parameters as resolution becomes coarser or finer. A typical Loess Hilly area in the Loess Plateau of China was taken as the study area. DEMs with resolutions of 2.5 m and 25 m were derived from topographic maps with map scales of 1:10,000 using ANUDEM software. The ISM, in which the semi-variogram was modeled as the sum of component semi-variograms, was used to model the measured semi-variogram of the elevation surface. Components were modeled using an analytic ISM model and corresponding landscape components identified using Kriging and filter bank analyses. The change in the spatial components as resolution became coarser was investigated by modeling upscaling as a low pass linear filter and applying a general result to obtain an analytic model for the scaling process in terms of semi-variance. This investigation demonstrated how topographic structures could be effectively characterised over varying scales using the ISM model for the semi-variogram. The loss of information in the short range components with resolution is a major driver for the observed change in derived topographic parameters such as slope. This paper has helped to quantify how information is distributed among scale components and how it is lost in natural terrain surfaces as resolution becomes coarser. It is a basis for further applications in the field of geomorphometry.


Introduction
Terrain variables such as slopes, curvatures and catchment structures are primary parameters for applications in hydrology, soil erosion, soil water and carbon budgets [1][2][3][4][5].In recent years, Digital Elevation Models (DEMs) have become the most common source of terrain data for such applications.However, the scale of available data and the scale at which process understanding and modeling is to be applied may not always match.For example, it has been widely reported that terrain variables are sensitive to the resolution of the DEM used [6,7] and the effect of resolution on terrain variables can be significant for these applications.Relationships between changes in terrain variables with DEM resolution have consequently been extensively (but mostly empirically) researched in recent years [8][9][10][11].When available DEM data are coarser than required, it is useful to study scaling effects and rescaling methods to know what resolution should be chosen for certain applications and what to do if no data at suitable resolutions are available.Some efforts have been put into how to rescale coarser resolutions to estimate terrain parameters in finer ones [8,12].As the resolution of available DEMs increases, the issue of matching DEM and processing scale has not become less important as DEM scale is not only a matter of grid cell resolution but also the information content of the DEM down to the finest scales is provided by the resolution.
To advance this area of research, it is useful to find out how topographic variables are affected by scale and what kinds of information are lost in the process of generating or generalizing a DEM.If this can be done, provided the primary terrain spatial structures making up a certain area can also be established, it may be possible to characterize loss of information in these terrain structures in (for example) upscaling.Conversely, it may be possible to characterize terrain type by the way information changes with scale [13].In this paper, the authors have taken a step towards this by modeling topographic structures in DEMs and the change in information as resolution becomes coarser or finer.Geostatistical analysis [14], multi-component semi-variogram methods [15] and analytical methods used in the theory of autocorrelation and regularization in digital images [16,17] were the primary tools in this study.The tools are applied using methods of digital image processing [3].
The semi-variogram is a basic tool in geostatistical analysis.Its primary development was in geological resource estimation and was later provided with a strong mathematical and statistical basis [18].It has developed to become a valuable tool in many geographical spatial studies such as remote sensing modeling [16,17], Land cover [19][20][21] and soil moisture [22].Semi-variogram analysis has also been applied successfully to topographic data in a number of ways.For example, to select an appropriate pixel size for DEMs interpolated from contour data [23], to compute the fractal dimension of a DEM [24], for interpolation using Kriging methods [25] and for investigating micro-scale surface structures [26].In this paper, semi-variograms are used to explore the topographic structures represented by DEMs at different scales and to help characterize landscapes in terms of the topographic structures found in the data.
It is clear that landscape processes have a range of characteristic scales and that it is usually not easy (or even possible) to model them by a single commonly defined (authorized) semi-variogram model.This issue has been discussed by Oliver and others in the field of image exploration [21,27].In this paper, an "Independent Structures Model" (ISM), commonly called the "Nested Model" e.g., [27], has been used to quantitatively analyze these scales.Using ISM, the topographic structures may be identified and the way they change with DEM resolution can be explored through the changes that occur in the components in the ISM using methods presented in this paper.Using terminology consistent with Serra suggests that the name "Nested" models is a more appropriate description for a hierarchical model rather than for a sum of independent component semi-variograms [28].Hence, the term "ISM" (based on Serra [28]) has been used in the following pages.
The paper outlines the form of geostatistical model most useful for this task.It then develops new analytic methods and a general result that allows the effects of scaling to be quantified.It also introduces tools to help identify components and account for trends then demonstrates the tools with some practical examples.Specifically, in Section 2, the materials used, including background information and DEMs, are described and the mathematical methods applied to the DEMs are also outlined as needed.In Section 3, the methods and new analytic models developed here are used to derive an ISM for typical data sets and demonstrate how it can model the semi-variograms of various DEMs of the same area at different scales using the same underlying ISM model.A tool called the "Filter Bank" is introduced to help identify the finest landscape components of the ISM model and Kriging analysis is used to create images of the scale components.In Section 4, some issues that arose are discussed.The conclusions are drawn in Section 5.This paper shows how when DEM resolution becomes coarser (in this case from 2.5 m and 25 m), the detailed components are more influenced by the change in scale than the broader scale features.In the example selected, for topographic structures within a range of 70 m the variance associated with the component was reduced by 31.1%;variance of topographic structures larger than 70 m were reduced much less and those above 1500 m not at all.The results of this paper help in the understanding of the upscaling process in topographic analysis, offer tools with which to model scale components of a DEM and help in defining an appropriate DEM resolution to match applications.

The Study Area
The study area is located in the Loess Hilly-Gully Region in the Loess Plateau (LP).It is a square of side 8.1 km.There is a significant watershed named the Jiu Yuan watershed crossing this study area in a northeast-southwest direction.It is a typical Loess-hilly area with complex terrain surface and intensive soil erosion.This type of region has complex terrain with many gullies, steep slopes and short slope lengths.The location and complexity of the terrain surface is as indicated in Figure 1.

The Study Area
The study area is located in the Loess Hilly-Gully Region in the Loess Plateau (LP).It is a square of side 8.1 km.There is a significant watershed named the Jiu Yuan watershed crossing this study area in a northeast-southwest direction.It is a typical Loess-hilly area with complex terrain surface and intensive soil erosion.This type of region has complex terrain with many gullies, steep slopes and short slope lengths.The location and complexity of the terrain surface is as indicated in Figure 1.

Source Datasets
A topographic map with map scale of 1:10,000 was used in this paper as base data.It was derived from surveys based on aerial photography and used updated modern map making facilities during the 1990s by the Chinese government.Based on the 1:10,000 topographic map, DEMs with resolution of 2.5 m and 25 m were built using the software ANUDEM5.3 developed by Hutchinson [29,30].They have been referred to as DEM2.5 and DEM25 for short in this paper.These two examples are used here to demonstrate that the methods work as the theory predicts and some current and future applications are provided in the Discussion section.The input datasets to the ANUDEM software include contours, spot heights, streams and lakes.The contours were the primary data and the excessive profile curvature penalty was 0.7 (high) in the settings of the ANUDEM software.The common projection for all data is the Gauss-Kruger based on the Krasovsky 1940 Datum as commonly used in China.The grid cell shapes and sizes are all square and in meters.The primary study area is contained in a rectangle with size of 8.1 km.A buffer distance of 500 m was applied to all data in order to reduce edge effects.That is, the extent for all data used in data processing in this paper is actually a square with size 9.1 km but results only generated for the smaller study area.
It is important to note that the DEM data developed and used here are different from data traditionally analyzed by geostatistics in geology and ecology.They are not point samples but are the result of processing various forms of data into dense grids.The pre-processing affects the nature of

Source Datasets
A topographic map with map scale of 1:10,000 was used in this paper as base data.It was derived from surveys based on aerial photography and used updated modern map making facilities during the 1990s by the Chinese government.Based on the 1:10,000 topographic map, DEMs with resolution of 2.5 m and 25 m were built using the software ANUDEM5.3 developed by Hutchinson [29,30].They have been referred to as DEM2.5 and DEM25 for short in this paper.These two examples are used here to demonstrate that the methods work as the theory predicts and some current and future applications are provided in the Discussion section.The input datasets to the ANUDEM software include contours, spot heights, streams and lakes.The contours were the primary data and the excessive profile curvature penalty was 0.7 (high) in the settings of the ANUDEM software.The common projection for all data is the Gauss-Kruger based on the Krasovsky 1940 Datum as commonly used in China.The grid cell shapes and sizes are all square and in meters.The primary study area is contained in a rectangle with size of 8.1 km.A buffer distance of 500 m was applied to all data in order to reduce edge effects.That is, the extent for all data used in data processing in this paper is actually a square with size 9.1 km but results only generated for the smaller study area.
It is important to note that the DEM data developed and used here are different from data traditionally analyzed by geostatistics in geology and ecology.They are not point samples but are the result of processing various forms of data into dense grids.The pre-processing affects the nature of the post-processing approach as it does remotely sensed data, and DEMs produced by digital photogrammetry, Lidar data and Radar data.The data typically also form huge data files where image processing is an important way to analyze them.Tools from geostatistics can, however, be successfully generalized to handle these types of data [21].

Semi-Variogram Statistics
Variation in spatial properties, including that in DEMs, can be described using the theory and methods of geostatistics [31].The semi-variogram is a central tool of geostatistical analysis.It measures the spatial correlation between locations based on the principle that a spatial property is more likely to be similar at nearby locations than at distant ones [31] and that the similarity can be modeled using properties of spatially random functions.
The general formula used here to calculate a semi-variogram γ, for a (1D or 2D) spatial function Z at a given distance h [31] is: where γ phq is the semi-variogram; x i is the location of a certain point; Z px i q and Z px i `hq are the observed values of Z at locations x i and x i `h; for 2D data, h is a vector with both direction and distance, known generally as the lag; N phq is the number of pairs of points with a separation distance (or lag) of h.Values of the semi-variogram obtained by changing h make up the experimental semi-variogram which is sometimes called the "semi-variance".
In this formulation, a DEM surface is modeled as the realization of a random 2D function Z px, yq which is at least assumed to be stationary in increments [14] and in most cases assumed to be differentiable as well.The 2D semi-variogram for Z px, yq is defined in terms of increments in the x and y directions.That is: In a DEM image, the spatial grid coordinates are defined by the row and column values.So, the increments (i.e., lags) in the x and y directions can also be defined by increments in pixel (or grid cell) units in the row and column directions.The result in this case is a 2D semi-variogram image which is the primary information set for the 2D image data.
For practical reasons, using the fact that the semi-variogram is an even function (i.e., diagonally-symmetric) (γ p´hq " γ phq) makes it convenient to only compute the semi-variograms in four primary shift directions.The four directions may be denoted East-West (E-W), North-South (N-S), Northeast-Southwest (NE-SW) and Northwest-Southeast (NW-SE).The diagonal directions are interpolated to equal radial steps with the length of steps in the square directions so that lags are in radial steps.The average of these semi-variograms in the four directions constitutes a good estimate for the overall average radial semi-variogram since, because of the symmetries; it is the same as averaging the eight principal radial directions.
The overall average radial semi-variogram is particularly useful for this work as the experimental semi-variograms from small areas often contain local directional anisotropy and bias due to specific local orientation of terrain features or regional trends.Characterizing terrain for scale analysis seeks to be independent of such local directional factors and the radial average semi-variogram provides a useful "rotation invariant" information measure.Using this radial average information, the semi-variogram statistics have been used to explore the rotation invariant spatial information in the Hc-DEMs with varying resolutions.
An experimental semi-variogram is often described in terms of two basic derived parameters.These are Range (the distance beyond which the values are no longer correlated) and Sill (the semi-variogram value at distances greater than the range).For experimental semi-variograms, direct methods such as those described in [13] can be used to estimate the overall range and sill parameters as part of general data exploration.When the semi-variogram is "stationary in increments", the increments (Z `x `hx, y `hy ˘´Z px, yq) in Equation ( 2) are stationary in px, yq for given `hx, h y ˘.In general, this is associated with the semi-variogram reaching a sill but the sill may not occur within the extent of the data.
For the kinds of surface used in this paper, the variance is finite and the spatial structure can be equivalently defined by a covariance function C phq where: Taking issues such as trends (see later) into account, the relationship between the semi-variogram and the covariance function appropriate for DEMs can be written as: C p0q is the data variance or sill.Further use will be made of this relationship in the following sections where in order to model the radial average semi-variogram we will use radially symmetric variogram and covariance functions.The application of this model tends to imply a stronger form of stationarity than stationary in increments as often used for semi-variograms of geological data.

The Independent Structures Model
Spatial structures can be organized hierarchically at scales differing by several orders of magnitude from a few meters to hundreds of kilometers.The Independent Structures Model (ISM) was developed and used here to model these spatial structures.In the ISM model, the DEM Z pxq is considered to be the sum of N "independent" fields; the covariance between the N fields is assumed to be zero so that the total variance of Z pxq can be written as the sum of variances of the components.This may be expressed as: In this form of the model, the individual covariance terms (γ j pxq) are assumed to have zero mean and be normalized to unit variance so that the overall spatial covariance functions of Z (C Z ph)) and the semi-variogram function for Z pxq can be conveniently expressed as: In this expression, the component covariance functions C j `Rj ; h ˘correspond with each of the component fields Y j pxq and C j `Rj ; 0 ˘" 1.The range values "R" of the components are assumed to be ordered as: Written in this way, the first component is the one with most "detail" or high spatial frequency information and the last component is the one with greatest low frequency or regional information.
Because we are using radial average semi-variograms, the ISM models are based on radially symmetric covariance functions, in which case ranges are radial 2D distances and "h" is radial distance.The Gaussian model has been found to be suitable for regular surfaces, such as are often found among DEMs [32], and it is well known that sums of Gaussians provide very flexible models for many functions.For these reasons, a radially symmetric 2D Gaussian model was selected here to fit the experimental semi-variograms of the DEMs.In common with semi-variogram components normally, or perhaps ideally, used in ISM modeling, the Gaussian model provides an "authorized" semi-variogram component in 2D and as a covariance it is also a totally positive kernel with properties that Schoenberg ascribed to "Polya Frequency Functions" [33].
The Gaussian is also convenient for the purpose of deriving analytical results.The model used for a Gaussian component is: The definition of range in Equation ( 9) is consistent with [31] as at this (radial) distance the component semi-variogram reaches close to 95% of the sill value.This definition of range also corresponds to what is called the "FWHM" in the section where the filtering model is introduced.FWHM strictly stands for "Full Width Half Maximum" but is often more practically defined as the integral of a filter divided by its maximum value, and it measures the spatial band width.In deriving these results, it needs to be recalled that C j phq is a radial symmetric 2D function and its total integral over R 2 (using cylindrical symmetry) is: If Gaussian components are used in the ISM model, the change in each component with resolution change can be modeled analytically.

An Upscaling Model for ISM Semi-Variograms and Covariance
The effect of upscaling on the covariance function C phq uses results that can be found in various forms in [14,16,17].However, as the forms vary between these sources and the result used here is not obvious, the basic result is derived in general form in Appendix A in a form suitable for both the present application and others to be used in the future.A general result derived in this paper is suitable for image processing and generalizes the well-known fundamental principles such as change of support and similar analyses used for many years in geology [34] and ecology [35] to analyze scale effects.The specific relationship between these methods and situations has been investigated in [17].In this document, the upscaling model is taken to be the action of a (linear) Low Pass 2D filter ϕ.Change of scale and of grid cell size can then be modeled as applying a pre-filter of this form followed by sampling in which the pre-filter is acting effectively as an anti-aliasing filter.In some cases, the basic filtering is provided by an instrument that collects data rather than by digital filtering of the data.In regards to the situation presented in this paper, the filter is assumed to be linear and radially symmetric and its action modeled as that of 2D image convolution.The covariance function for the upscaled data (C F phq) can be expressed as the convolution of a filter based covariance and the original data covariance (see Appendix A): C F phq " ppϕ ˚Zq ˚pϕ ˚Zqq phq " rpϕ ˚ϕq ˚Cs phq " `Cϕ ˚C˘p hq The process can be illustrated analytically using Gaussian covariance and a Gaussian filter.Let the respective (radial) covariance (C j prq) and radially symmetric filter (ϕ pS; rq) functions be: The term S defines the effective spatial bandwidth of the filter in a similar way that the range R j defines the spatial length scale of the component C j prq.Using the previous result for the 2D integral of radial Gaussian covariance functions, the filter ϕ has (2D) integral 1 and FWHM of ?πS.In addition, as defined above it can take advantage of a convenient analytic convolution property [3] for Gaussian functions that: That is, the convolution of two Gaussians is again a Gaussian [36].It follows that if C j prq is a component of an ISM with variance (sill) σ 2 j and range R j , the result of applying the iterated filter (C ϕ prq) to the covariance function and using the above convolution property for Gaussians is: This is a function of the same type as before but with reduced variance and increased range in each component.The values of the modified (upscaled) semi-variogram component parameters are: It follows from this that the reduction in component variance due to upscaling is smaller if the original range (R j ) of the component is large relative to S (the spatial band width of the filter) and larger if the original range (R j ) of the component is small relative to S. On the other hand, the range ( r R j ) of the filtered data will increase and the increase will be greatest for the high frequency components.Consequently, if the high frequency components only make up a small part of the original variance then the filter will apparently not have a great effect.
Downscaling also involves a linear filter complementary to the upscaling filter.However, it also involves an additional step that will not be covered in this paper.The additional step arises because upscaling loses information.Consequently, complete downscaling involves both filtering and addition of the lost information.It is described in more detail in texts such as [3].

Taking Account of Trends
In the ISM model, each of the components has a range and a sill.The sum therefore also eventually reaches a sill and provides a model that is both stationary in increments and "authorized".The composite sill occurs at a (radial) range where the correlation lengths of the components are all reached and the sill is the sum of the individual model sills.In DEM studies, a data semi-variogram that does not reach a sill is commonly observed due to finite extent of the data.In this paper, it is assumed that large scale regional variations (such as regionally sloping basins or anisotropic mountain chains) create the "trend" as an observed effect in the data so that the situation where a sill is not reached in the extent of the data has been modeled by including components that have correlation lengths (ranges) beyond the extent of the data.It is assumed all of these can be combined into one effect realized as the initial rising section of a composite semi-variogram with range well beyond the data extent.If a Gaussian component with range beyond the extent of the data is used to model the trend within the extent of the data, it can be shown that the separate range and sill parameter values are not individually well defined, but rather a parameter like σ 2 {R is well defined.That is, the way the model rises from the origin to the data extent represents the trend.Therefore, significance will not be ascribed to the individual trend component range and sill values -rather only to the ratio.
Reaching a composite sill within the extent of the data is quite important for some geostatistical methods.It has been established for some time [37] that if the sill is not reached then Ordinary Kriging (OK) is biased or not possible.In this situation, it is traditional to assume that the lack of sill within the extent of the data indicates the presence of a systematic trend.Typically, the trend is estimated as a polynomial and the residuals from a least squares fit are assumed to be random variables and able to be estimated by Kriging using semi-variogams that reach a sill.However, people have also found that the correlation structure of residuals obtained in this way may be hard to model by "authorized" semi-variogram components; and polynomial models have little meaning beyond the extent of the data.The unsatisfactory nature of such polynomial models led Matheron to discuss extensions of geostatistics to functions that are stationary in generalized increments [38] and involve corresponding generalized covariance functions.This approach has reached a high level of sophistication as described by Chilès and Delfiner [14].However, it is rather complex, and simpler methods have also been attempted in practice such as subtracting a moving average mean which includes the Kriging method described by Oliver [27] for estimating ISM components.These still, however, lead to modifications of the covariance structure that are not modeled.
Having a semi-variogram model that does not reach a sill is of no consequence for the upscaling studies undertaken in this paper.However, Kriging will be used to identify scale components in the data and it will affect the use of that tool.Methods introduced in the next section can handle trends and ensure correct identification of components by the Kriging tool as well as provide other useful applications.

The Filter Bank Model for ISM Semi-Variograms and Covariance
The upscaling model has been used in the application of a low pass filter to the finest scale of DEM to reach one of coarser "scale".The results section will demonstrate how this can provide a model for mapping products where "scale" has been chosen by various survey mapping level or data collection methods.A second example uses a method called here the "Filter Bank".It has many application in Remote Sensing analysis [39].The simplest case corresponds to applying a basic high pass filter to the data.It is the result of subtracting an upscaled image of the same resolution from the original image, or: This form of filter is well known in digital image processing [3] to remove low frequency components so it is a "high pass filter".The "I" here is symbolic of the image being untouched by its application.However, even the finest scale (smallest pixel) image can be regarded as an upscaled version of an underlying function with an "S" value (s 0 ) based on its pixel size.Denoting the filter with effective spatial bandwidth s 0 by ϕ 0 , the operation can be written: In general, such operations are examples of "band pass" filters [3].They remove low frequencies defined by the FWHM (S) of ϕ and also remove high frequency effects as defined by the FWHM (s 0 ) of ϕ 0 .From the previous derivations, s ą s 0 .When ϕ 0 ˚f is the data set, this filter corresponds to removing a running "mean" (i.e., low pass filter) from the data.B f can therefore be looked at as the "residual" image after such a running mean is extracted.This type of approach is a common way to remove trends prior to Kriging but, as shown below, it is possible to use the theorem of Appendix A to establish how it affects the covariance structure of the data.
If an image has an ISM model then the corresponding ISM model for the resulting image after applying the filter bank can be derived and for Gaussian components the result is again analytic.Specifically, when subtracting low pass components from data, it follows from the result of Appendix A that: The pairs of convolutions can all be modeled and if the spatial covariance is Gaussian, the result is analytic.In that case, in a similar way as was done for the upscaling and using the Gaussian convolution formula with a single Gaussian ISM component with range R j and sill σ 2 j it follows: This method has also been used in the Results section to identify ISM components in the landscape by modeling the semi-variograms of DEMs and filtered DEMs with different S values together.Thus, the details in a DEM surface will be better identified and modeled.In addition, the filter bank method extinguishes trends and so it also allows Kriging to be used correctly and without the problems it has with trends in the detailed components.

Indicative Goodness of Fit
In this paper, the Indicative Goodness of Fit (IGF) [40] is used to evaluate the result of fitting the semi-variogram models.IGF is defined as follows.
where: N V is the number of semi-variograms included in the fit (indexed by k); n pkq is the number of lags in semi-variogram k (indexed by i); nl ik is the number of pairs used in lag i of semi-variogram k; d k piq is the lag distance for lag i of semi-variogram k; σ 2 k is the variance for the semi-variogram k (not sill but an image estimate for it); ω ik is the combined weight for lag i of semi-variogram k; γ k piq is the semi-variogram for lag i of semi-variogram k; r γ k piq is the estimated (model) semi-variogram for lag i of semi-variogram k.
For a single semi-variogram, the IGF simplifies to: The value for IGF is the RMS error of fit for the semi-variogram as a percentage of the sill value but weighted to give greater importance to fit near the origin and less to long range components.A small value for IGF will indicate a more accurate model.The value of the semi-variogram at the first lag step at finest scale could be considered as noise from the ground.If the IGF value is lower than the percentage value at one lag step (at finest scale) from the origin relative to the sill value for the data, it would mean the model is fitting within the noise level.We will use this criterion to define what we mean by a "good" fit to data.
The IGF provides a rational way to combine data from different semi-variograms into a balanced objective function.For example, semi-variograms based on the finest lag step may be computed to one range and a semi-variogram with larger lag step be computed out to a much farther range.The IGF provides a way to rationally combine them.The same applies for semi-variograms based on the same image at different scales or to semi-variograms from an image and a transformation of the image (including Filter Bank transformations or radial slope and curvature transformations) if the effect on the semi-variogram is known.

Estimating Components of the ISM Model Using Kriging
The ISM model provides a set of semi-variogram component models at different spatial scales.However, it is also important to identify the landscape components corresponding to the ISM components.A method to do this which is suitable for digital images and based on ordinary Kriging is described in [27].This method has been used here in a number of ways to extract estimates of the individual components of the independent components model of Equation (5).When applied directly, the method seems to have some issues (such as semi-variograms of short-range components usually not being the same as or similar to the component semi-variograms) but has often proven useful to help identify the nature of the detailed components of an ISM in the landscape and is taken here for reference to the extension based on a Filter Bank.
This extended method involves Kriging so it needs a semi-variogram model that reaches a sill.It also solves the Kriging equations over subsets of the image so that the matrix inversion is more practicable if only relatively small range components are to be used.The use of subsets also means the processing can be conveniently and efficiently expressed as a set of digital image filters.As the semi-variogram components need to reach sills within the extent of the filters (subsets) long range components pose a problem for all of these requirements.They are handled in the original method [27] by the implicit subtraction of a local mean over the Kriging window in the application.When discussing this later in Appendix C, it will be shown how this is an implicit Filter Bank transformation and as a result the component semi-variograms are changed.In the Results section, it will be clear how the filter bank extinguishes trends and removes very low frequency components.Because the modified semi-variograms are known analytically, it is fully valid to use Kriging to identify the components after application of the Filter Bank.This is demonstrated in the Results section and Appendix C.

Evidence for Upscaling Effects in the Example Selected
Terrain parameters extracted from the selected DEM were influenced by resolution because some detailed scale information was lost as DEM resolution became coarser.For example, Table 1 showed how in the selected example, as resolution became coarser from 2.5 m to 25 m, the max slope value decreased from 70 ˝to 46 ˝, mean slope decreased from 29 ˝to 21 ˝, and STDev of slope decreased from 12 ˝to 8 ˝.Such behavior was well known, as outlined in the Introduction.However, the reason for the decrease of slope as resolution became coarser was that detailed (or finest scale) terrain structures were lost or reduced in magnitude via the upscaling effect.Characterizing terrain by the identified scale components was the aim of the ISM analysis.

2D Anisotropy Analysis
The experimental semi-variograms were computed from the digital numbers of DEMs collected for this study.The semi-variograms were computed to a maximum lag of 4 km (about half of the size of the images).The graphical representations of the semi-variograms in the four principal directions were as shown in Figure 2 which was calculated using the DEM2.5.There was little evidence of anisotropy to lag distances of about 500 m since the semi-variogram from the four directions was similar.The sill variances diverged beyond 500 m.The variation in directions NW-SE and N-S became flat at a lag distance of 3.5 km, whereas variation in directions NE-SW and E-W continued to increase.This could be explained by the direction of the main channel in the area as well as the catchment shape and orientation.Since the differences between variances in different directions did not diverge greatly up to a lag distance of 2 km, the radial semi-variogram was only computed to a lag distance of 2 km in the following as it seemed beyond this point the behavior was dominated by a trend.The effects of anisotropic trends identified here were due to the specific overall catchment shape and orientation and were combined into a radial 'trend" semi-variogram component in the model as it was only the "rotation free" structures within the catchment that was interesting for this study.For the present purposes, a "trend" would be defined as a composite broad scale structure whose semi-variogram component did not reach a sill within lags up to at least half the extent of the DEM. was only the "rotation free" structures within the catchment that was interesting for this study.For the present purposes, a "trend" would be defined as a composite broad scale structure whose semivariogram component did not reach a sill within lags up to at least half the extent of the DEM.

Semi-Variograms for DEMs with Different Resolutions
The radial average experimental semi-variograms calculated based on the selected DEMs were shown in Figure 3.The semi-variograms for these DEMs did not reach a sill before the maximum lag (2 km) which was likely due to the consistent regional structural trends in this area.Radial average semi-variograms of DEM2.5 and DEM25 m were similar in shape.The radial average semi-variogram of DEM2.5 was slightly larger than that of DEM25 which indicatesd that in the DEM25 surface the variance associated with some terrain information was lost compared with DEM2.5.While at first sight the semi-variograms of the DEM might not look very different the upscaling had produced changes in the slopes as shown in Table 1.

Semi-Variograms for DEMs with Different Resolutions
The radial average experimental semi-variograms calculated based on the selected DEMs were shown in Figure 3.The semi-variograms for these DEMs did not reach a sill before the maximum lag (2 km) which was likely due to the consistent regional structural trends in this area.Radial average semi-variograms of DEM2.5 and DEM25 m were similar in shape.The radial average semi-variogram of DEM2.5 was slightly larger than that of DEM25 which indicatesd that in the DEM25 surface the variance associated with some terrain information was lost compared with DEM2.5.While at first sight the semi-variograms of the DEM might not look very different the upscaling had produced changes in the slopes as shown in Table 1.

Semi-Variograms for DEMs with Different Resolutions
The radial average experimental semi-variograms calculated based on the selected DEMs were shown in Figure 3.The semi-variograms for these DEMs did not reach a sill before the maximum lag (2 km) which was likely due to the consistent regional structural trends in this area.Radial average semi-variograms of DEM2.5 and DEM25 m were similar in shape.The radial average semi-variogram of DEM2.5 was slightly larger than that of DEM25 which indicatesd that in the DEM25 surface the variance associated with some terrain information was lost compared with DEM2.5.While at first sight the semi-variograms of the DEM might not look very different the upscaling had produced changes in the slopes as shown in Table 1.

Results for ISM Modeling of the Radial Average Semi-Variograms of DEM
Semi-variograms based on the DEMs with different resolutions were modeled using an ISM in the study area.In order to make the details in the DEM2.5 surface modeled better, a filter bank with band widths of 15 m, 70 m, 200 m and 500 m was used as described in Section 2.7.After a process to determine the smallest number of components needed to fit the data in the finest scale data to an effective level of IGF (see Section 2.8), six components (Y1, Y2, Y3, Y4, Y5 and Y6) were finally used to model the DEMs in the study area.In the following, the range values for the components were numbered as R1, R2, R3, R4, R5 and R6 and the sill values were numbered s1, s2, s3, s4, s5 and s6.The component Y6 represents an arbitrary very long range "trend" component which models broader

Results for ISM Modeling of the Radial Average Semi-Variograms of DEM
Semi-variograms based on the DEMs with different resolutions were modeled using an ISM in the study area.In order to make the details in the DEM2.5 surface modeled better, a filter bank with band widths of 15 m, 70 m, 200 m and 500 m was used as described in Section 2.7.After a process to determine the smallest number of components needed to fit the data in the finest scale data to an effective level of IGF (see Section 2.8), six components (Y1, Y2, Y3, Y4, Y5 and Y6) were finally used to model the DEMs in the study area.In the following, the range values for the components were numbered as R1, R2, R3, R4, R5 and R6 and the sill values were numbered s1, s2, s3, s4, s5 and s6.The component Y6 represents an arbitrary very long range "trend" component which models broader scale features not resolved in the part of the image where the study was focused.Its range value was not critical, only needing to be significantly larger than 2 km, and its ratio of sill to range (its best resolved feature) is consistent for all data sets.When DEMs at a number of resolutions were being compared, what were needed to compute the ISM models using the theory developed in Section 2.5 were the values of S representing the different effects of generalization involved in the various DEMs.The S value for DEM2.5 was fixed to be S 0 " 2.5{ ?π, which was the theoretical minimum value corresponding to the grid cell size.The other S value (in this case there is only one other DEM) was varied as part of the fitting process but constrained so that its value was greater than S 0 .The fitted S value for DEM25 was 28.9 m compared with a theoretical value based on the nominal resolution of 14.1 m indicating the DEM25 is more generalized than implied by the size of the resolution cell.The percentage value at one lag step (at finest scale) from the origin relative to the sill value for the data was 1.54%.The IGF values shown in Table 2 were all smaller than 1.54%, which meant the model was fitting within the noise level.The results of fitting in this way were shown in Figure 4 and Table 3.By using Gaussian models for the semivariogram components as a model for scaling by filtering as described in Section 2.5, an analytic model for the changes in the semi-variogram with scale was derived and then fitted to the data.As a result, the ranges of the components increased and the sills decreased in a direct and predictable relationship with scale generalization.In this process, there was only one underlying model for the ISM and only the parameter S was used to upscale it to match the DEM25.The S values define the relationship between the base model and the model that fits the experimental semi-variograms.Between the upscaled semivariogram models, the sill and range values in the short-range ISM components changed more rapidly and predictably than those in the long-range components as expected from equation (Equation ( 16)).Table 3 showed the range and sill values of each component at resolutions of 2.5 m and 25 m.sill (i.e., variance) of component Y1 decreased by 31.1% as resolution becomes coarser from 2.5 m to 25 m and range increased by 45.2%.The decrease of sill and the increase of range of component Y2 were within 10%.For components Y3, Y4, Y5 and Y6, the decrease of sill and increase of range were very little or not at all.In the upscaling process as the resolution becomes coarser from 2.5 m to 25 m, topographic information was lost.In this case, the loss of the information was mainly up to the scale of about 70 m.As change was almost all associated with component 1 it was possible to infer that it is this change that has brought about the changes in slope distribution observed in Table 1.The results of fitting in this way were shown in Figure 4 and Table 3.By using Gaussian models for the semivariogram components as a model for scaling by filtering as described in Section 2.5, an analytic model for the changes in the semi-variogram with scale was derived and then fitted to the data.As a result, the ranges of the components increased and the sills decreased in a direct and predictable relationship with scale generalization.In this process, there was only one underlying model for the ISM and only the parameter S was used to upscale it to match the DEM25.The S values define the relationship between the base model and the model that fits the experimental semivariograms.Between the upscaled semivariogram models, the sill and range values in the short-range ISM components changed more rapidly and predictably than those in the long-range components as expected from equation (Equation ( 16)).Table 3 showed the range and sill values of each component at resolutions of 2.5 m and 25 m.sill (i.e., variance) of component Y1 decreased by 31.1% as resolution becomes coarser from 2.5 m to 25 m and range increased by 45.2%.The decrease of sill and the increase of range of component Y2 were within 10%.For components Y3, Y4, Y5 and Y6, the decrease of sill and increase of range were very little or not at all.In the upscaling process as the resolution becomes coarser from 2.5 m to 25 m, topographic information was lost.In this case, the loss of the information was mainly up to the scale of about 70 m.As change was almost all associated with component 1 it was possible to infer that it is this change that has brought about the changes in slope distribution observed in Table 1.

Mapping the Topographic Components by Kriging Analysis
Kriging was used in this study to identify the components of the model as described in Section 2.9 and based on the ideas provided in [27].A filter bank with filter band width of 200 m was applied to the DEM and the resulting image used in the Kriging process principally to identify the first two high frequency components (Y1 and Y2).The filter bank was used to improve the modeling of the detailed components as described in more detail in Appendix C. The Y6 component was not included as the range of Y6 was larger than the extent of the area, thus making it unstable in the process of Kriging.The other components were mapped based on the rural DEM surface using Kriging.The trend was also modeled by the Kriging as a local mean estimate but for the filter bank the modeled trend should be effectively zero providing a test for the settings used.
Figure 5 showed an image of the Kriged estimates of the topographic surfaces corresponding to each component at resolutions of 2.5 m and 25 m.The surface for the shortest range component Y1 was significantly affected by resolution.The surfaces derived from longer range Y3, Y4, Y5 did not change very much through the process of upscaling as they were relatively stable low frequency components of the landscape at both resolutions of 2.5 m to 25 m.There was some change in Y2 component, with range changed 8.3% in the upscaling process, however, the differences in Y2 component surface was hard to be recognized visually.

Discussions
(1).The use of geostatistics in a digital image analysis framework In this paper, a fundamental relationship for change of semi-variance with change of scale in gridded digital images has been derived.It has then been applied in the context of the ISM model with the help of the filter bank method to quantitatively describe changes in information as DEM resolution changes.The existence of changes of variance of DEM and variance of slope with scale are well known in practice but using methods involving an analytic ISM model and a general result relating semi-variograms and linear digital filtering, the way they change with resolution has been quantified and proven with support from a simple example.The methods have combined tools from both geostatistics and digital image processing to reach this state.This combination was important as gridded DEM data are a modern data source that involves huge data sets and massive processing loads.The key items derived from traditional geostatistics have been the semi-variogram statistic, which allows data to be stationary in increments, the principles of change of support which are natural in the geostatistical framework and kriging which can be implemented as a linear filter.Kriging was used in this paper in combination with a digital filter bank.
(2).Sources for the DEM data Particular attention needs to be given to the primary sources and pre-processing of the base data.In this paper, the data were generated from a topographic map with software that ensured the grid cell size was small enough to represent the variation in the data and as large as possible to correctly represent the smoothing involved in the original topographic map [29,30].The degree of smoothing in this case has been expressed in the S 0 term for the model.In other cases, gridded DEM data can also be derived from Radar (SAR) data, Lidar data and digital photogrammetry, the effects of processing and degree of implied smoothing as well as the existence of coherent noise has also been investigated [41][42][43].
(3).Relationships between derived terrain parameters and topographic structures An important aim of the analysis of terrain structures is to find the reasons for the changes in a range of terrain factors with resolution and to model the behavior.Among these, slope is one of the most important and practical factors.According to the results of this paper, the properties contributed by short range components will decrease or disappear in the upscaling process, terrain structures will be simpler, and variance from small valleys will decrease.Because slope gradient is affected most by the fine detail components, it also reduces because structures of the long range components dominate in a DEM at coarser resolution.The change of slope as resolution becomes coarser from 2.5 m to 25 m can be observed in Table 1 and is almost all associated with component 1 (range value of about 70 m) as described in the results of this paper.Moving to describe this behavior quantitatively involves application of the general relationship to derivative filters and is the subject of another paper.

Future Work
This paper has provided useful information for understanding the scaling process.However, its major applications are much wider than the demonstration here and there is a great deal of work to be done to apply the full extent of the methods presented.These applications include: (1).Prediction of the terrain with an intermediate spacing In this paper, resolutions of 2.5 m and 25 m were analyzed.As described in Table 2 in the Result section, the S values for 25 m DEM in this study area is 20.7.S values for series of intermediate resolutions could also be detected and modeled.Terrain structures and variances of DEMs and slopes with an intermediate spacing could be obtained assuming similar terrain types (ISM components) and using Kriging with Filter Bank to accommodate trends where they occur.
(2).Upscaling and downscaling topographic parameters Among the topographic parameters, the models for semi-variograms have immediate implications for slope gradients.The variance of slope gradient at each component and the way it changes with resolution can be mathematically derived from the sills and ranges of DEM components.Modeling slope gradient variances at varied resolutions will play a key role in the process of upscaling and downscaling of slope [8].
(3).The application of the methodology to land classification In this paper, terrain structures at varied spatial frequencies were explored.It may be possible to use the ISM components for a more detailed terrain classification.A terrain type is here described on the way the slope or other terrain parameters change with resolution as detected in the former studies by the authors [13].For example, in some places, there may be more variance at high frequency components.In those places, the terrain surface may be more complex than those with less variance in high frequency components.Also, the loss of information as resolution becomes coarser will be more intensive in those areas.Using the method discussed in this paper, we hope to go further with this research into terrain classification.(4).The application of the methodology to DEM data quality assessment The work can be extended to study relative data quality of a certain DEM surface at each component by using the methods discussed in this paper.The filter bank has previously been used for this purpose [39] but not the ISM.If a certain type of DEM with a specified resolution has better or worse data quality (measured by presence or absence of important spatial frequency components) than it is a significant statement of the relative quality.This can improve the assessment and application of different types of DEMs or DEMs based on different sources.

Conclusions
(1) The semi-variogram is a useful tool to characterize topographic spatial structure.However, as there are often different wavelength components in terrain, these components can be usefully modeled by an ISM model.For the data sets used here, the ISM model was able to model the ranges and sills of model semi-variograms to fit the data semi-variograms of various maps with IGF values of about 0.1%; (2) Using general methods describing the way DEM semi-variograms change with scale and an analytic model, it is shown how modeled ranges and sills of the topographic components change quantitatively as DEM scale changes.For example, the sill and range values of the short-range ISM components change more rapidly with upscaling than those of the long-range components.
In the examples provided, the variance of topographic structures that were within 70 m reduced by 31.1%, and topographic structures that were larger than 70 m were affected much less and those above 1500 m not at all.The ISM model and the theory introduced were able to precisely explain this behavior mathematically and measure relative scales of the different DEMs; The ISM components can be identified in the landscape using a Kriging analysis modified by use of the analytic Filter Bank method.By using the filter bank, the short-range components were enhanced and long-range components (including trends) reduced or extinguished.The changes in topographic parameters with resolution such as slope can be explained and quantified as the loss of short-range information in terrain structures at the finest scales.(4) A model based comparison between DEMs with varied resolution clarifies the differences in properties that will be observed when using DEMs in different applications and at different scales.Further, the model for the changes in specified components with resolution also makes it clear how the short range components change with resolution reduction.This can provide estimates for the effects of loss of semi-variogram related information in elevation and other terrain parameters (such as slope) when using coarser resolution DEMs and can also help to choose an appropriate DEM resolution for a particular application.also be found in Brown [44].Covariance and convolution operations are not the same.Convolution is often used in Digital Image Processing due to its having nice mathematical properties, a convenient notation and some convenient but difficult to derive results such as the convolution of two Gaussians that is used in this paper.A general expression for convolution of two functions f and g (where the independent variable t is multivariate (2D) value) is: f `t1 ˘g `t1 ´t˘d ˇˇt 1 ˇŤhe alternative expression involved in the (cross) covariance and the semi-variogram functions is: f `t1 ˘g `t1 `t˘d ˇˇt 1 ˇŤhe relationship between convolution and covariance can be expressed as: cov p f , gq " f ptq ˚g p´tq Among 2D functions, an even function is one for which f ptq " f p´tq and an odd function is one for which f p´tq " ´f ptq.Every function can be expressed as the sum of an even function and an odd function.The Fourier Transformation of an even function is purely real and the covariance and the convolution operations are equivalent.The image semi-variogram can be expressed in terms of the covariance (assuming functions have mean zero and finite variance) as: γ ptq " cov p f , f q p0q ´cov p f , f q ptq The expression cov p f , f q is called (here) the auto-covariance of the function f , and cov p f , gq is called the cross-covariance of f and g.Many texts call these expressions auto-correlation and cross-correlation (e.g., Castleman, 1996 [3]) but it can be confusing when "correlation" can mean normalization by the variance.The Fourier Transforms of the auto-covariance and cross-covariance functions are called the Power Spectrum and the Cross-Power Spectrum.When a random function with mean zero and finite variance has a finite underlying covariance function it is related to the auto-covariance as: The objective is to derive an expression for the underlying covariance of a function after a linear filter has been applied to it.This filtering operation can be expressed as either a convolution or a covariance and we will choose to use covariance.In this case, what we need to prove is: That is, the underlying covariance of the function obtained by filtering the original function is the convolution of the original function's covariance function with the auto-covariance of the filter.The auto-covariance function is always even and real which explains the use of convolution in this expression.
The simplest proof is obtained by using Fourier Transforms.If F psq represents the Fourier Transform of f and writing to represent the action of taking a Fourier Transformation, it follows that: p f ˚gq " F psq G psq p f p´tqq " F psq pcov p f , gqq " F psq G psq where F indicates the complex conjugate of F. That is: The convolution appears because the auto-covariance functions are both even functions and convolution and covariance are equivalent.If the filter ϕ is also an even function then, taking expectations of each side, the result can be written as: This is the form used in the paper to derive the results for upscaling when upscaling is modeled by an even (and radial symmetric) low pass filter and for the Filter Bank when the band pass functions are modeled as differences of differently upscaled DEMs.

Appendix B Illustration of Upscaling Model for ISM Semi-Variograms and Covariance
The effects of upscaling on ISM semi-variograms and covariance can be illustrated using information recorded in Figure A1 and Table A1. Figure A1 plots the data in Table A1.Table A1 lists the parameters for an ISM model with 3 components (Y j ) at a base scale as well as the parameters of the filtered components (Y j 1 ).The filter was chosen to have S " 10m.Table A1 indicates how the sills drop and ranges increase as resolution becomes coarser and the changes in the semi-variogram with resolution, based on Equation ( 16), are illustrated in Figure A1.The behavior is similar to that previously discussed empirically in terms of various terrain types and resolutions by [13].Figure A1 and Table A1 also show how different semi-variogram components will scale differently as modeled by Equation ( 16), with Y1 changing the most and Y3 the least.
ISPRS Int.J. Geo-Inf.2016, 5, 107 17 of 20 , = (, ) That is, the underlying covariance of the function obtained by filtering the original function is the convolution of the original function's covariance function with the auto-covariance of the filter.The auto-covariance function is always even and real which explains the use of convolution in this expression.
The simplest proof is obtained by using Fourier Transforms.If () represents the Fourier Transform of and writing ℑ to represent the action of taking a Fourier Transformation, it follows that: ℑ( * ) = ()() ℑ((−)) =  ̅ () ℑ((, )) = () ̅ () where  ̅ indicates the complex conjugate of .That is: ((, ), (, )) = (, ), (, ) The convolution appears because the auto-covariance functions are both even functions and convolution and covariance are equivalent.If the filter  is also an even function then, taking expectations of each side, the result can be written as: This is the form used in the paper to derive the results for upscaling when upscaling is modeled by an even (and radial symmetric) low pass filter and for the Filter Bank when the band pass functions are modeled as differences of differently upscaled DEMs.

Appendix B. Illustration of Upscaling Model for ISM Semi-Variograms and Covariance
The effects of upscaling on ISM semi-variograms and covariance can be illustrated using information recorded in Figure A1 and Table A1. Figure A1 plots the data in Table A1.Table A1 lists the parameters for an ISM model with 3 components (  ) at a base scale as well as the parameters of the filtered components (  ′ ).The filter was chosen to have S = 10m.Table A1 indicates how the sills drop and ranges increase as resolution becomes coarser and the changes in the semi-variogram with resolution, based on Equation ( 16), are illustrated in Figure A1.The behavior is similar to that previously discussed empirically in terms of various terrain types and resolutions by [13].Figure A1 and Table A1 also show how different semi-variogram components will scale differently as modeled by Equation ( 16), with Y1 changing the most and Y3 the least."Sum" refers to the total modeled semi-variogram; "Yn" refers to the modeled semi-variogram of a component."Sum 1 " (Sum prime) refers to the total upscaled modeled semi-variogram; "Yn 1 " (Yn prime) refers to the modeled semi-variogram of an upscaled component.An S value of 10 can be matched by a change in grid cell size if the filtered data are then sampled at steps of 1 FWHM of the filter or ?πS " 17.17m.This would be the new pixel size after scaling that represents low levels of aliasing and good levels of ability to reconstruct the filtered data at the original pixel size.Conversely, if the pixel size is h, then we can define an equivalent S equal to h{ ?π to be the theoretical maximum S value (Se) that may have been applied to the underlying data before sampling it to that pixel size.That is, it provides a measure of the (minimum) "theoretical grid cell size".

Appendix C Illustration of Importance of Filter Bank in the Process of Kriging in This Paper
The use of a Filter Bank or a designed filter bank pair S 0 and S 1 in the process of Kriging helps to more accurately identify detail components in the landscape.It also leads to more practical use of the Kriging method (smaller filter sizes) and valid handling of trends.Following is an illustration.Figure A2 shows the Component 1 of DEM2.5 used in the paper mapped using Kriging before and after filter bank.The S 0 value equals to 1.41 and S 1 equals to 200.The detailed structures of topographic surface can be better mapped by using filter bank methods in the process of Kriging.
ISPRS Int.J. Geo-Inf.2016, 5, 107 18 of 20 "Sum" refers to the total modeled semi-variogram; "Yn" refers to the modeled semi-variogram of a component."Sum′" (Sum prime) refers to the total upscaled modeled semi-variogram; "Yn′" (Yn prime) refers to the modeled semi-variogram of an upscaled component.An S value of 10 can be matched by a change in grid cell size if the filtered data are then sampled at steps of 1 FWHM of the filter or √ = 17.17.This would be the new pixel size after scaling that represents low levels of aliasing and good levels of ability to reconstruct the filtered data at the original pixel size.Conversely, if the pixel size is h, then we can define an equivalent S equal to ℎ √ ⁄ to be the theoretical maximum S value (Se) that may have been applied to the underlying data before sampling it to that pixel size.That is, it provides a measure of the (minimum) "theoretical grid cell size".

Appendix C. Illustration of Importance of Filter Bank in the Process of Kriging in This Paper
The use of a Filter Bank or a designed filter bank pair S0 and S1 in the process of Kriging helps to more accurately identify detail components in the landscape.It also leads to more practical use of the Kriging method (smaller filter sizes) and valid handling of trends.Following is an illustration.Figure A2 shows the Component 1 of DEM2.5 used in the paper mapped using Kriging before and after filter bank.The S0 value equals to 1.41 and S1 equals to 200.The detailed structures of topographic surface can be better mapped by using filter bank methods in the process of Kriging. Figure A3 shows the modeled and calculated semi-variograms of the highest frequency component before and after applying the filter bank.The two modeled semi-variograms are similar, that is because the most detailed component (The range value is about 70 m) is enhanced by the filter bank used here.The estimated semi-variogram from the original data using Kriging with no filter bank shows a shift in range and a decrease in sill.For the estimated semi-variogram from the original data using Kriging with the filter bank, the sill also decreased, but it shows better properties than that with Figure A3 shows the modeled and calculated semi-variograms of the highest frequency component before and after applying the filter bank.The two modeled semi-variograms are similar, that is because the most detailed component (The range value is about 70 m) is enhanced by the filter bank used here.The estimated semi-variogram from the original data using Kriging with no filter bank shows a shift in range and a decrease in sill.For the estimated semi-variogram from the original data using Kriging with the filter bank, the sill also decreased, but it shows better properties than that with no filter bank.The range is similar with the modeled ones which indicated the Kriging components are more consistent with the modeled components after using the filter bank in the Kriging process.

Figure 1 .
Figure 1.The location and terrain surface of the study area.

Figure 1 .
Figure 1.The location and terrain surface of the study area.

Figure 2 .
Figure 2. Semi-variograms from the four principal directions calculated based on Hc-DEM2.5 m.

Figure 2 .
Figure 2. Semi-variograms from the four principal directions calculated based on Hc-DEM2.5 m.

Figure 4 .
Figure 4. Modeling result for semi-variograms based on the theory of effect of upscaling on semivariograms and Filter bank.(a) DEM2.5;(b) DEM25.

Figure 4 .
Figure 4. Modeling result for semi-variograms based on the theory of effect of upscaling on semi-variograms and Filter bank.(a) DEM2.5;(b) DEM25.

Figure 5 .
Figure 5. Surfaces of components at resolutions of 2.5 m and 25 m Low.

Figure A1 .
Figure A1.Indication of effects of upscaling on semi-variogram components.Figure A1.Indication of effects of upscaling on semi-variogram components.

Figure A1 .
Figure A1.Indication of effects of upscaling on semi-variogram components.Figure A1.Indication of effects of upscaling on semi-variogram components.

Figure A2 .
Figure A2.Highest detail component before and after filter bank.(a) Before filter bank; (b) After filter bank.

Figure A2 .
Figure A2.Highest detail component before and after filter bank.(a) Before filter bank; (b) After filter bank.
ISPRS Int.J. Geo-Inf.2016, 5, 107 19 of 20no filter bank.The range is similar with the modeled ones which indicated the Kriging components are more consistent with the modeled components after using the filter bank in the Kriging process.

Figure A3 .
Figure A3.Semi-variograms of the highest component before and after filter bank.

Table 1 .
Statistics of slope based on DEMs with resolution of 2.5 m and 25 m.

Table 2 .
Modeling results for S and IGF values.

Table 3 .
Parameters of ISM for DEMs based on the theory of effect of upscaling on semi-variograms.

Table A1 .
Parameters for each component of the example.

Table A1 .
Parameters for each component of the example.