Adaptive Multiscale Estimation for Fusing Image Data

1. Adaptive Multiscale Estimation for Fusing Image Data

2. Outline • Introduction • Physical modeling • Data fusion methodology • Results • Conclusions

3. Research Motivation • Develop general capability for mapping and updating topography from multiple sensors • Historical data acquired from different sensors • Acquisitions often have different extents and resolutions • Sensor-dependent height uncertainty • Develop methodology for computing vegetation heights and surface topography in low relief environments • Some applications require ≤10 m horizontal resolution and ≤1 m vertical accuracy • Hydrology: shallow water runoff channels (~0.1 m vertical accuracy) • Seismology: active faults (~1 m vertical accuracy) • Best technologies for low-relief topographic mapping • Interferometric synthetic aperture radar (INSAR) • Covers large area, but insufficient accuracy • Laser altimeter (LIDAR) • Excellent accuracy, but covers small area

4. zg < zS < zv Measuring Topography with INSAR • Problem: no direct measurement of zg in presence of vegetation • INSAR data provide height of phase scattering center zS • Cannot distinguish surface elevation zg from vegetation elevation zv • Neglecting noise, zS = zg for bare surfaces • Proposed solution: • Estimate zg and Dzv from INSAR data using electromagnetic scattering model • Incorporate additional high-resolution measurements (LIDAR) zg = ground height Dzv = vegetation height zS = scattering center height (measured height) Dzv

5. Contributions (1) Combine physical modeling with multiscale Kalman filter to accommodate nonlinear measurement-state relations • Previous multiscale Kalman filters applied to problems in which • Observations are linearly related to state • Closed-form inverse exists for nonlinear measurement-state relationship • No closed-form inverse for multisensor estimation of topography in general • Physical-model inversion yields new set of observations • Enables application of Kalman-based fusion methods to estimating topography (2) Develop spatially-adaptive multiscale Kalman filter • Kalman filter is optimal (in mean squared sense) if filter parameters are correct • Adaptive implementation can update incorrect process noise variance (3) Improve estimates of zg and Dzv for remote sensing applications • True estimate errors are smaller than those obtained by other SAR methods (excluding studies that use controlled training sites with a priori information) (4) Develop framework to update historical data with newer or complementary data • Robust fusion with independent weighting for each data type

