Mikhail L Uss
email: uss@xai.edu.ua

Benoit Vozel
email: benoit.vozel@univ-rennes1.fr

Vladimir V Lukin
email: lukin@ai.kharkov.com

Kacem Chehdi
email: kacem.chehdi@univ-rennes1.fr

Estimation of Variance and Spatial Correlation Width for Fine-Scale Measurement Error in Digital Elevation Model

Keywords: Advanced Land Observing Satellite (ALOS) World 3D (AW3D), Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) GDEM2, blind noise parameter estimation (BNPE), digital elevation model

come

I. INTRODUCTION

G RIDDED digital elevation models (DEMs) have found applications for disaster and crisis-management support [START_REF] Voigt | Satellite image analysis for disaster and crisis-management support[END_REF], urban growth monitoring and planning [START_REF] Chaabouni-Chouayakh | Automatic urban area monitoring using digital surface models and shape features[END_REF], and remote sensing (RS) image processing [START_REF] Zhang | DEM-assisted RFM block adjustment of pushbroom nadir viewing HRS imagery[END_REF]. A DEM is subject to error that may impair its quality. This error originates from various sources, including elevation measurement method (traditional optical stereo matching, radar interferometry (IfSAR), or light detection and ranging (LiDAR) [START_REF] Nelson | DEM production methods and sources[END_REF]), instrument, terrain structure, vegetation cover, and DEM interpolation [START_REF] Chen | Robust interpolation of DEMs from LiDAR-derived elevation data[END_REF], making both theoretical and experimental analyses of its properties a complicated task [START_REF] Fisher | Causes and consequences of error in digital elevation models[END_REF]. Only vertical DEM error is of interest in this article. In what follows, the multifactor nature of DEM error is referred as multivariate contextdependence.

DEM error leads to uncertainty in the calculation of terrain attributes such as terrain slope, aspect, or roughness. Even partial knowledge of DEM error properties and its spatial pattern is valuable [START_REF] Fisher | Causes and consequences of error in digital elevation models[END_REF]- [START_REF] Hutchinson | Adding the Z-dimension[END_REF] in such areas as DEM fusion [START_REF] Papasaika | Fusion of digital elevation models from various data sources[END_REF], [START_REF] Schultz | Error detection and DEM fusion using self-consistency[END_REF], filling voids [START_REF] Yue | Highquality seamless DEM generation blending SRTM-1, ASTER GDEM v2 and ICESat/GLAS observations[END_REF], DEM filtering [START_REF] Hutchinson | Adding the Z-dimension[END_REF], [START_REF] Takaku | Adaptive filter for improving quality of ALOS PRISM DSM[END_REF]- [START_REF] Reuter | Preparation of DEMs for geomorphometric analysis[END_REF], DEM interpolation [START_REF] Chen | A robust interpolation method for constructing digital elevation models from remote sensing data[END_REF], modeling of DEM error propagation [START_REF] Temme | Geostatistical simulation and error propagation in geomorphometry[END_REF] in hillslope erosion/failure analysis, land slide risk estimation, and hydrological modeling. Improving knowledge of DEM error was identified as a major research direction in digital terrain modeling domain [START_REF] Wilson | Digital terrain modeling[END_REF].

The basic DEM error characteristic [START_REF] Hutchinson | Adding the Z-dimension[END_REF] is its RMSE = ((1/n) n i=1 (z DEM.iz Ref.i ) 2 ) 1/2 , where z DEM.i is an elevation measurement from the DEM, z Ref.i is the reference elevation measurement of significantly higher accuracy, and n is the number of available measurements. RMSE can be decomposed as RM S E = (σ 2 e + M 2 e ) 1/2 , where mean error or bias M e = (1/n) n i=1 (z DEM.iz Ref.i ) and standard deviation (SD) σ e = ((1/n) n i=1 (z DEM.iz Ref.i -M e ) 2 ) 1/2 [START_REF] Fisher | Causes and consequences of error in digital elevation models[END_REF]. The σ e and M e terms characterize the random and systematic error components, respectively. Spatial properties of DEM error are characterized by its spatial correlation function [START_REF] Temme | Geostatistical simulation and error propagation in geomorphometry[END_REF].

To measure DEM error characteristics, a high-quality reference DEM (e.g., created using LiDAR) or point elevation measurements (obtained by geodetic surveys) are typically utilized [START_REF] Holmes | Error in a USGS 30-meter digital elevation model and its impact on terrain modeling[END_REF], [START_REF] Tadono | Calibration and validation of PRISM onboard ALOS[END_REF], [START_REF] Carabajal | SRTM C-band and ICESat laser altimetry elevation comparisons as a function of tree cover and relief[END_REF]. Comparison between DEMs is a complicated task. Both analyzed and reference DEMs should have the same spatial resolution, cover the same area, and share approximately the same acquisition date. The reference DEM should have significantly better accuracy than the analyzed DEM. Acquisition of such data can be expensive and is not always possible. Direct DEM comparison is further complicated for vegetation cover that can yield different elevation measurements for different instruments (optical, IfSAR, and LiDAR) [START_REF] Toutin | Comparison of stereo-extracted DTM from different highresolution sensors: SPOT-5, EROS-A, IKONOS-II, and QuickBird[END_REF]- [START_REF] Becek | Assessing global digital elevation models using the runway method: The advanced spaceborne thermal emission and reflection radiometer versus the shuttle radar topography mission case[END_REF].

A c c e p t e d M a n u s c r i p t

Abovementioned tight requirements on the reference DEM can be relaxed for featureless flat terrain. For a local DEM patch representing a flat terrain, z Ref = const and DEM SD can be directly accessed as σ e = ((1/n) n i=1 (z DEM.i -(1/n) n i=1 z DEM.i ) 2 ) 1/2 . In this case, reference DEM is not needed to estimate σ e (however, it is still needed to estimate the DEM bias). Having enough flat patches allows collecting many local σ e estimates for DEM accuracy characterization. Importantly, using flat DEM patches, SD of DEM error can be estimated in the no-reference fashion. This approach was leveraged by Becek [START_REF] Becek | Assessing global digital elevation models using the runway method: The advanced spaceborne thermal emission and reflection radiometer versus the shuttle radar topography mission case[END_REF]. He proposed to use runways as flat patches and estimated the bias and SD of the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) global DEM2 [START_REF] Tachikawa | Characteristics of ASTER GDEM version 2[END_REF], [START_REF] Fujisada | Technical methodology for ASTER global DEM[END_REF] (later referred to as GDEM2) and Shuttle Radar Topography Mission (SRTM) DEM [START_REF] Rabus | The shuttle radar topography mission-A new class of digital elevation models acquired by spaceborne radar[END_REF] error.

Usage of flat terrain encounters several problems. First, flat areas are unknown beforehand (with narrow exception exploited by Becek [START_REF] Becek | Assessing global digital elevation models using the runway method: The advanced spaceborne thermal emission and reflection radiometer versus the shuttle radar topography mission case[END_REF]). Second, absolutely flat areas are hard to find. Let us consider the estimation of σ e value locally from a DEM patch of 15 × 15 pixels. For GDEM2, one pixel corresponds to 30 m on the ground; the patch is, thus, 450 × 450 m 2 . The mean σ e value for GDEM2 has been shown in [START_REF] Becek | Assessing global digital elevation models using the runway method: The advanced spaceborne thermal emission and reflection radiometer versus the shuttle radar topography mission case[END_REF] to vary from 2 to 6 m as mean stacking number reduces from 25 to about 5. SD of a DEM patch representing a mixture of terrain elevation with SD σ tr and random noise with SD σ e is σ noisy.terrain = (σ 2 e + σ 2 tr ) 1/2 . To use σ noisy.terrain as an unbiased estimate of σ e , the value of σ tr should be a magnitude smaller than σ e . For example, for σ tr < 0.2σ e , the value of σ noisy.terrain < 1.02σ e . For a low value of σ e , such as σ e = 2 m for GDEM2, σ tr should be of the order of 0.2 • 2 m = 0.4 m. This implies maximum elevation variation ±1.2 m (±3sigma) within the patch. This restriction severely reduces the number of flat patches available for no-reference DEM error characterization, as well as their variability.

The process of estimation of DEM error SD from featureless flat terrain has a direct counterpart in image processing domain, namely, blind noise parameter estimation (BNPE) approach [START_REF] Vozel | Blind determination of noise type for spaceborne and airborne remote sensing[END_REF], where the goal is to estimate image sensor noise characteristics (variance or autocorrelation function) from a mixture of true signal and noise. The common solution is to find and use image homogeneous areas (HAs), where true signal is negligible as compared to sensor noise [START_REF] Van Der Meer | Imaging Spectrometry: Basic Principles and Prospective Applications[END_REF]. The link between no-reference DEM error characterization and BNPE problem is that flat terrain for DEM corresponds to HA for images, DEM error to sensor noise, and true signal to error-free DEM. The problem of automatic HA search is the core of BNPE approach, and many efficient solutions have been proposed. Furthermore, BNPE methodology moved beyond searching HA areas only: moderately heterogeneous areas could be used for noise parameters estimation by assuming spectral differences between noise and the true signal [START_REF] Vozel | Blind determination of noise type for spaceborne and airborne remote sensing[END_REF]. Therefore, BNPE methodology provides suitable approaches to two main barriers to no-reference DEM error characterization discussed above: automatic search of HA (flat DEM patches) and usage of moderately heterogeneous patches (DEM patches corresponding to moderately undulating terrain).

Going through obvious similarities underlined above, adapting key ideas from the BNPE methodology so as to make them fit the issues and problems to solve for gridded DEM is the main contribution of this article. The main difficulty lies in multivariate context-dependence of DEM error compared to the univariate dependence of sensor noise on image intensity widely considered in the BNPE domain. (This issue is discussed more in detail in Section II). To the best of authors' knowledge, the only blind noise parameter estimator designed to deal with the multivariate context-dependent error was proposed by Uss et al. [START_REF] Uss | Analysis of signaldependent sensor noise on JPEG 2000-compressed Sentinel-2 multispectral images[END_REF] and called multivariate, vector estimator of spatially correlated noise using noise informative (NI) areas and fractal Brownian motion (mvcNI+fBm). Therefore, in this article, we modify and evolve mvcNI+fBm estimator to make it adjusted to DEM data and apply it to two global photogrammetric gridded DEMs: GDEM2 and Advanced Land Observing Satellite (ALOS) Global Digital Surface Model "ALOS World 3D" (AW3D) with a 5-m spacing [START_REF]RESTEC Distribute AW3Dtm Digital 3D Topographic Data[END_REF] and 30-m spacing (AW3D30) [START_REF][END_REF], and IfSAR DEM TanDEM-X-DEM [START_REF] Rizzoli | Generation and performance assessment of the global TanDEM-X digital elevation model[END_REF]. In regard to the complexity of DEM error as compared to pure sensor noise, we devote significant efforts to explain and interpret the meaning of estimates obtained with the proposed approach.

The remainder of this article is structured as follows. In Section II, we start by discussing possible sources of error in DEM. Then, we establish that sensitivity to fine-scale error offered by BNPE methodology is fully beneficial to fine-scale elevation measurement error. Section III introduces the mathematical model of the DEM measurement error and several bivariate regression models for the measurement error variance and spatial correlation width. Section IV details the necessary information on the mvcNI+fBm estimator. In Section V, the experimental section, the mvcNI+fBm is used to select the most relevant DEM measurement error model and estimate its coefficients using GDEM2, AW3D30, AW3D, and TanDEM-X-DEM data. It is demonstrated that obtained DEM error models are physically adequate and in good agreement with the respective DEM accuracy analysis published in the literature. On the basis of these promising results, concluding remarks and future work are given.

II. DEM ERROR MODEL

An instrument (sensor) measures elevation at discrete, possibly irregular spaced points on the earth surface. To form a gridded DEM, discrete measurements are interpolated at nodes on a regular grid: Ẑgrid (x t , y s ) = Ẑgrid (t, s), where x t = x 0 + t • r DEM , y s = y 0 + s • r DEM , r DEM is the grid step or DEM spatial resolution, and (x 0 , y 0 ) is the grid origin.

The measured elevation is a coarsened representation of the actual terrain because of the finite instrument spatial resolution (correlation window support in photogrammetric DEM and sounding beam footprint in IfSAR and LiDAR). Elevation measurement process is subject to measurement errors (Fig. 1). For bare soil, the instrument introduces error related to its positioning/orientation accuracy, disparity estimation errors, higher or lower correlation between stereo images, stacking number, epipolar line error in photogrammetric DEM, phase unwrapping errors in IfSAR, and the time-of-delay measurement accuracy in LiDAR [START_REF] Fisher | Causes and consequences of error in digital elevation models[END_REF], [START_REF] Holmes | Error in a USGS 30-meter digital elevation model and its impact on terrain modeling[END_REF], [START_REF] Becek | Assessing global digital elevation models using the runway method: The advanced spaceborne thermal emission and reflection radiometer versus the shuttle radar topography mission case[END_REF], [START_REF] Shortridge | Spatial structure and landscape associations of SRTM error[END_REF], [START_REF] Liu | Investigating DEM error patterns by directional variograms and Fourier analysis[END_REF]. This error source was called instrument-induced in [START_REF] Becek | Assessing global digital elevation models using the runway method: The advanced spaceborne thermal emission and reflection radiometer versus the shuttle radar topography mission case[END_REF]. The bias and SD of this error could depend on terrain morphological attributes such as slope, aspect, and roughness [START_REF] Hutchinson | Adding the Z-dimension[END_REF]. Vegetation canopy is additional source of DEM error [START_REF] Shortridge | Spatial structure and landscape associations of SRTM error[END_REF]. For example, optical sensors measure elevation at the canopy top, LiDAR could measure both the vegetation canopy top (the first return) and ground surface (the last return), and IfSAR measures intermediate elevation between the vegetation canopy top and the ground surface [START_REF] Hutchinson | Adding the Z-dimension[END_REF], [START_REF] Carabajal | SRTM C-band and ICESat laser altimetry elevation comparisons as a function of tree cover and relief[END_REF]. In [START_REF] Nelson | DEM production methods and sources[END_REF], it was demonstrated that both bias and SD of SRTM DEM error significantly increase in areas covered by natural forests as compared to grassland pastures and agricultural areas. The error related to terrain surface properties was called environment-induced in [START_REF] Becek | Assessing global digital elevation models using the runway method: The advanced spaceborne thermal emission and reflection radiometer versus the shuttle radar topography mission case[END_REF]. Finally, interpolation and quantization errors are added. Gridded DEM error may exhibit spatial correlation that should be taken into account as well [START_REF] Fisher | Causes and consequences of error in digital elevation models[END_REF], [START_REF] Holmes | Error in a USGS 30-meter digital elevation model and its impact on terrain modeling[END_REF], [START_REF] Liu | Investigating DEM error patterns by directional variograms and Fourier analysis[END_REF]. In what follows, we refer to factors/attributes influencing DEM error as error predictors and denote them by p 1 , p 2 , . . . , p N p , where N p is number of predictors. In the vector form, predictor vector is defined as p = ( p 1 , p 2 , . . . , p N p ). Dependence of DEM error parameters (SD in particular) on multiple predictors is the essence of multivariate context-dependence.

Therefore, the values of gridded DEM can be represented as Ẑgrid (t, s) = Z grid (t, s) + e grid.mes (t, s, p(t, s)), where Z grid (t, s) is the error-free DEM and e grid.mes (t, s, p(t, s)) is the DEM measurement error. In this article, we treat the difference Z true (t, s) -Z grid (t, s) between the actual elevation Z true (t, s) and error-free DEM not as error, but as uncertainty [START_REF] Fisher | Causes and consequences of error in digital elevation models[END_REF] that depends on the relation between DEM spatial resolution and terrain characteristic size and DEM interpolation method. By DEM error analysis, we understand the study of e grid.mes (t, s, p(t, s)) measurement error term.

To evaluate the DEM error, many researchers have used residuals r (t, s) = Ẑgrid (t, s) -Ẑgrid.ref (t, s) between the analyzed DEM and a more accurate reference DEM. The drawback of this approach is that the reference DEM has its own errors that not always could be neglected [START_REF] Burns | Analysis of elevation changes detected from multi-temporal LiDAR surveys in forested landslide terrain in western Oregon[END_REF]. (t,s))] comprises two error terms. For the first term-difference between error-free DEMs-to be negligible, both DEMs should have the same spatial resolution and perfect spatial alignment (the absence of planimetric errors) [START_REF] Müller | Accuracy assessment of airborne photogrammetrically derived high-resolution digital elevation models in a high mountain environment[END_REF]. The second term includes errors of both analyzed and reference DEMs. To minimize the influence of the reference DEM error, this DEM should be obtained with an order of magnitude accurate instrument free from vegetation canopy and terrain surface parameters influence. Analysis of residuals complicates further if analyzed and reference DEMs were collected with a time lag because of terrain dynamics, vegetation change (leaf-on, leaf-off conditions [START_REF] Burns | Analysis of elevation changes detected from multi-temporal LiDAR surveys in forested landslide terrain in western Oregon[END_REF]), among other factors. The problem remains how to characterize error of an accurate DEM when finding a sufficiently accurate reference DEM is impossible.

In this article, we investigate the possibility of characterizing DEM error without the use of any reference DEM, but relying on noisy measurements Ẑgrid themselves. This problem is known in image processing domain as BNPE problem, and it basically aims at estimating random noise characteristics (variance and autocorrelation function parameters) from a mixture of noise-free image and noise, i.e., noisy image [START_REF] Vozel | Blind determination of noise type for spaceborne and airborne remote sensing[END_REF]. From image processing point of view, the DEM Ẑgrid is a noisy single-component image, Z grid is the noise-free image, and e grid.mes is the noise term. In what follows, we analyze what terms of DEM error are accessible by BNPE. Such a preliminary analysis is essential to correctly interpret the results obtained by the mvcNI+fBm estimator derived for real DEMs.

To illustrate better the connection between DEM error characterization and BNPE problems, let us discuss more in detail the "runway" method proposed by Becek [START_REF] Becek | Assessing global digital elevation models using the runway method: The advanced spaceborne thermal emission and reflection radiometer versus the shuttle radar topography mission case[END_REF], [START_REF] Becek | Investigating error structure of shuttle radar topography mission elevation data product[END_REF]. According to it, if we assume that for a particular DEM patch, the reference DEM is flat, the variance of residual between analyzed and reference DEMs is composed of instrumental (σ 2 I ) and environmental (σ 2 e ) components and undesirable target-induced component σ 2 T that is caused by the unaccounted terrain roughness:

σ 2 Z = σ 2 I + σ 2 e + σ 2 T . The value σ 2
T can be approximated as σ 2 T = (1/12)r 2 DEM tan 2 (α) + (1/12)q 2 , where α is the terrain root mean square slope and q is the quantization step [START_REF] Becek | Assessing global digital elevation models using the runway method: The advanced spaceborne thermal emission and reflection radiometer versus the shuttle radar topography mission case[END_REF]. For a flat terrain (e.g., runway), α ≈ 0 and target-induced term becomes negligible. In the above introduced terms, for flat terrain, Z grid (t, s) = Z grid = const, Ẑgrid. (t,s,p(t,s

))) = Var( Ẑgrid (t, s) -Z grid.ref ) = Var( Ẑgrid (t, s)). ( 1
)
This idea directly corresponds to the so-called HA approach, which is the simplest BNPE approach [START_REF] Van Der Meer | Imaging Spectrometry: Basic Principles and Prospective Applications[END_REF]. The limiting factors of the "runway" method is that: 1) it still relies on reference data to locate flat terrain patches; 2) number of flat patches provided by this method is limited; 3) these patches represent only one particular terrain class-concrete surfaces, and cannot, in principle, be extended to other terrain classes; and 4) (1) requires absolutely flat patches that is hard to satisfy. These limitations have been overcome in advanced BNPE methods that provide means for automatic search of HA and ability to deal with moderately heterogeneous patches. Thus, for global DEMs, BNPE could potentially rely on a larger amount of data for analyzing DEM error and cover larger variety of DEM error predictors.