6. Outline • Introduction • Physical modeling (contribution #1) • Characterizing INSAR measurements • Obtaining estimates of zg and Dzv • Data fusion methodology • Results • Conclusions

7. INSAR and LIDAR Imaging • INSAR (nominal) • Side-looking • Single or repeat pass • Airborne or space-based • Fixed illumination • C-band - 6 cm wavelength • Vertical accuracy ≥ 2 m • 5-25 m pixel spacing Large coverage area: primary sensor • LIDAR (nominal) • Downward-looking • Airborne (until GLAS 2002) • Scanning illumination • 1 mm wavelength • Vertical accuracy ≥ 0.1 m •  1-5 m pixel spacing complementary sensor

8. INSAR Processing • Antennas 1 and 2 • Transmit modulated chirp pulses • Receive backscattered energy • Quadrature demodulation of received signals yields two complex-valued images and • [n1, n2] are pixel coordinates • Compute normalized interferometric cross correlation (NICC) Single-pass INSAR b = baseline length hS = INSAR altitude a = baseline angle qS = incidence angle zS = height rS1, rS2 = path lengths to yS = horizontal distance antennas 1 and 2 • Compute height of phase scattering center zS[Calculating height][H.O.T.] • Magnitude of NICC function of scattering coherence of target • Phase of NICC function of differential path length to target

9. Addressing Vegetation Effects in INSAR • Empirical statistical relationships • Most common approach to date • Relate NICC to backscattering coefficient so using regression [Wegmuller & Werner, 1995] • Use regression equations to distinguish forest types • Results apply only to training sites used in regression • Relate volume scattering to vegetation height • Calculate tree heights from NICC phase [Hagberg, Ulander, & Askne, 1995] • Assume nearly opaque canopy so phase scattering center at tree tops • Heights 50% underestimated when forest not extremely dense • Relate zg and Dzv directly to INSAR measurements • Use interferometric scattering model M[scattering model] [Treuhaft, Madsen, Moghaddam, & van Zyl, 1996] • No assumptions on vegetation density (t = extinction coefficient) required • Nonlinear optimization (iterative)

10. x = parameter vector d = observation vector M= nonlinear mapping (scattering model) Contribution #1 • Scattering modelMis nonlinear mapping from parameters to observations • Need M-1 to solve for terrain parameters x, but no closed-form inverse exists • Invert M numerically using nonlinear constrained optimization [SQP] • Use two INSAR observations (2 baselines) to solve for three parameters [baseline sensitivity] • Estimate zg and Dzv from INSAR data by solving the objective function[vegetation sensitivity] • Improvements: • Formulate the problem in a constrained optimization framework • Solve the problem at the pixel level (no spatial averaging required) • Fuse estimates with LIDAR to improve accuracy obtained from model inversion

11. Measuring Topography with LIDAR • LIDAR measures zv directly • Optical wavelengths do not penetrate vegetation, except in gaps • Process raw data to obtain zg and Dzv using first and last returns [Weed & Crawford, 2001] • Simplified algorithm • N is the set of all LIDAR pixels • Now have zg and Dzv from both INSAR and LIDAR data

12. Outline • Introduction • Physical modeling • Data fusion methodology (contribution #2) • Kalman filter for data fusion • Adaptive estimation • Multiscale Kalman filter • Results • Conclusions

13. Kalman Filter Model • Use state-space approach • Can model any random process having rational spectral density function with finite state dimension [Brown & Hwang 1997] • Can estimate internal variables not directly observed [block diagram] • Able to track non-stationary data • Use discrete formulation • Data from sampled (imaged) continuous process • Noise sequences w and v have white autocorrelation and zero cross-correlation

14. Kalman Filter Algorithm • Kalman filter is widely used to estimate stochastic signals • Linear, time-varying filter • Implemented in time domain by a recursive algorithm • Requires prior model for filter parameters {F, Q, H, R} • Bounded estimate error covariance Pk|k • Reach steady state if {F, H} are constant and {w, v} are WSS [Grewal and Andrews 1993] input output

15. Determining Model Parameters • R determined by sensor characteristics and data • RS: standard deviation of heights (meters) due to phase noise [Madsen, Martin, & Zebker, 1995] • RL: standard deviation of heights (meters) due to pulse distribution • H is a binary indicator function • INSAR and LIDAR data transformed into estimates of zg and Dzv prior to fusion • H=1 where observations are available • State parameters F and Q determined by signal (data) model • Robust techniques exist for adaptively estimating Q [Mehra, 1972] • Corrections for F require additional knowledge of process dynamics • R, H, and F assumed correct • Assume any model errors reside in Q • implement adaptive correction for Q

16. Detecting Errors in Q • Innovations represent prediction error uk = yk - Hxk|k-1 = Hek + vk • Where ek = xk - xk|k-1 denotes error in estimate • Sequence uk is Gaussian, white sequence for optimal filter[statistics] • Model errors cause assumptions of uncorrelated noise to be violated • Yield correlation in uk , E{uk uTj} ≠ 0 in general • Detect correlation uk using autocorrelation function (ACF) • Non-white ACF(uk) implies model errors • Relate model parameters to ACF(uk) to update Q [estimating Q] • Innovation-correlation method [Mehra, 1970] ^ f(Q)

17. Effect of Model Errors in Kalman Filter • 1-D simulation, dimensionless • Autocorrelation function of innovations not white for non-adaptive filter [setup] • Non-white innovations indicate error in Q Autocorrelation function of innovations with correct Q Autocorrelation function of innovations with incorrect Q

18. Multiscale Data Fusion • Motivation • Captures multiscale character of natural processes or signals • Combines signals or measurements having different resolutions • Common methods • Fine-to-coarse transformations of spatial models • Direct modeling on multiscale data structures, e.g. quadtree [MKS model] [Chou, Willsky, & Benveniste, 1994] • Multiscale signal modeling has been heavily studied in recent years • Multiscale Kalman Smoother (MKS) algorithm • Use fractional Brownian motion data model for self-similar processes like topography[Fieguth, Karl, Willsky, & Wunsch, 1995][stochastic data model]

19. Contribution #2: Adaptive MKS • Q (height variance) not constant for INSAR images • Terrain is non-stationary in general, e.g. forest changing to grassland • MKS assumes Q is uniform at each scale • No mechanism to make Q spatially varying in MKS • 1-D adaptive estimation of Q • Innovation-correlation method used to estimate constant but unknown Q[Mehra, 1970] • Extended to estimate Q locally in sliding window [Noriega and Pasupathy, 1997] • Develop adaptive MKS algorithm (AMKS) [Slatton, Crawford, & Evans, 2001] • Use innovation-correlation method in sliding window of Ns pixels • Assume separable dynamics • Apply 1-D Kalman filters on coarse-resolution image • Incorporate into multiscale framework [algorithm] • Update Q locally by applying spatial Kalman filters to image data in quadtree • Apply to data with dense coverage (INSAR) for physically meaningful innovations

20. filter in scale QAMKS at coarse scale (non-uniform) Low-resolution data (dense) filter in space High-resolution data (sparse) Q at leaf nodes (uniform) Adaptive Multiscale Estimation • Benefits • Filter compensates for modeling errors in Q • Spatially-varying Q better represents images of topography • Updated Q images provide insight about the state process m j i

21. Decreasing support Evolution of Q(m, i, j) • 2-D simulation showing evolution of process noise variance Q(m, i, j) • For QMKS < Qtrue, Adaptive Multiscale Kalman Smoother revises Q upward Root node variance (m2) Qtrue = 32 (uniform) QMKS = 0.3 (uniform) avg(Qtrue – QAMKS) = 25 avg(Qtrue – QMKS) = 35 Leaf nodes

22. Fused Estimates • Simulated observations are ground truth plus measurement noise • Achieve smaller mean squared error (MSE) than MKS Actual error (m2) MSEcoarse obs = 102 MSEMKS = 93.4 MSEAMKS = 79.1 D% (AMKS, MKS) = 15%

23. Computational Complexity • Scattering model inversion • Iterative • Operations per pixel: ~ 15 (nominal), 300 (max) • MKS • Non-Iterative • Operations: • Let N2 = number of leaf nodes, then have floor[4/3 N2] total nodes • Number of operations grows linearly with N2 • Innovation-correlation component of AMKS • Non-Iterative • Operations: • If non-white innovations detected, solve Nlag equations (at one scale only) • Additional complexity versus estimate improvement • Implemented at scale m=M-1, have (N /2)2 nodes • 1.25 times more complex than MKS • Achieves reduction in MSE up to 15% (data dependent) • Current implementation is a development tool • CPU time = 7 min for 256 x 256 on 400 MHz Sun Ultra Sparc, single processor • Computes data model in separate operations and 2-D power spectra • Written in MATLAB with “for loops” (for easier debugging)

24. Outline • Introduction • Physical modeling • Data fusion methodology • Results (contributions #3 and #4) • Application 1: Combine INSAR and LIDAR data adaptively • Colorado River, Austin, Texas • Application 2: Combine Data from multiple INSAR platforms • Finke River, Australia • Conclusions

25. Colorado River Trees Grass fields Texas 2 km Application 1: Combining INSAR and LIDAR • Exploit advantages of both sensors (data fusion) • INSAR is best option for large-scale coverage, LIDAR is best for highly-accurate small-scale mapping (e.g. urban areas) [Wang and Dahman, 2001] • Austin, TX: typical urban floodplain area with trees and grassland • Fuse INSAR @ 10 m with high resolution LIDAR @ 1.25 m • Color infrared aerial photograph • Acquired during winter • Deciduous trees appear gray • Grass and understory appear red

26. no data no data Contribution #3: Improved Estimates Dzv • AMKS algorithm improves estimates of Dzv • Significant improvement over scattering inversion alone • Slight improvement over MKS • Impact of AMKS lessened because implemented 3 levels above LIDAR no data MSEinversion = 9.33 m2 MSEMKS = 5.56 m2 MSEAMKS = 5.55 m2 D% (inversion, MKS) = 40.4 % D% (inversion,AMKS) = 40.5 % no data

27. Improved Estimates zg • AMKS algorithm improves estimates of zg • Sparse LIDAR acquisition MSEraw = 6.1 m2 MSEinversion = 1.3 m2 MSEAMKS = 0.1 m2 D% (raw,AMKS) = 98 % D% (inversion,AMKS) = 92 %

28. ERS MacDonnell Ranges Missionary Plain Finke River James Range Latitude = 132o35’E Longitude = 23o45’S N TOPSAR 50 km Contribution #4: : Updating INSAR Data • Update topographic map with data fusion • Spaceborne ERS (1996) provides large coverage area, but poor vertical accuracy • Airborne TOPSAR (1996 & 2000) provides higher resolution, • Phase unwrapping errors and radar shadowing can produce erroneous data • Finke River study area in Australia • Semi-arid region of geological interest ERS = European Remote Sensing Satellite (European Space Agency) TOPSAR = Topographic SAR (NASA/Jet Propulsion Laboratory)

29. 40 30 20 10 Fuse Three Data Sets • Combine ERS @ 20 m, TOPSAR @ 10 m, and TOPSAR @ 5 m • Phase Unwrapping error in ERS is compensated by fusion

30. Fused Estimates • Differences in resolution apparent [quadtree estimation] • Transects show details for estimator performance • Estimates track ERS when TOPSAR not available, • High RERS prevents filter from tracking noisy ERS data too closely • Estimated process lies between the TOPSAR observations

31. Perspective View [ estimate uncertainty]

32. Outline • Introduction • Physical modeling • Data fusion methodology • Results • Conclusions

33. Conclusions • Combining physical modeling and spatially-adaptive multiscale estimation yields near-optimal estimates of zg and Dzv (MSE sense) • Under assumption of separable dynamics and variations in Q with extent >2Ns pixels • Physical modeling and target-dependent measurement errors allow proper fusion of dissimilar data • Adaptive framework accommodates errors in Q • Significantly smaller mean squared error than MKS in simulations • Modest improvement with data from study site • INSAR 8 times lower resolution than LIDAR • Spatial variations in true Q are small • Algorithm does not track rapid changes (<Ns pixels) in Q • Roughly 15% additional computation time • Contributions • Combined physical modeling with multiscale estimation to accommodate nonlinear measurement-state relations • Developed spatially-adaptive multiscale Kalman filter • Improved estimates of zg and Dzv for remote sensing applications [comparison] • Updated and improved topographic imagery using multiple INSAR data sets

34. Future Work • Non-separable 2-D Kalman filters • Implement spatially-adaptive component of AMKS using a reduced-update Kalman filter [Woods & Radewan, 1977] • Accommodates 2-D spatial correlation in images • Requires extending innovation-correlation method to 2-D • Non-iterative methods for estimating zg and Dzv from INSAR • Current nonlinear optimization is most computationally intensive part of the data fusion process • Investigate methods that relate INSAR decorrelation to volume scattering from vegetation [Rodriguez & Martin, 1992] • These methods are approximate but non-iterative • Multi-model (filter bank) approaches • Develop sets of Kalman model parameters for generalized terrain types, e.g. forest and grassland • Implement a different filter for each terrain, or a linear combination of filters • Combining with image segmentation could avoid Ns latency of innovation-correlation approach • Change detection • Flexible incorporation of new information from different sensors • Weight recent observations more heavily (R)