A c c e p t e d M a n u s c r i p t BNPE has seen fast development in the last decade and now provides a mature set of methods applicable to a variety of scenarios [START_REF] Vozel | Blind determination of noise type for spaceborne and airborne remote sensing[END_REF]: single channel [START_REF] Zoran | Scale invariance and noise in natural images[END_REF], multispectral (or color) [START_REF] Liu | Automatic estimation and removal of noise from a single image[END_REF], and hyperspectral images [START_REF] Uss | Local signaldependent noise variance estimation from hyperspectral textural images[END_REF], [START_REF] Alparone | Signal-dependent noise modelling and estimation of new-generation imaging spectrometers[END_REF]; optical [START_REF] Aiazzi | Unsupervised estimation of signal-dependent CCD camera noise[END_REF] and radar [START_REF] Abramov | Methods for blind evaluation of noise variance in multichannel optical and radar images[END_REF] images; signal-independent (additive) [START_REF] Zoran | Scale invariance and noise in natural images[END_REF], Poisson [START_REF] Foi | Practical Poissonian-Gaussian noise modeling and fitting for single-image rawdata[END_REF], multiplicative noise [START_REF] Abramov | Methods for blind evaluation of noise variance in multichannel optical and radar images[END_REF], or their mixture in the form of general model with signal-dependent noise variance [START_REF] Liu | Automatic estimation and removal of noise from a single image[END_REF], [START_REF] Azzari | Gaussian-Cauchy mixture modeling for robust signal-dependent noise estimation[END_REF]- [START_REF] Colom | Nonparametric multiscale blind estimation of intensity-frequency-dependent noise[END_REF]; methods with ability to characterize only noise variance [START_REF] Zoran | Scale invariance and noise in natural images[END_REF]- [START_REF] Uss | Image informative maps for component-wise estimating parameters of signal-dependent noise[END_REF] as well as noise spatial correlation properties [START_REF] Colom | Nonparametric multiscale blind estimation of intensity-frequency-dependent noise[END_REF], [START_REF] Uss | Maximum likelihood estimation of spatially correlated signal-dependent noise in hyperspectral images[END_REF]. Accuracy of these methods is high enough to consider them as an alternative to direct sensor calibration [START_REF] Liu | Automatic estimation and removal of noise from a single image[END_REF]- [START_REF] Alparone | Signal-dependent noise modelling and estimation of new-generation imaging spectrometers[END_REF], [START_REF] Foi | Practical Poissonian-Gaussian noise modeling and fitting for single-image rawdata[END_REF], [START_REF] Uss | Image informative maps for component-wise estimating parameters of signal-dependent noise[END_REF].

The BNPE problem is ill-conditioned and requires additional a priori information on noise and noise-free image properties to be solved reliably. Such additional information could originate from either spatial [START_REF] Van Der Meer | Imaging Spectrometry: Basic Principles and Prospective Applications[END_REF] or spectral [START_REF] Danielyan | Noise variance estimation in nonlocal transform domain[END_REF] domains leading to two main groups of BNPE methods working in spatial or spectral domain, respectively. Spatial a priori information assumes that noise and noise-free image could be separated in image HA, where noise-free image level variation is negligible as compared to noise. Spectral a priori information assumes that noise-free image and noise difference in spectral domain could help their separation. The stateof-the-art methods typically utilize both spatial and spectral information; the mvcNI+fBm method also belongs to this group. Irrespectively from a priori information used, BNPE methods operate by finding image patches where noise-free image and noise can be separated in the best way and locally estimating noise variance from such patches [START_REF] Vozel | Blind determination of noise type for spaceborne and airborne remote sensing[END_REF].

To illustrate how moderately heterogeneous patches can be used for BNPE problem and which components of DEM error can be estimated by this approach, it is useful to represent BNPE as a high-pass filter. The value of Var(

Ẑgrid (t, s)) = (1/N 2 ) i, j ≤N i, j =1 ( Ẑgrid (t, s) -(1/N 2 ) i, j ≤N i, j =1 Ẑgrid (t, s)) 2 is the local variance of residuals Ẑgrid (t, s) -(1/N 2 ) i, j ≤N i, j =1 Ẑgrid (t, s) between DEM Ẑgrid (t, s) and its low-pass filtered version (1/N 2 ) i, j ≤N i, j =1 Ẑgrid (t, s).
Here, N is the linear patch size. The difference between signal and its low-pass filtered version is equivalent to a high-passed filter applied to the signal. The error-free elevation has multiscale property and possesses fine-, meso-, as well as macroscale-variations [START_REF] Hengl | Mathematical and digital models of the land surface[END_REF]. In contrast, instrument-induced error has slow-varying bias component caused by instrument calibration errors (orientation error and jitter noise) and exhibits random fluctuations at the fine scale. Environmentinduced error is related to the terrain surface and could borrow at some extent its multiscale structure. High-pass filter suppresses constant term and large-scale elevation variations and keeps only fine-scale details. This procedure partly removes error-free DEM variations making some of moderately heterogeneous patches almost homogeneous and suitable for estimating DEM error through BNPE; the error-free DEM and error bias terms are also removed as well as low-frequency calibration error; large-scale environmentinduced variations of DEM error (bias between bare soil DEM and terrain with vegetation cover) are likewise filtered-out. In contrast, fine-scale random variations of DEM error (both instrument-and environment-induced) pass through the high-pass filter. This is illustrated in Fig. 2 where a terrain surface and DEM error are simulated. After high-pass filtering, HA areas represent fine-scale component of DEM error. Heterogeneous areas (shown by gray in Fig. 2) do not represent DEM error, should be detected, and removed from further consideration. In other words, BNPE can be viewed as a filtering process that separates a DEM error from the actual DEM using difference in their autocorrelation functions.

Let us summarize BNPE properties important for understanding effect of their application to gridded DEM: 1) it does not require any reference image (DEM); 2) it characterizes image noise (DEM error) locally; 3) it employs spatial (HAs) and spectral (different autocorrelation functions of the true signal and noise) separability of true signal (error-free DEM) and noise (DEM error); 4) it performs automatic search of image patches where the true signal (error-free DEM) and noise (DEM error) are separable at an extent allowing noise (DEM error) parameters estimation; and 5) it estimates only fine-scale random noise (DEM error) component, and bias and large-scale components cannot be estimated. Therefore, BNPE methods are able to estimate local SD of fine-scale component of DEM measurement error. This error is similar to the residual error or leveled error introduced in [START_REF] Torlegård | A comparative test of photogrammetrically sampled digital elevation models[END_REF] that indicates the best possible DEM error with no instrument orientation error affecting the result.

In the majority of BNPE methods, noise can be represented as the sum of two terms: a signal-independent or additive and signal-dependent one. In RS applications, the additive noise is sensor noise (e.g., thermal noise), and it does not depend on the sensed image. The signal-dependent noise is due to the physical nature of measured quantity (e.g., photon-counting or Poisson noise in optical images, and coherent speckle noise in synthetic aperture radar images), and its variance is a linear or nonlinear function of sensed image intensity [START_REF] Tsin | Statistical calibration of CCD imaging process[END_REF], [START_REF] Aiazzi | A robust method for parameter estimation of signal-dependent noise models in digital images[END_REF]. Therefore, this observation noise model is univariate. In contrast, DEM error depends on many predictors both instrument-and environment-related. To the best of authors' knowledge, the only BNPE method able to deal with the multivariate noise model is the recently proposed mvcNI+fBm (see Section IV for details). This property has determined our choice in favor of this estimator for characterizing the DEM error.

Including all significant predictors selected for analysis in predictor vector p, one expects DEM error to exhibit A c c e p t e d M a n u s c r i p t the same statistical properties for all patches with the same predictor vector. However, not all predictors could be observed, identified, or considered in a particular study. Therefore, predictor vector is split into observable and not-observable, latent, parts: p = (p observed , p latent ). The mvcNI+fBm is designed to estimate measurement error SD σ e for patches with similar observable predictors p observed . However, latent predictors are not controlled leading to variations of the estimates SD: σ 2 e.i = σ 2 e (p observed , p latent.i ), where i = 1 . . . N pt , N pt is the number of patches. In the simplest case, final SD estimates are obtained as average over all patches:

σ 2 e = N -1 pt N pt i=1 σ 2 e.i (p observed , p latent.i ) ≈ σ 2 e (p observed , platent ), where platent = N -1 pt N pt i=1 p latent.i .
The mvcNI+fBm estimates error parameters for a fixed value of observable predictors and an average value of the latent ones.

Significant limitation of BNPE resides in its operation principle. As DEM error SD is estimated from homogeneous or moderately heterogeneous DEM patches representing low relief terrain, predictor vector values that correspond to high relief terrain are inaccessible. Primarily, DEM error parameters for patches with high-slope SD cannot be estimated using the no-reference approach. Homogeneity level required for reliable DEM error SD estimation is related to the DEM error SD itself: the lower the DEM error is, the more homogeneous patches are needed to access it. The DEM error component accessible by the BNPE methods can be finally formulated as fine-scale measurement error SD in low-relief areas averaged over latent, unobservable predictors. In what follows, we refer to this DEM error component simply as measurement error or even error, omitting other qualifiers.

III. ELEVATION MEASUREMENT ERROR MODEL

A. Mathematical Model of DEM Error

Following the discussion in Section II, for an (2N h + 1) × (2N h + 1) DEM patch, we assume the following model:

Ẑgrid (t, s) = Z grid (t, s) + e grid.mes (t, s, p) (2)
where t, s = -N h , . . . , N h denote pixel coordinates within the patch, N h is the patch half size, N = 2N h + 1 is the patch linear size, Ẑgrid (t, s) and Z grid (t, s) are affected by errors and error-free DEMs, and e grid.mes (t, s, p) is a zero-mean spatially correlated random error that models the fine-scale measurement error. For a particular patch, the predictor vector p is assumed constant but changing its value from patch to patch. Model (2) resembles DEM modeling in geostatistics when a DEM patch is split into zero-mean error term and spatially autocorrelated random field modeling error-free DEM [START_REF] Hutchinson | Adding the Z-dimension[END_REF].

It is assumed that error is uncorrelated with the errorfree DEM. By this, we mean that only error parameters are dependent on DEM structure and not error realization itself, i.e., error and error-free elevation values in the same pixel have no linear relationship with each other. The Gaussian distribution is the most natural assumption for DEM error, which is generated by a complex mixture of many factors. The validity of this distribution has been reported in many sources [START_REF] Nelson | DEM production methods and sources[END_REF], [START_REF] Fisher | Causes and consequences of error in digital elevation models[END_REF], [START_REF] Holmes | Error in a USGS 30-meter digital elevation model and its impact on terrain modeling[END_REF], [START_REF] Schultz | Error detection and DEM fusion using self-consistency[END_REF], [START_REF] Liu | Investigating DEM error patterns by directional variograms and Fourier analysis[END_REF].

In addition, the semivariogram analysis of DEM error indicates the presence of spatial autocorrelation [START_REF] Holmes | Error in a USGS 30-meter digital elevation model and its impact on terrain modeling[END_REF], [START_REF] Liu | Investigating DEM error patterns by directional variograms and Fourier analysis[END_REF]. At the lag distances considered in this article (below 100-200 m), isotropy of DEM errors is also typically assumed. (Anisotropy of DEM errors was reported in [START_REF] Liu | Investigating DEM error patterns by directional variograms and Fourier analysis[END_REF] for lag distances exceeding 500 m.) With this, no strong evidence supporting a particular spatial autocorrelation function shape is present in the available literature. We performed experiments with two widely used correlation function shapes: expo

- nential R e (t 1 , s 1 , t 2 , s 2 ) = σ 2 e (p) exp(-(|d|/σ Corr (p))) [16] and Gaussian R e (t 1 , s 1 , t 2 , s 2 ) = σ 2 e (p) exp(-(d 2 /2σ 2 Corr (p))), where d = ((t 2 -s 1 ) 2 + (t 2 -s 1 ) 2 ) 1/2 , σ 2
e (p) is the noise variance, and σ 2

Corr (p) is the parameter that characterizes unknown error spatial correlation width (σ Corr is the distance in pixels where error correlation drops to a certain value, 0.6 for Gaussian shape, and 0.37 for exponential shape). It was found that the Gaussian shape is more adequate than exponential one (in terms of regression quality measured by R-square, see Section V) for GDEM2 and AW3D30 global DEMs, but exponential model is more suitable for TanDEM-X-DEM.

The DEM error properties accepted in the mvcNI+fBm can be summarized as follows: normal distribution, stationary fine-scale error with Gaussian isotropic spatial correlation function, and uncorrelated with the error-free DEM. The goal of the mvcNI+fBm estimator is to evaluate types and coefficients of regression functions σ 2 e (p) and σ 2 Corr (p) using the large number of Ẑgrid patches.

The error-free DEM is modeled as a fractal surface, namely, nonstationary 2-D fBm [START_REF] Pesquet-Popescu | Stochastic fractal models for image processing[END_REF]. The fBm model is suitable for describing the earth surface that was shown to fulfill the multiscale, fractal properties [START_REF] Hengl | Mathematical and digital models of the land surface[END_REF], [START_REF] Barton | Fractals in the Earth Sciences[END_REF], [START_REF] Feder | Fractals[END_REF]. According to the fBm model, the error-free DEM increments Z (t, s) = Z (t, s) -Z (0, 0) with respect to the patch central pixel Z (0, 0) represent a nonstationary random Gaussian process with covariation matrix

R Z (t 1 , s 1 , t 2 , s 2 ) = Z (t 1 , s 1 ) • Z (t 2 , s 2 ) = 0.5σ 2 x t 2 1 + s 2 1 H q + t 2 2 + s 2 2 H q -((t 1 -t 2 ) 2 + (s 1 -s 2 ) 2 ) H q .
where H q ∈ (0,1) is the Hurst exponent describing fBm-field roughness (H q → 0 relates to a rougher terrain, H q → 1 to a smoother one), and σ x describes the fBm amplitude as SD of elevation increments on unit distance. The covariation matrix of measured increments Ẑ (t, s) = Ẑ (t, s) -Ẑ (0, 0) can be shown to have the following form:

R Ẑ (t 1 , s 1 , t 2 , s 2 ) = Ẑ (t 1 , s 1 ) • Ẑ(t 2 , s 2 ) = R Z (t 1 , s 1 , t 2 , s 2 ) + R e (t 1 -t 2 , s 1 -s 2 ) + R e (0, 0) -R e (t 1 , s 1 ) -R e (t 2 , s 2 ). (3)

B. Bivariate Model of Measurement Error Variance and Correlation

Among the multitude of factors influencing the DEM error, in this article, we have considered only two instrument-induced factors valid for photogrammetric DEMs: the stacking number and epipolar line error. e should converge to zero with N stk approaching infinity, which is not realistic due to finite correlation kernel (window) size [START_REF] Uss | A precise lower bound on image subpixel registration accuracy[END_REF]. Therefore, we consider more adequate model σ 2 e = 1+ N -1 stk (in what follows, we represent regression models in simplified form by setting model coefficients to unity).

In stereovision and photogrammetry, finding elevation of a terrain object consists in 1-D search of the object image on the left(right) image along the corresponding epipolar line at the right(left) image. The epipolar line on one stereo pair image is formed by projecting all 3-D points having the same projection in the other stereo pair image [START_REF] Hartley | Multiple view geometry in computer vision[END_REF]. The object displacement (with respect to a reference elevation, typically Z = 0) along the epipolar line is called disparity D, which is related to elevation as Z = D(r/B), where B denotes the base-to-height ratio and r is an instrument spatial resolution [START_REF] Fujisada | Advanced methodology for ASTER DEM generation[END_REF]. In Fig. 3, elevation Z = 0 corresponds to point A and line A-B is the corresponding epipolar line. Object E has higher elevation (and disparity) than object C and is placed farer from point A along the epipolar line. Position of the reference point A and orientation of the epipolar line are determined from the instrument calibration data and are subject to errors. The deviation of the estimated epipolar line from the true one is called epipolar line error [START_REF] Fujisada | Advanced methodology for ASTER DEM generation[END_REF]. It has been shown that camera errors propagate to epipolar line error in a way dependent on disparity [START_REF] Shen | Error propogation from camera motion to epipolar constraint[END_REF]. As a result, SD of distance between points on the true and estimated epipolar lines is increasing with disparity: in average

A -A < C -C < D -D < B -B .
The object disparity is found by maximizing correlation between left and right stereo images within the search window moving along the epipolar line. As estimated epipolar line deviates from its true position, the value of correlation drops, correlation profile maximum smears, and disparity estimation accuracy decreases. In this manner, σ 2 e becomes dependent on D and consequently on elevation Z . We, therefore, select elevation Z as the second predictor of DEM error: p 2 = Z .

Obviously, the influence of epipolar line error depends on the accuracy of an instrument onboard calibration data. In Section V, we show that ASTER GDEM2 with calibration data of relatively low accuracy [START_REF] Fujisada | Advanced methodology for ASTER DEM generation[END_REF] demonstrates the strong influence of elevation on DEM error. On the contrary, for AW3D and AW3D30 DEMs with calibration data of relatively high accuracy [START_REF] Tadono | Precise global DEM generation by ALOS PRISM[END_REF], this dependence is not observed.

In the absence of exact model of epipolar line error influence on σ 2 e , we propose to model it in the form 1 + Z m , where m takes values 1 or 2 in this study (formalizing linear and quadratic dependence hypotheses of measurement error parameters on elevation).

Combining the two univariate terms 1 + N -1 stk and 1 + Z m , we get for σ 2 e a bivariate model in the form (1

+ Z m ) + (1 + Z m )N -1 stk = 1 + N -1 stk + Z m + Z m N -1
stk , leading to the following expression:

σ 2 e (p, θ e ) = σ 2 Var01 +σ 2 Var02 N -1 stk +c VarZ1 Z m + c VarZ2 Z m N -1 stk ( 4
)
where p = (N stk , Z ), N p = 2 is a 2-D predictor vector, and

θ e = (σ 2 Var01 , σ 2 Var02 , c VarZ1 , c VarZ2
) is the DEM error variance parameter vector.

In (4), σ 2 Var01 + c VarZ1 Z m component is shared by all DEM realizations in the stack and related to the underlying surface properties. This error component can be viewed as measurement error of DEM obtained from noise-free stereo images. The term σ 2

Var02 + c VarZ2 Z m is the variance of component that is uncorrelated between DEM realizations in the stack (due to different sensor noise realizations in each stereo pair and different viewing geometries).

For photogrammetric DEM, measurement error spatial correlation width is related to the correlation peak sharpness and should not exceed correlation window size. The value of σ 2 Corr is supposed not to change from one stereo pair to another one. The averaging over stereo pairs does not affect σ 2 Corr of random components and could slightly increase it as a result of spatial smoothing due to horizontal registration error of stereo pairs. Therefore, σ 2 Corr could slightly increase with N stk . As epipolar line error smears correlation peak, σ 2 Corr is expected to increase with elevation.

We propose to use the same model for σ 2 Corr as for

σ 2 e σ 2 Corr (p, θ Corr ) = σ 2 Corr01 + σ 2 Corr02 N -1 stk + c CorrZ1 Z m + c CorrZ2 Z m N -1 stk . (5) where θ Corr = (σ 2 Corr01 , σ 2 Corr02 , c CorrZ1 , c CorrZ2
) is the DEM error correlation width parameter vector.

In addition to full models (4) and ( 5), we have also considered two simpler submodels for both σ 2 e and σ 2 Corr . The first submodel uses only predictor N stk , and the second one employs only Z . Taking into account two possible values for m (1 or 2), we obtain a set of six alternative models to test for the DEM measurement error parameters: 1, 1 + N -1 stk , 1 + Z , 1 + Z 2 , 1 + N -1 stk + Z + Z N -1 stk , and 1

+ N -1 stk + Z 2 + Z 2 N -1 stk .
In what follows, we use the mvcNI+fBm both to select the most relevant model for each considered DEM and estimate the corresponding model coefficients.

IV. NO-REFERENCE ESTIMATION OF MEASUREMENT ERROR PARAMETERS: THE MVCNI+FBM ESTIMATOR

A. Generalized Description of the Adapted mvcNI+fBm Estimator of DEM Error Parameters

The main processing stages of the mvcNI+fBm method include: 1) initialization; 2) patches homogeneity estimation; 3) grouping of image patches; 4) processing of the formed groups; and 5) multivariate estimation of context-dependent error parameters (see Fig. 4) [START_REF] Uss | Analysis of signaldependent sensor noise on JPEG 2000-compressed Sentinel-2 multispectral images[END_REF]. In the mvcNI+fBm method, homogeneity estimation stage depends on DEM error parameters; therefore, iterative repeating of stages (2)-( 5) is required. Originally, the mvcNI+fBm was proposed for vector RS images. In this article, it is adapted to the case of gridded DEM images that represents single-channel data. Let us consider each stage of the mvcNI+fBm in more detail.

To estimate each error parameter, one specific call of sequence of stages (2)-( 5) is performed. During each call, one parameter is fixed while the other ones are estimated, and vice-versa. For the purpose of clarity, we describe below the proposed estimation procedure for an arbitrary error parameter s, meaning either noise SD, s = σ e , or spatial correlation width, s = σ Corr .

At the initialization stage, the image to process is split into nonoverlapping patches. Each patch represents a group of adjacent pixels of size N × N pixels. The patches with known unreliable data are removed from further consideration. This stage is application-specific and is described in Section V. Overall number of retained patches is referred to as N pt .

In the mvcNI+fBm, homogeneity index is measured as a ratio between lower bound on parameter s estimation SD (Cramér-Rao lower bound, CRLB), σ s , and parameter value: r HA = σ s /s. For homogeneous patches, it takes low values. For example, for SD parameter, very homogeneous patches represent pure DEM error, σ σ 2 e approaches to SD of sample variance: σ σ 2 e ≈ σ e (2/(N 2 -1)) 1/2 , and r HA ≈ (2/(N 2 -1)) 1/2 . For N = 11 used in Section V, (2/(N 2 -1)) 1/2 = 0.129. For heterogeneous patches, where actual DEM and DEM errors cannot be separated, both σ s and r HA ratio increase. In extreme cases, for very rough patches r HA → ∞. The value of r HA indicates how much information about error parameter s can be extracted from a given patch with patch homogeneity increasing as r HA decreasing. Correspondingly, in the mvcNI+fBm, HAs with low r HA index were called NI. Homogeneity estimation stage involves the estimation of error-free DEM parameters, value of σ s , and, finally, the value of r HA for each patch (see Section IV-B for details). Notice that a specific value of r HA is defined for each considered error parameter. Therefore, homogeneity estimation stage (and subsequent stages of mvcNI+fBm) is repeated for each error model parameter.

Let us consider two patches with the same error parameter s and CRLBs σ s1 and σ s2 . The two patches considered together provide more information on s than each of them considered individually. By joint processing of two patches, CRLB value could be potentially reduced to σ s = σ s1 σ s2 /(σ s1 + σ s2 ) < min(σ s2 , σ s1 ). Therefore, for two patches, r HA < min(r HA1 , r HA2 ). From this point of view, several patches considered together are "more homogeneous" and more "NI" than each of them. For a group of N group patches with similar DEM properties, homogeneity index decreases approximately as 1/(N group ) 1/2 . The mvcNI+fBm utilizes this property by grouping homogeneous and moderately heterogeneous patches with similar predictor values into even more homogeneous groups. Each group is further processed as a whole. In this manner, even moderately heterogeneous patches could form a homogeneous group, thus relaxing requirement to actual DEM variability and increasing number of individual error parameter estimates.

The grouping stage splits all available patches into cells with similar predictor values (implying similar error properties). Similarity of predictors is measured as

d(p 1 , p 2 ) = ( N p j =1 w 2 j • (p 1 ( j ) -p 2 ( j )) 2 ) 1/2
, where w j is the associated weight. Weights w j , j = 1, . . . , N p , serve two interdependent goals: the first one is to limit the error parameter variation within the group and the second one is to assure the desired resolution with respect to predictors. To provide resolution A c c e p t e d M a n u s c r i p t IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING p j with respect to j th predictor, the weight should be set as w j = p -1 j (the used weight values for the DEM error parameter estimation are specified in Section V). In each cell, patches are sorted in the ascending order of r HA value (from the most homogeneous to the most heterogeneous ones). Patches are sequentially assigned to groups. The groups are considered homogeneous when r HA drops below the threshold value of 0.15 and N group < 15. (A larger size of a homogeneous group means that the set of patches in the group is more heterogeneous and consequently extracting error parameters from such a group may become less reliable; in the tests carried out later in Section V, we found 15 to be a reasonable group size.) We define by N NI the number of the found NI groups. For the i th group, mean value of predictor vector pi is calculated.

At the next stage (4) of processing the formed groups, for an i th group, a single error parameter s i (σ 2 e.i or σ 2 Corr.i ) is estimated using all patches in this group. Thus, each group provides one observation ŝi of the error parameter s i for a given value of predictor vector pi . Note that the lower bound on ŝi estimation error, denoted by σ 2 s i , is derived as an additional output of the estimation procedure (see Section III-B).

The final stage (5) of estimation of multivariate contextdependent error parameter s mathematically corresponds to the associated robust, multivariate, heteroscedastic regression problem so as to estimate s as a function of its predictors [START_REF] Uss | Analysis of signaldependent sensor noise on JPEG 2000-compressed Sentinel-2 multispectral images[END_REF]. Multivariate linear regression (nonlinear regression can be applied in a similar way) is, thus, used in the mvcNI+fBm rather than simple univariate linear regression typically considered in the past for signal-dependent noise parameter estimation. Robustness is needed to cope with gross DEM errors (outliers and blunders) [START_REF] Fisher | Causes and consequences of error in digital elevation models[END_REF].

B. Mathematical Details of the mvcNI+fBm Estimator

Let us formally introduce the mvcNI+fBm estimator. Using the definitions of Section III-A, the log-likelihood function of image intensities within a single patch is given by ln

L( Ẑ; θ) = - 1 2 ẐT R -1 Ẑ Ẑ + ln |R Ẑ| ( 6
)
where Ẑ is an N 2 × 1 sample composed of all Ẑ (t, s) within the patch, kth Ẑ element and patch pixel with coordinates (t k , s k ) are related to each other by k = t k + Ns k + N h (N + 1) + 1, and θ = (σ 2

x , H q , σ 2 e , σ 2 Corr ) denotes the full parameter vector characterizing both error-free DEM and DEM error.

The information contained in the sample Ẑ is characterized by the Fisher information matrix (FIM) I θ with elements 4 [61].

I θ i θ j = (1/2)tr(R -1 Ẑ(∂ R Ẑ/∂ θ i )R -1 Ẑ(∂ R Ẑ/∂ θ j )), i, j = 1, . . . ,
Before estimating error parameters, we first need to determine error-free DEM parameters, σ 2

x and H q , using the following estimator:

σ 2 x , Ĥq = arg min σ 2 x ≥0,0<H q <1,σ 2 e = σ 2 e ,σ 2 Corr = σ 2 Corr [ln L( Ẑ; θ )]. (7)
In [START_REF] Holmes | Error in a USGS 30-meter digital elevation model and its impact on terrain modeling[END_REF], we specify that error parameters are fixed at this stage at their current estimated values σ 2 e and σ 2 Corr . Error-free DEM parameter estimates are only reliable for patches with ratio σ 2

x /σ 2 e that can be understood as local signal-to-noise ratio (SNR) within the given DEM patch. High SNR patches with σ 2

x /σ 2 e > 2 are called texture informative (TI). For other patches, Hurst exponent estimation is not reliable and is evaluated using the inverse distance interpolation method [START_REF] Shepard | A two-dimensional interpolation function for irregularlyspaced data[END_REF] based on reliable estimates Ĥq from TI patches. In this manner, the interpolated Hurst exponent values Ĥq.interp are obtained for each patch. To characterize interpolation error, we sequentially exclude one TI patch and interpolate Hurst exponent value for this patch using remaining TI patches. The interpolation error SD σH q.interp is calculated as SD of difference Ĥq.interp -Ĥq for all TI patches.

Two estimators for the error variance σ 2 e and spatial correlation width σ 2 Corr can now be defined in the following form:

σ 2 x , Ĥq , σ 2 e = arg min σ 2 x ≥0,0<H q <1,σ 2 e ≥0,σ 2 Corr = σ 2 Corr × ln L( Ẑ; θ ) - (H q -Ĥq.interp ) 2 2 σ 2 H q.interp (8) σ 2 x , Ĥq , σ 2 Corr = arg min σ 2 x ≥0,0<H q <1,σ 2 e = σ 2 e ,σ 2 Corr ≥0 × ln L( Ẑ; θ ) - (H q -Ĥq.interp ) 2 2 σ 2 H q.interp . ( 9
)
In ( 8) and ( 9), we have taken into account the prior on the Hurst exponent as the additional term in the optimized loglikelihood function. This prior has maximal effect for the very homogeneous patches with low SNR, for which reliable Hurst exponent estimation is impossible. The effect becomes less pronounced with SNR increase.

FIMs I (σ 2

x ,H q ,σ 2 e ) and I (σ 2

x ,H q ,σ 2

Corr ) for the estimators ( 8) and ( 9), respectively, are obtained from general matrix I θ by eliminating row and column corresponding to the missing parameter and adding prior information on H q by substituting I H q H q component with I H q H q + σ -2 H q.interp . CRLBs on σ 2 e and σ 2 Corr estimate error are obtained by inverting I (σ 2

x ,H q ,σ 2 e ) and I (σ 2

x ,H q ,σ 2 Corr ) . The last elements of these matrices, σ 2 Taking into account error-free DEM and error independence between patches, estimators for the group of N group patches can be straightforwardly formulated using [START_REF] Hutchinson | Adding the Z-dimension[END_REF] and [START_REF] Papasaika | Fusion of digital elevation models from various data sources[END_REF] Dx , Ĥq , σ 2 e = arg min

D x ≥0,0<H q <1,σ 2 e ≥0,σ 2 Corr = σ 2 Corr × ⎡ ⎣ N group v=1 ln L( Ẑv ; θ v ) - (H q.v -Ĥq.interp.v ) 2 2 σ 2 H q.interp ⎤ ⎦ (10)
D x ≥0,0<H q <1,σ 2 e = σ 2 e ,σ 2 Corr ≥0 × ⎡ ⎣ N group v=1 ln L( Ẑv ; θ v ) - (H q.v -Ĥq.interp.v ) 2 2 σ 2 H q.interp ⎤ ⎦ ( 11
)
where Ẑv is the sample for vth patch,

θ v = (σ 2 x.v , H q.v , σ 2
e , σ 2 Corr ), σ 2 x.v , H q.v , and Ĥq.interp.v are σ 2

x , H q , and Ĥq.interp for vth patch, D x = (σ 2

x.1 , σ 2 x.v , . . . , σ 2

x.N gr ), and H q = (H q.1 , H q.2 , . . . , H q.N gr ). Note that the same error parameters are assumed for all patches in the group.

For nonoverlapping patches, error parameter estimates are independent. Therefore, CRLBs on σ 2 e and σ 2 Corr using a group of patches are given by σ 2

σ 2 e = 1/ N group v=1 (σ -2 σ 2 e.v
) and σ 2

σ 2 Corr = 1/ N group v=1 (σ -2 σ 2 Corr.v
), respectively.

Let us provide R Ẑ derivatives with respect to

θ i , i = 1, . . . , 4 ∂R Ẑ(k 1 , k 2 ) ∂σ 2 x = σ -2 x R Z (t k 1 , s k 1 , t k 2 , s k 2 ) (12) ∂R Ẑ(k 1 , k 2 ) ∂ H q = 0.5σ 2 x log t 2 k 1 + s 2 k 1 t 2 k 1 + s 2 k 1 H q + log t 2 k 2 + s 2 k 2 t 2 k 2 + s 2 k 2 H q -log((t k 1 -t k 2 ) 2 + (s k 1 -s k 2 ) 2 ) ×((t k 1 -t k 2 ) 2 + (s k 1 -s k 2 ) 2 ) H q (13) ∂R Ẑ(k 1 , k 2 ) ∂σ 2 e = ⎛ ⎜ ⎜ ⎝ 1 + exp - d 2 0 2σ 2 Corr -exp - d 2 1 2σ 2 Corr -exp - d 2 2 2σ 2 Corr ⎞ ⎟ ⎟ ⎠ ( 14
)
∂R Ẑ(k 1 , k 2 ) ∂σ 2 Corr = σ 2 e 2σ 4 Corr ⎛ ⎜ ⎜ ⎝ d 2 0 exp - d 2 0 2σ 2 Corr -d 2 1 exp - d 2 1 2σ 2 Corr -d 2 2 exp - d 2 2 2σ 2 Corr ⎞ ⎟ ⎟ ⎠ (15)
where

d 0 = ((t k 1 -t k 2 ) 2 + (s k 1 -s k 2 ) 2 ) 1/2 and d 1 = (t 2 k 1 + s 2 k 1 ) 1/2 , d 2 = (t 2 k 2 + s 2 k 2 ) 1/2 .
Finally, DEM measurement error parameters estimation algorithm is detailed in Algorithm 1.

V. APPLICATION OF ADAPTED MVCNI+FBM ESTIMATOR

TO REAL DEMS ERROR ANALYSIS In this section, we apply the proposed mvcNI+fBm method to two global DEMs with 1 arc second (30 m at equator) spatial resolution: GDEM2 and AW3D30. More attention is

Algorithm 1 Estimation of DEM Measurement Error Parameters

Input: A number of a raster DEM tiles; Output: DEM context-dependent error parameters (error variance σ 2 e and spatial correlation width σ 2 Corr dependence on predictor vector); 1. Split DEM tiles into non-overlapping patches; Reject unreliable patches (as described in subsection 5.1); 2. Set iteration number k = 1; 3. For each of N pt selected patches, calculate predictor vector (e.g. stacking number, elevation) p i , i = 1 . . . N pt ; 4. For each patch, calculate fBm field parameters vector according to [START_REF] Holmes | Error in a USGS 30-meter digital elevation model and its impact on terrain modeling[END_REF]. Interpolate Hurst exponent for low SNR patches (NI) using high SNR patches (TI); 5. For each patch, calculate Cramer-Rao lower bounds on σ 2 e or σ 2 Corr estimate error using ( 12) -( 15); Calculate patch homogeneity index; 6. Group patches into NI groups using homogeneity index (as described in subsection 4.1); For each NI group calculate mean value of predictor vector p j , j = 1 . . . N NI , N NI is number of NI groups; 7. For each NI group, estimate error parameter s j using [START_REF] Schultz | Error detection and DEM fusion using self-consistency[END_REF] or ( 11 paid to GDEM2 as it reveals the most complex measurement error behavior. To prove the mvcNI+fBm scalability, we test it on a sample granule of AW3D DEM with a 5-m spatial resolution. In addition, radar DEM TanDEM-X-DEM is analyzed to compare mvcNI+fBm results with theoretically predicted random noise SD. The obtained results for each DEM are compared with respective DEM accuracy analysis available in the literature.

A. Experimental Settings

For GDEM2 and AW3D30 DEMs, 59 tiles were used to estimate σ 2 e and σ 2 Corr model coefficients with the following lower left (southwest) corner pixels: N27-28E086-087, N30-32E035, N33-35E076, N47-48W001-002, N48-49E001-003, N48-49E031-032, N49-50E036-037, N50E099, N50E103, N51-52E099, N52E101-102, N53E103, and S23-28W067-070. Here, hyphen defines range of tiles with respect to latitude, longitude, or both; N and S denote the north and south latitudes, respectively; W and E denote the west and east longitudes, respectively.

Each tile has fixed latitudinal and longitudinal spatial resolution of 1 arc second (approximately 30 m at the equator) but expressed in meters, and resolution with respect to image rows and columns differs. To compensate for this effect, tiles were resampled to coarser 90-m resolution cell with respect to both spatial coordinates. For each tile, 10 000 nonoverlapping patches of 11 × 11 pixels were selected on a regular grid. The selected size of the patch is an experimentally found compromise between accuracy of error parameters estimates (that increase with patch size) and DEM patch homogeneity (that decreases with patch size). The stacking number for each pixel of these patches was obtained using quality assessment (QA) files supplied for each tile (the QA file has the same number of rows and columns as DEM tile with each pixel representing the corresponding N stk value. QA file has extension.num for GDEM2 and suffix "STK" for AW3D). QA files were also interpolated to 90-m grid to match DEM tiles.

We set two criteria to decide whether a patch is reliable for measurement error parameter estimation: 1) it contains data provided by the "normal" workflow for the particular DEM and 2) patch is not significantly affected by GDEM2 stacking procedure artifacts (see the next paragraph). The first criterion for photogrammetric DEM excludes pixels where disparity measurement procedure fails, primary from lack of data (e.g., caused by clouds), water body (that are typically masked out), and low correlation areas (e.g., deserts). Those pixels are substituted from different sources (e.g., other DEMs) or remain as voids and do not characterize measurement error. Such pixels can be identified by a negative stacking number for GDEM2 and from mask information file (marked by "MSK" suffix) for AW3D, and from water indication mask ("WAM") files of Tan-DEM-X DEM.

The second criterion is for handling another problem causing gross elevation error is related to GDEM2 stacking procedure: it produces false elevation discontinuities if stacking number changes severely. This problem is illustrated in Fig. 5, where a GDEM2 DEM patch, the corresponding stacking number map, and the AW3D30 patch from the same location are shown. The false elevation discontinuity correlated with the stacking number discontinuity is clearly visible, while at AW3D30 patch, this feature is missing. To avoid this source of gross errors, the minimum N stk.min and maximum N stk.max numbers of stereo pairs were calculated for each patch. The patch was rejected if N stk.max /N stk.min > 2, thus limiting stacking number variability. All patches that were found reliable are further processed by the mvcNI+fBm estimator.

Predictor weights were set as w 1 = 1 and w 2 = 0.01 m -1 for N stk and Z , respectively. Required homogeneity index for NI groups was experimentally set to 0.125. We found that in all cases, mvcNI+fBm converged in less than 15 iterations. Computational complexity of mvcNI+fBm is moderate: single iteration of mvcNI+fBm for 59 tiles takes about 5 h on Intel Core i5-7200U CPU. The algorithm can be significantly sped up by parallel processing of groups of fragments on GPU.

We use classical tools of linear regression to characterize regression models: generalized determination coefficient, partial residual plot, and significance test for predictor values using t-statistic [START_REF] Seber | Linear Regression Analysis[END_REF], [START_REF] Altman | Statistics with Confidence: Confidence Intervals and Statistical Guidelines[END_REF]. The generalized determination coefficient [START_REF] Magee | R 2 measures based on wald and likelihood ratio joint significance tests[END_REF]: R 2 = 1exp(-(2/N NI )(l Rl U )), where l R is the log-likelihood for regression model restricted only with intercept, l u is the log-likelihood of the unrestricted regression model, and N NI is the number of estimates (number of NI groups). Assuming normal distribution of error parameter estimate errors provided by the mvcNI+fBm and omitting constants shared by both l R and l u , the log-likelihoods for restricted and unrestricted models take the following quite simple form:

l R = -(1/2) N NI i=1 ((ŝ i -s const ) 2 /σ 2 s i ) and l u = -(1/2) N NI i=1 ((ŝ i -s pr.i ) 2 /σ 2 s i ),
where s pr.i is the parameter prediction for i th measurement.

B. ASTER GDEM2 Experiment

The ASTER global DEM version 2 was released by the National Aeronautic and Space Administration (NASA), Washington, DC, USA, and the Ministry of Economy, Trade, and Industry (METI) of Japan in 2011 [START_REF] Tachikawa | Characteristics of ASTER GDEM version 2[END_REF], [START_REF] Fujisada | Technical methodology for ASTER global DEM[END_REF]. Experimental assessment of the GDEM2 has been published in numerous articles including [START_REF] Becek | Assessing global digital elevation models using the runway method: The advanced spaceborne thermal emission and reflection radiometer versus the shuttle radar topography mission case[END_REF], [START_REF] Tachikawa | Characteristics of ASTER GDEM version 2[END_REF], [START_REF] Goncalves | Accuracy analysis of DEMs derived from ASTER imagery[END_REF] among others.

The ASTER GDEM2 was produced from data collected by the ASTER sensor on board of the NASA Terra spacecraft, which is capable of collecting in-track stereo pairs using nadir-and aft-looking near-infrared cameras [START_REF] Fujisada | ASTER stereo system performance[END_REF], [START_REF] Iwasaki | ASTER geometric performance[END_REF]. The procedure for the GDEM2 generation is described in [START_REF] Fujisada | Advanced methodology for ASTER DEM generation[END_REF]. According to this procedure, elevation estimates are obtained via correlation-based registration with 5 × 5 pixel window. Curve fitting approach [START_REF] Debella-Gilo | Sub-pixel precision image matching for measuring surface displacements on mass movements using normalized cross-correlation[END_REF] is used to reach the subpixel registration accuracy. The search is performed in the along-track direction neglecting a small cross-track shift component appearing due to the contribution of earth rotation effects during the stereoscopic observation period. Multiple stereo pairs (stacking number up to 50) are used to obtain elevation estimates for a given area.

1) Measurement Error Variance Analysis:

The results obtained for error parameter σ 2 e are presented in Table I. It was found that models with quadratic dependence on elevation are more significant according to R 2 coefficient (R 2 is 0.4335 and 0.6298 for models 1 + Z and 1 + Z 2 ; 0.7717 and 0.8954 for models 1 + N -1 stk + Z + Z N -1 stk and 1 + N -1 stk + Z 2 + Z 2 N -1 stk , respectively). Using only N stk predictor, we get a low determination coefficient R 2 = 0.1221 (model 1 + N -1 stk ). The second Z predictor notably increases R 2 to 0.6298 (model 1 + Z 2 ). Using all terms in the model (4) increases R 2 even further to the value of 0.8954. Therefore, the influence of both predictors on σ 2 e is significant. The estimated model for the GDEM2 elevation measurement error variance is given as follows: ). Let us validate further the selected model ( 16) both qualitatively and quantitatively in comparison with GDEM2 analysis results provided in the available literature. For the ASTER instrument, B = 0.6 and r = 15 m. The ratio (r/B) = 15 m/pixel is conversion ratio between pixel and meter units. It is systematically used in the analysis below where conversion between measurement units is involved.

The value of σ e at the sea level (Z = 0) and for one stereopair (N stk = 1) is given as (σ 2

Var01 + σ 2 Var02 ) 1/2 = 5.1668 m. This value was estimated by Fujisada et al. [START_REF] Fujisada | Advanced methodology for ASTER DEM generation[END_REF] (see [START_REF] Fujisada | Advanced methodology for ASTER DEM generation[END_REF]Fig. 10], kernel size 5 × 5 pixels) to vary from 0.2 to 0.3 pixels or, equivalently, from 5 to 7.5 m, depending on the terrain. One can conclude and underline good agreement between these two estimates.

The value of σ e introduced in this article corresponds to the SD of elevation bias in [START_REF] Becek | Assessing global digital elevation models using the runway method: The advanced spaceborne thermal emission and reflection radiometer versus the shuttle radar topography mission case[END_REF]. Therefore, we can directly compare model [START_REF] Temme | Geostatistical simulation and error propagation in geomorphometry[END_REF] with the results by Becek [START_REF] Becek | Assessing global digital elevation models using the runway method: The advanced spaceborne thermal emission and reflection radiometer versus the shuttle radar topography mission case[END_REF]. The reported results of the dependence of SD of bias on N stk make this comparison even more informative. The majority of 96 runways considered in [START_REF] Becek | Assessing global digital elevation models using the runway method: The advanced spaceborne thermal emission and reflection radiometer versus the shuttle radar topography mission case[END_REF] are situated at low altitudes near the sea level with the mean elevation of 234.7 m. We simplify model ( 16) by substituting Z = 234.7 m getting the reduced model σ 2 e = 1.0563 + 26.0031N -1 stk . The results obtained with the reduced model were compared to those shown in [START_REF] Barreiro-Fernández | Accuracy assessment of LiDAR-derived digital elevation models in a rural landscape with complex terrain[END_REF]Fig. 10]. Both estimates are shown in Fig. 7 for N stk varying from 0 to 50.

The results of the mvcNI+fBm and runway methods are highly consistent for N stk > 20. For N stk < 20, σ 2 e estimates obtained by the mvcNI+fBm are up to two times smaller than those provided by the runway method. The possible reason of this difference is the better robustness of mvcNI+fBm to gross GDEM2 errors and outliers that affect the GDEM2 for low stacking numbers as indicated by Becek [START_REF] Becek | Assessing global digital elevation models using the runway method: The advanced spaceborne thermal emission and reflection radiometer versus the shuttle radar topography mission case[END_REF].

The estimates of c VarZ1 and c VarZ2 related to epipolar line error show very significant influence of elevation on σ 2 e : σ e at Z = 6000 m increases by 1.9…3.4 times as compared to this value at Z = 0 m making elevation influence the predominant error source at high elevations. The epipolar line error SD is measured in [START_REF] Fujisada | Advanced methodology for ASTER DEM generation[END_REF] to be about 0.15 pixels or 3.75 m for elevation about 4000 m (Mount Elbert). According to model ( 16), σ 2 e at elevation Z = 4000 m and typical number of stereo pairs N stk = 10 increase by 4.899 • 10 -7 (4000 m) 2 + 6.107 • 10 -6 (4000 m) 2 /10 = 17.61 m 2 . This additional SD of (17.61 m 2 ) 1/2 = 4.2 m is very close to the value reported in [START_REF] Fujisada | Advanced methodology for ASTER DEM generation[END_REF].

2) Measurement Error Correlation Width Analysis: For the correlation width parameter (Table I), significant dependence on elevation can be observed (determination coefficient is about 0.85-0.9 for the models 1 + Z and 1 + Z 2 ) revealing epipolar line error influence on measurement error correlation function. Similar to what occurs for modeling σ 2 e , models depending on Z are less relevant than models depending on Z 2 for modeling σ 2

Corr . Joint usage of two predictors, N stk and Z , cannot definitely confirm or discard hypothesis of N stk influence on correlation width: t-stats for terms dependent on N stk have rather small values (less than 10; see partial regression plots in Fig. 6). Therefore, we consider model 1 + Z 2 as the most relevant one for the GDEM2 measurement error spatial correlation width (numeric values are in 90-m pixels)

σ 2 Corr = 0.1937 + 1.7786 • 10 -8 Z 2 . ( 17
)
The constant term of correlation width estimate obtained by the mvcNI+fBm method for Z = 0 is about 0.44 pixels at the resolution of 90 m that corresponds to 2.

C. ALOS Experiment (5-m Resolution)

AW3D is a recent global DEM that is derived from the data generated by Panchromatic Remote-sensing Instrument for Stereo Mapping (PRISM), one of the onboard sensors carried in the ALOS. It has a high spatial resolution of 0.15 arc second (approximately 5 m at equator) and target vertical accuracy of 5 m (RMSE). The low-resolution version (30-m spacing) is generated by averaging the original one [START_REF] Takaku | Generation of high resolution global DSM from ALOS PRISM[END_REF] and is provided free of charge. Elevation is estimated by correlation coefficient maximization according to the triplet image matching algorithm with correlation window size optimization [START_REF] Takaku | High resolution DEM generation from ALOS PRISM data[END_REF].

For the PRISM instrument, B = 0.5 (between nadir and forward/aft-looking sensors) and r = 2.5 m. The ratio (r/B) = 5 m/pixel is conversion ratio between pixel and meter units.

Sample AW3D tile for the area of Yushan mountain, Taiwan, (central coordinates are 23 • 28 08.4 N 120 • 57 36.0 E, 4336 × 4702 pixels) were kindly provided by NTT DATA. As for the GDEM2, the AW3D tile was downscaled three times to a 15-m spatial resolution. For the AW3D, experiments do not reveal significance of elevation predictor on measurement error parameters (both σ 2 e and σ 2 Corr , see Table II). This result indicates negligible epipolar line error for AW3D that can be related to higher accuracy of ALOS satellite ephemeris data [START_REF] Tadono | Precise global DEM generation by ALOS PRISM[END_REF] as compared to Terra satellite [START_REF] Fujisada | Advanced methodology for ASTER DEM generation[END_REF]. Similar to GDEM2, measurement error variance is strongly dependent on stacking number. Restricted amount of data available covers only the limited range of N stk change from about 4 to 7. Therefore, σ 2 e for infinite number of N stk cannot be estimated with high precision. Surprisingly, σ 2 e for AW3D and GDEM2 for sea level elevation is of the same order: σ e is about 5 m for N stk = 1. Converted to disparity error in sensor pixels, 5-m SD corresponds to about 0.

D. ALOS Experiment (30-m Resolution)

For the AW3D30 (Table III) as for AW3D, the experiments reveal the significance of stacking number predictors for σ 2 e . For σ 2 Corr , both predictors are not significant. Measurement error variance for the AW3D30 is lower than the one for the GDEM2 even in the best settings (large stacking number and sea level elevation): σ Var01 = 0.6140 m for the AW3D30 as compared to σ Var01 = 1.0145 m for GDEM2. Taking into account that the AW3D30 is obtained by averaging AW3D pixels by six times with respect to both coordinates, σ e for the AW3D30 should be six times lower than for the AW3D. This is actually observed for N stk = 1: σ e.AW3D (1)/6 ≈ 0.83 m as compared to σ e.AW3D30 (1) ≈ 0.78 m. For large values of N stk , measurement error SD for the AW3D30 is comparable to the quantization error. In this extreme case, the mvcNI+fBm might not reveal decay of σ 2 e with stacking number. Measurement error correlation width for the AW3D30 is about 0.2556 pixels at the resolution of 90 m.

E. TanDEM-X-DEM Experiment (90-m Resolution)

In the next experiment, we apply the mvcNI+fBm to the open 90-m resolution version of TanDEM-X-DEM-DEM created by interferometric processing of data from the two twin SAR satellites TerraSAR-X and TanDEM-X [START_REF] Rizzoli | Generation and performance assessment of the global TanDEM-X digital elevation model[END_REF]. The unique feature of TanDEM-X-DEM important for this study is that it contains estimates of the DEM random error SD called height error map (HEM). The HEM is derived by rigorous error propagation on the basis of TanDEM-X platform characteristics, interferometric coherence, and geometrical considerations [START_REF] Rizzoli | Generation and performance assessment of the global TanDEM-X digital elevation model[END_REF], [START_REF] Wessel | TanDEM-X Ground Segment-DEM Products Specification Document[END_REF]. HEM is provided for each DEM tile as auxiliary files with "HEM" suffix and contains pixel-wise estimate of TanDEM-X-DEM random error SD. For TanDEM-X-DEM, random error is understood in the same manner as in this work: "random errors are high-frequency errors with low spatial correlation contributing to both the pointto-point relative vertical accuracy and the absolute vertical accuracy" [START_REF] Wessel | TanDEM-X Ground Segment-DEM Products Specification Document[END_REF]; therefore, we consider HEM as σ e estimate.

Random error for a DEM generated by a radar platform is dependent on many factors. The goal of this experiment is not to analyze these factors, but to verify agreement of HEM and mvcNI+fBm estimates. For this, we consider regression model in the form 1 + HEM 2

) 18
Assuming HEM is a correct estimate of TanDEM-X-DEM random error SD and mvcNI+fBm provides accurate estimates of σ e , one expects σ 2 e = HEM 2 and coefficients of model ( 18) should be σ 2 Var01 = 0 and c VarHEM = 1. For TanDEM-X-DEM, we used the same DEM tiles and the same experiments settings as for the ASTER GDEM2 experiment. The complexity of TanDEM-X-DEM analysis exceeds that of ASTER GDEM2 because random noise SD is two orders of magnitude lower. For the considered tiles, HEM varies from 0.018 (0.01% quantile) to 3.5 m (99.99% quantile) with the mode equal to 0.05 m as compared to 2 . . . 6 m for 18) are significant. For values of HEM exceeding 0.2 m, σe estimated with the mvcNI+fBm is close to HEM and exceeds it by about 1.68 times. For HEM < 0.2, ratio between σe and HEM increases and reaches 2.2 for HEM equal 0.05. This result means that estimates provided by the mvcNI+fBm are close to theoretically predicted ones in the wide range of HEM. Discrepancy between the mvcNI+fBm estimates and HEM might be in part attributed to the fact that HEM is an optimistic estimate of TanDEM-X-DEM error; random noise with SD higher than HEM was reported in [START_REF] Rizzoli | Relative height error analysis of TanDEM-X elevation data[END_REF]. Significance of the HEM component is illustrated by the partial residual plots (Fig. 8).

VI. CONCLUSION

In this article, we have proposed and investigated application of the blind noise parameter estimator to characterize gridded DEM vertical error, specifically fine-scale elevation measurement error.

While BNPE is a well developed area, the new application scenario has not been covered by existing methods: the elevation measurement error SD is dependent on several predictors, while existing methods deal with signal-dependent noise model with only one predictor, namely, image intensity. Therefore, the recently proposed mvcNI+fBm estimator that is able to deal with multivariate noise signal-dependence has been modified and used to estimate both SD and spatial correlation width of DEM measurement error.

In Section V, the mvcNI+fBm has been applied to build bivariate models of ASTER GDEM2 and AW3D (both with 5-and 30-m spatial resolution) elevation measurement error. The two predictors in these models are the number of stereo pairs (responsible for stacking procedure influence) and elevation itself (responsible to epipolar line error influence). These models have been found consistent with the accuracy analysis results published in the available literature for GDEM2 and AW3D data. Not previously reported in the literature, our analysis reveals the epipolar line error as a very important factor responsible for the GDEM2 quality degradation for high elevations. The derived regression models for GDEM2, AW3D30, and AW3D can be used to predict elevation measurement error parameters for low relief areas (with elevation SD from 0 to about 15 m).

Experiment with TanDEM-X-DEM (90-m resolution) shows that estimates of the mvcNI+fBm closely follows theoretically predicted SD (HEM) provided by the German Aerospace Center (DLR). It is important to mention that the mvcNI+fBm was validated in very wide range of measurements error SD, from about 5 m for ASTER GDEM2 down to 0.05 m for TanDEM-X-DEM.

While the proposed no-reference approach to DEM accuracy analysis cannot provide information on DEM bias, the detailed study of DEM measurement error SD and correlation width can see variety of applications including the analysis of LiDAR or RADAR-derived raster DEMs, direct efficiency comparison of different approaches for stereo matching or interferometric phase reconstruction, and the analysis of sensor-related (e.g., jitter causing the epipolar line error of GDEM2) components in elevation measurement error. Future work is seen in the direction of using additional observable predictors (especially related to vegetation cover). Another direction is to overcome the inability of characterizing DEM error in high relief areas.

Fig. 1 .

 1 Fig. 1. Illustration of the DEM error structure.

Fig. 2 .

 2 Fig. 2. Illustration of the local DEM error structure. Heterogeneous areas are marked by gray fill.

  A c c e p t e d M a n u s c r i p t

Fig. 3 .

 3 Fig. 3. Illustration of epipolar line error.

Fig. 4 .

 4 Fig. 4. Processing pipeline of the mvcNI+fBm method.

A c c e p t e d M a n u s c r i p t Dx , Ĥq , σ 2 Corr

 2

); 8 .

 8 Apply robust heteroscedastic regression to ( p j , s j ) pairs to refine DEM context-dependent error parameter vector:θ e.k or θ Corr.k ; 9. Repeat 4-8 for both error variance σ 2 e and spatial correlation width σ 2 Corr ; 10. Increase iteration number: k = k + 1; 11. Repeat the steps from 4 to 10 till convergence defined as ||θ e.k -θ e.k-1 || < ε and ||θ Corr.k -θ Corr.k-1 || < ε, where ε is small constant;

Fig. 5 .

 5 Fig. 5. Illustration of GDEM2 gross errors related to the stacking procedure. (a) GDEM2 patch, (b) GDEM2 stacking number, and (c) AW3D30 patch. Black corresponds to the elevation of 65 m and N stk = 0, and white corresponds to 290 m and N stk = 11.

σ 2 e

 2 = 1.0293 m 2 + 25.6667 m 2 N -1 stk + 4.8991 • 10 -7 Z 2 + 6.1069 • 10 -6 Z 2 N -1 stk . (16)The partial residual plots for components of the model (16) with the highest R 2 (Fig.6) reveal significance of

A c c e

  p t e d M a n u s c r i p t

  64 pixels at the ASTER resolution of 15 m. This value is about the half-size of 5 × 5 pixels correlation window, what is reasonable. For Z = 5000 m, correlation width increases almost twice to σ Corr = 0.8 or 4.8 ASTER pixels. The value of correlation width estimate close to the correlation window size indicates that the measurement error correlation function shape could deviate from Gaussian one at least for high elevation values. Additional research is needed to clarify this hypothesis.

Fig. 6 .

 6 Fig. 6. Partial regression plots for (a), (c), and (e) σ 2 e and (b), (d), and (f) σ 2 Corr model coefficients.

Fig. 7 .

 7 Fig. 7. Comparison of the mvNI+fBm results for the GDEM2 with the results obtained by Becek [22].

  2 pixels for GDEM2 [5 m/25(m/pixel)] and 1 pixel for AW3D [5 m/5(m/pixel)]. Such an elevated error SD may be related to peculiarities of stereo-matching algorithm implementation for AW3D. Measurement error correlation width for the AW3D is about 0.39 pixels at the resolution of 15 m.

  A c c e p t e d M a n u s c r i p t

σ 2 e = σ 2

 2 Var01 + c VarHEM HEM 2 . (

Fig. 8 .

 8 Fig. 8. Partial regression plots for σ e HEM component. ASTER GDEM2. Measurement error estimation results are summarized in TableIV. Both σ 2 Var01 and c VarHEM coefficients of the model (18) are significant. For values of HEM exceeding 0.2 m, σe estimated with the mvcNI+fBm is close to HEM and exceeds it by about 1.68 times. For HEM < 0.2, ratio between σe and HEM increases and reaches 2.2 for HEM equal 0.05. This result means that estimates provided by the mvcNI+fBm are close to theoretically predicted ones in the wide range of HEM. Discrepancy between the mvcNI+fBm estimates and HEM might be in part attributed to the fact that HEM is an optimistic estimate of TanDEM-X-DEM error; random noise with SD higher than HEM was reported in[START_REF] Rizzoli | Relative height error analysis of TanDEM-X elevation data[END_REF]. Significance of the HEM component is illustrated by the partial residual plots (Fig.8).

  s c r i p t

  The residual r (t, s) = [Z grid (t, s) -Z grid.ref (t, s)] + [e grid.mes (t, s, p(t, s))-e grid.ref.mes (t, s, p ref

ref (t, s) == Z grid.ref (t, s) = Z grid.ref = const, e grid.ref.mes (t, s, p ref (t, s)) = 0, and r (t, s) simplifies to [Z grid -Z grid.ref ] + e grid.mes (t, s, p(t, s)). Even not knowing the actual evaluation Z grid.ref , variance of DEM error can be directly estimated as follows:

	Var(e grid.mes

TABLE I RESULTS

 I OF THE ELEVATION MEASUREMENT ERROR PARAMETER ESTIMATION OBTAINED WITH THE mvcNI+fBm FOR THE ASTER GDEM2 each component. Note the practical absence of outliers among mvcNI+fBm estimates (outliers percentage is about 0.1% for σ 2 e

TABLE II RESULTS

 II OF THE ELEVATION MEASUREMENT ERROR PARAMETER ESTIMATION OBTAINED WITH THE mvcNI+fBm FOR THE ALOS 3D WORLD, 5 m TABLE III RESULTS OF ELEVATION MEASUREMENT ERROR PARAMETER ESTIMATION OBTAINED WITH THE mvcNI+fBm FOR THE ALOS 3D WORLD, 30 m

TABLE IV RESULTS

 IV OF THE ELEVATION MEASUREMENT ERROR PARAMETER ESTIMATION OBTAINED WITH THE mvcNI+fBm FOR THE TANDEM-X-DEM, 90 m