Unifying Experiment Design and Convex ... - Semantic Scholar

82

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 48, NO. 1, JANUARY 2010

Unifying Experiment Design and Convex Regularization Techniques for Enhanced Imaging With Uncertain Remote Sensing Data—Part I: Theory Yuriy V. Shkvarko, Senior Member, IEEE

Abstract—This paper considers the problem of high-resolution remote sensing (RS) of the environment formalized in the terms of a nonlinear ill-posed inverse problem of estimation of the power spatial spectrum pattern (SSP) of the wavefield scattered from an extended remotely sensed scene via processing the discrete measurements of a finite number of independent realizations of the observed degraded data signals [single realization of the trajectory signal in the case of synthetic aperture radar (SAR)]. We address a new descriptive experiment design regularization (DEDR) approach to treat the SSP reconstruction problem in the uncertain RS environment that unifies the paradigms of maximum likelihood nonparametric spectral estimation, descriptive experiment design, and worst case statistical performance optimization-based regularization. Pursuing such an approach, we establish a family of the DEDR-related SSP estimators that encompass a manifold of algorithms ranging from the traditional matched filter to the modified robust adaptive spatial filtering and minimum variance beamforming methods. The theoretical study is resumed with the development of a fixed-point iterative DEDR technique that incorporates the regularizing projections onto convex solution sets into the SSP reconstruction procedures to enforce the robustness and convergence. For the imaging SAR application, the proposed DEDR approach is aimed at performing, in a single optimized processing, adaptive SAR focusing, speckle reduction and RS scene image enhancement, and accounts for the possible presence of uncertain trajectory deviations. Index Terms—Descriptive experiment design, radar imaging, regularization, spatial spectrum pattern (SSP), synthetic aperture radar (SAR).

L IST OF A CRONYMS ASF DEDR EO ML MSF MVDR POCS RASF RMVDR

Adaptive spatial filter. Descriptive experiment design regularization. Equation of observation. Maximum likelihood. Matched spatial filtering. Minimum variance distortionless response. Projections onto convex sets. Robust ASF. Robust MVDR.

Manuscript received May 15, 2008; revised October 16, 2008 and February 26, 2009. First published September 29, 2009; current version published December 23, 2009. The author is with the Department of Electrical Engineering, Centro de Investigación y de Estudios Avanzados del Instituto Politécnico Nacional, Guadalajara 45015, Mexico (e-mail: [email protected]). Digital Object Identifier 10.1109/TGRS.2009.2027695

RS SAR SFO SSP WCSP

Remote sensing. Synthetic aperture radar. Signal formation operator. Spatial spectrum pattern. Worst case statistical performance. I. I NTRODUCTION

S

PACE–TIME signal processing for imaging with sensor arrays and synthetic aperture radar (SAR) systems has been an active research area in the environmental remote sensing (RS) for several decades, and many sophisticated techniques are now available (see among others [1]–[5] and the references therein). In this paper, we address a new approach to enhanced RS imaging stated and treated as an ill-posed nonlinear inverse problem with model uncertainties. The problem at hand is to perform high-resolution reconstruction of the power spatial spectrum pattern (SSP) of the wavefield scattered from the extended remotely sensed scene via space–time processing of finite recordings of the RS data distorted in a stochastic uncertain measurement channel. The SSP is defined as a spatial distribution of the power (i.e., the second-order statistics) of the random wavefield backscattered from the remotely sensed scene observed through the integral transform operator [4]– [6]. Such an operator is explicitly specified by the employed radar signal modulation and is traditionally referred to as the signal formation operator (SFO) [4]–[6]. The classical imaging with an array radar or SAR implies application of the method called “matched spatial filtering” (MSF) to process the recorded data signals [7]–[9]. Stated formally, the MSF method implies application of the adjoint SFO to the recorded data, squared detection of the filter outputs, and their averaging over the actually recorded samples (the so-called snapshots [4]) of the independent data observations. A number of approaches had been proposed to design the constrained regularization techniques for improving the resolution in the SSP obtained by ways different from the MSF, e.g., [9]–[17] but without aggregating the descriptive regularization principles with the minimum risk statistical estimation strategy. Although the existing theory offers a manifold of statistical and deterministic regularization techniques for reconstructive imaging [3], [9]–[15], adaptive image despeckling [3], [18], and problem-oriented image enhancement [12], [15], [19]–[23] in many application

0196-2892/$26.00 © 2009 IEEE

SHKVARKO: UNIFYING EXPERIMENT DESIGN AND CONVEX REGULARIZATION TECHNIQUES—PART I

areas, there still remain some unresolved crucial theoretical and processing problems related to large-scale radar/SAR reconstructive imaging in uncertain operational scenarios. The predominant challenge of this paper is to solve the inverse problem of high-resolution SSP estimation in the uncertain RS environment. The operational scenario uncertainties are associated with the unknown statistics of multiplicative signal-dependent noise and random SFO perturbations in the turbulent medium, imperfect array calibration, finite dimensionality of measurements, uncontrolled antenna vibrations, and random carrier trajectory deviations in the case of SAR. We address a new descriptive experiment design regularization (DEDR) approach to treat such SSP reconstruction problems that unifies the paradigms of maximum likelihood (ML) nonparametric spectral estimation, descriptive experiment design, and worst case statistical performance (WCSP) optimizationbased regularization. Pursuing such an approach, we establish a family of the DEDR-related SSP estimators that encompass a manifold of algorithms ranging from the traditional matched spatial processing to the new robust adaptive spatial filtering and minimum variance distortionless response (MVDR) beamforming methods. On one hand, we found that both the robust MVDR (RMVDR) and the developed here DEDR can be unified under some specifications of the controllable degrees of freedom of the DEDR method. On the other hand, the DEDRrelated estimators allow computational implementation scheme that does not involve the direct inversion of the estimated data correlation matrix. The latter property makes the DEDRrelated techniques applicable to the scenarios with the lowrank uncertain estimated correlation matrices, e.g., rank-1 in the case of a single-look SAR, in which the only single random realization of the trajectory data signal is processed to form the scene image [2], [11], [18]. In this paper, we present the mathematical foundation of the DEDR method. The theoretical study is resumed with the development of a fixed-point iterative DEDR technique that incorporates the regularizing projections onto convex sets (POCS) into the SSP reconstruction procedures to enforce the robustness and convergence. More practical aspects, performance issues and extensive simulation results are treated elsewhere [1]. This companion paper presents the robust numerical versions of the unified DEDR method particularly adapted to the imaging SAR application, in which the DEDR approach is aimed at performing, in a single optimized processing, adaptive SAR focusing, speckle reduction and RS scene image enhancement, and accounts for the possible presence of uncertain trajectory deviations. The rest of this paper is organized in the following manner. The theoretical background and problem formulation are outlined in Section II. In Section III, the ML and minimum risk statistical motivations for the DEDR method are presented. The unified DEDR method is developed in Section IV with the detailed specification of the family of the DEDRrelated techniques presented in Section V. Section VI addresses the regularization approach to iterative implementation of the DEDR techniques performed via incorporating into the solution procedure the projections onto convex solution sets ruled by the employed fixed-point iterative process. Section VII summarizes this paper.

83

II. B ACKGROUND A. Problem Statement The generalized mathematical formulation of the problem at hand presented here is similar in notation and structure to that in [5], [6], [19], and [21], and some crucial elements are repeated for convenience to the reader. Consider a coherent RS experiment in random medium and the narrowband assumption [2], [20], [21] that enables us to model the extended object backscattered field by imposing its time invariant complex scattering (backscattering) function e(r) in the scene domain (scattering surface) R r. The measurement data wavefield u(p) = s(p) + n(p) consists of the echo signals s and additive noise n and is available for observations and recordings within the prescribed time–space observation domain P = T × P, where p = (t, ρ) defines the time(t)–space(ρ) points in P ; t ∈ T , ρ ∈ P; p ∈ P . The model of the continuous-form observation field u in the conventional baseband representation format [2], [4], [11] is defined by specifying the stochastic equation ˜ + n; of observation (EO) of an operator form [23]: u = Se ˜ e = e(r) ∈ E(R), u = u(p) ∈ U(P ), n = n(p) ∈ U(P ); S: E(R) → U(P ), in the function signal spaces E and U with the induced by the scalar (inner) standard L2 metrics structures , u ] = u (p)u∗2 (p)dp, and [e1 , e2 ]E = products [8], [u 1 2 U 1 P ∗ R e1 (r)e2 (r)dr, respectively, where the asterisk indicates the complex conjugate. In the conventional integral form [21], [23], such stochastic EO is given by ˜ u(p) = (( Se(r))(p) + n(p) ˜ r)e(r)dr + n(p) S(p,

=

(1)

R

where the complex zero-mean scattering function e(r) represents the random scene reflectivity over the scene frame R; r ∈ R is a coordinate vector of the scan parameters over the illuminated scene, usually the polar, cylindrical, or Cartesian coordinates. In general, the scattering function e(t, r) may vary slowly in time but these variations are negligible during the time interval in which the particular point-type scattering element dr ∈ R (generic point target) remains within the illuminated patch (standard spatial narrowband assumption [2], [4], [11]). That is why, following the radar imaging applications, we adopt the time invariant model of the random scattering function e(r) at a particular observed data realization (1), while for different independent observed realizations of the data {u(j) (p)}, the scattering functions {e(j) (r)} randomly vary from realization to realization {j = 1, . . . , J}. In (1), the functional kernel ˜ defines the signal ˜ r) of the stochastic integral SFO S S(p, ˜ r) is wavefield formation model. Its mean S(p, r) = S(p, referred to as the nominal SFO in the RS measurement channel specified by the time–space modulation of signals employed in a particular imaging radar/SAR system [2], [4], [11], and ˜ r) − S(p, r) the variation about the mean δS(p, r) = S(p, pertains to the random distortion component in the SFO; it models stochastic perturbations of the wavefield at different propagation paths (ρ, r). In the Section II-C, we will detail the propagation model-based characterization of such the disturbed SFO. In addition, in all RS imaging applications, the antenna

84


elements are considered to be sparse aperture antennas (identical or different but never omnidirectional sensors). Thus, the theory that we develop in this paper is applicable for an arbitrary model of the RS system with one transmit aperture antenna and L ≥ 1 receive aperture antenna elements, e.g., a conventional sparse aperture array imaging radar or a SAR (in the latter case, the array is constructed by one moving antenna). In the simulations reported in the companion paper [1], we detail the adopted practical model of the SFO specified by the particular employed signal modulation and SAR configuration/ geometry. It is conventional in the RS applications to assume that due to the integral signal formation model (1), the central limit theorem conditions hold [3], [8], [24], hence the fields e, n, u in (1) are considered to be the zero-mean complex-valued random Gaussian fields. Note that pursuing here the general regularization approach to constructing the robust space–time processing method for RS imaging, we assume the resulting observation field u to be inhomogeneous [24], [27]. Next, since in all RS applications the regions of high correlations of e(r) are always small in comparison with the element of resolution on the probing scene [2], [18], [24], the scattered signals e(r) from different directions r, r ∈ R are assumed to be uncorrelated, i.e., characterized by the correlation function Re (r, r ) = e(r)e∗ (r ) = b(r)δ(r − r )) ,

r, r ∈ R

(2)

where δ(·) defines the delta-function and · is the averaging operator. The ensemble average r∈R (3) b(r) = e(r)e∗ (r) = |e(r)|2 ) , of the squared modulus of the random scattering field e(r) as a function over the analysis domain (scene frame) R r has a statistical meaning of the average power scattering function and is traditionally referred to (in the RS and radar imaging literature, e.g., [4], [5], [6], and [16]) as the SSP of the scattered field. Representing the spatial distribution of the average power of the random scatterers, the SSP b(r) characterizes in an explicit statistical sense the brightness reflectivity of the scene being mapped (for this reason, b is adopted in the notations as an abbreviation from brightness reflectivity). Note that other terms such as spatial power spectrum, power reflectivity, or differential radar cross section are sometimes used instead of SSP [2], [9]–[11], [13], [16], [24]–[26]. In this paper, we generalize all these definitions by the term SSP. The radar/SAR imaging problem is to derive an estimate ˆb(r) of the SSP distribution (3) over the scene R r by processing whatever available measurements of the actual data recordings (1) collected with a particular system operating in the uncertain RS environment. B. Experiment Design Formalism for Data Representation The formulation of the data discretization and sampling in this paper follows the unified formalism given in [5]–[7] and [19] that enables one to generalize the finite-dimensional approximations of (1), (3) independent of the particular sys-

tem configuration and the method of data measurements and recordings employed. Following [5] and [19], consider an array composed of L antenna elements characterized by a set of complex amplitude-phase tapering functions {κ∗l (ρ); ρ ∈ P; l = 1, . . . , L} (the complex conjugate is taken for convenience). In practice, the array (synthesized array) is composed of identical antenna elements distanced in space, i.e., the tapering functions {κ∗l (ρ)} have the disjoint supports in P ρ (do not overlap). Hence, {κl (ρ)} compose an orthogonal set: [κl , κl ]U = κl 2 δll ∀l, l = 1, . . . , L, where δll is the Kronecker operator and κl 2 = [κl , κl ]U represents the squared norm of the corresponding function. For identical array elements, all κl 2 take the same value κl 2 = κ 2 , thus the array can be easily calibrated by scaling κ 2 = 1 that results in the orthonormal regular spatial decomposition set {[κl , κl ]U = δll }. Consider, next, that the sensor output signal in every spatial measurement channel is then converted to I samples at the outputs of the identical temporal sampling filters defined by their impulse response functions {υi∗ (t); i = 1, . . . , I}. Without loss of generality [4], [6], [19] the set {υi } is also assumed to be orthonormal (e.g., via proper filter design and calibration [7], [19]), i.e., [υi , υi ]U = δii ∀i, i = 1, . . . , I. The composition {hm (p) = κl (ρ)υi (t); m = (l, i) = 1, . . . , M = L × I} ordered by multi-index m = (l, i) composes a set of orthonormal spatial-temporal decomposition (base) functions that explicitly determine the vector of outcomes m = 1, . . . ,M} u = vec{um = [u, hm ]U = u(p)h∗m (p)dp; m

P

(4) of such an M -dimensional (M -D in our notation) data recording channel, in which the employed base functions {hm (p)} span the relevant M -D data representation subspace U(M ) = P U(M ) U = Span(M ) {hm (p)} specifying the corresponding projection operator P U(M ) . In analogy to (4), one can define now the K-D vector-form approximation of the scene random scattering function e = vec{ek = [e, gk ]E = e(r)gk∗ (r)dr; k = 1, . . . , K} k

R

(5) with the elements composed of the decomposition coefficients {ek } with respect to some chosen normalized orthogonal set of expansion functions {gk (r)} that span such K-D signal approximation subspace E(K) = P E(K) E = Span(K) {gk (r)} specifying the corresponding projection operator P E(K) (see [3], [5], [22], [23] for mathematical and signal processing details related to the construction of such projectors P U(M ) and P E(K) ). The descriptive experiment design aspects of the SSP reconstruction problem involving the analysis of how to choose the base functions {gk (r)} that span the signal representation subspace E(K) = Span(K) {gk (r)} for a given observation subspace U(M ) = Span(M ) {hm (p)} were investigated in more details in the previous studies [19], [22], [23]. In this paper, we consider the conventional representation basis formed by a Kx × Ky regular pixel-formatted lattice with


unitary amplitudes and the spacing between lattice points normalized to one pixel width (i.e., ordinary rectangular pixels [3], [19]), where Kx defines the dimension of the rectangular pixel grid over the horizontal (azimuth) coordinate x and Ky defines its dimension over the orthogonal (range) coordinate y. Such regular lattice of pixels {pixk (r)} specified by the ordered multi-index k = (kx , ky ); kx = 1, . . . , Kx ; ky = 1, . . . , Ky ; k = 1, . . . , K = Kx × Ky is practically motivated in all RS imaging applications [2], [3], [19], [26] because the ordinary pixels coincide with their squares {pixk (r) = pix2k (r)}, which makes the same pixel grid {gk (r) = pixk (r)} applicable for discrete-form representation of both the complex scattering function e(r) and the SSP b(r). Using the decompositions (4), (5), we now proceed from the integral continuous-form EO (1) to its discrete (vector-form) approximation +n u = Se

(6)

where the u, n, and e define the vectors composed of the coefficients um , nm , and ek of the finite-dimensional approx imations (4), (5) of the fields u, n, and e, respectively, and S represents the matrix-form approximation of the uncertain SFO with elements [23] ˜ k (r) h∗m (p)drdp; ˜ k , hm ]U = Sg {Smk = [Sg R×P

k = 1, . . . , K;

m = 1, . . . , M }.

(7)

C. Propagation Uncertainties In this section, we present the basic model considerations related to random fluctuation of the fields propagating through the turbulent medium. Note that our treatment is problem oriented, i.e., the statistical model of the uncertain random medium effects is formalized following classical electromagnetic propagation theory [24], [25], [27] in the format required further on for solving the corresponding inverse SSP reconstruction problems in the context of the uncertain operational scenario. First, from the propagation theory [24], [25], [27], we invoke the general-form representation of the noise-free complex distorted field (1), u(p) = s(p), i.e., considering first, for simplicity, no additive noise in the EO (1), i.e., n(p) = 0. The field u(p) is observed in the prescribed time–space domain P = T × P (whereas previously, p = (t, ρ) defines the time(t)–space(ρ) points in P ) and is represented in the “complex phase fluctuation” terms as follows [24]: u(p) = u0 (p) exp {ψ(p)}

(8)

where u0 (p) is a reference field associated with the field propagated through “free space” when the random medium is removed, while all medium perturbations effects are modeled by the complex multiplicative term exp{ψ(p)} with the complex phase given by ψ(p) = χ(p) + iϕ(p).

(9)

85

In terms of the real-valued amplitudes A and A0 and the phases Φ and Φ0 , one can rewrite (8) as u(p) = A(p) exp {iΦ(p)}

u0 (p) = A0 (p) exp {iΦ0 (p)} (10)

χ(p) = ln(A/A0 )

ϕ(p) = Φ(p)−Φ0 (p)

(11)

respectively. Thus, the real part χ(p) of ψ represents the fluctuation of the logarithm of the amplitude A (log-amplitude fluctuation), and the imaginary part ϕ(p) of ψ represents the fluctuation of phase. The multilayer electromagnetic explanation of the model (10) and (11) is as follows [24], [27]. The wave experiences amplitude and phase fluctuations as it propagates through a thin layer. One can write this field as u0 u1 , where u0 is the reference (distortion free) wave and u1 represents the fluctuations. This field u0 u1 is incident on the second layer to produce u0 u1 u2 . In this way, the total output field becomes the product u0 u1 u2 u3 u4 · · ·. Expressing each fluctuation un by an exponent exp(ψn ), the total field is given by u = u0 exp(ψ1 + ψ2 + ψ3 + · · ·), i.e., the exponent ψ in (8) represents the sum of all fluctuations ψn along the propagation path. It is clear that in the imaging radar and/or SAR observation models, ψ = ψdown + ψup , i.e., accumulates the down-path ψdown and up-path ψup complex phase fluctuations [27]. In the commonly used Rytov model [24], the resulting field is approximated by u(p) = u(p) exp [uf (p)/u(p)] ;

u(p) = u0 (p) (12)

where uf (p) represents the so-called first-order multiple scattering solution to the corresponding wave equation for a lineof-sight propagation model [24], and the statistical averaging is performed over the randomness of ψ. Next, assuming that the magnitude of uf /u0 is small compared with unity and expanding, the exponent in a series yields the Rytov approximation u(p) = u0 (p) [1 + uf (p)/u0 (p) + · · ·] = u0 (p) [1 + μ(ρ)] (13) where the time invariant model μ(ρ) = uf (p)/u0 (p) + · · · of multiplicative distortions is usually adopted for a particular small observation time interval. Note that (13) becomes identical to the simplest so-called first-order scattering solution-based approximation [24], if only the first two terms u0 + uf in the series are kept. The perturbation term μ(ρ) in the Rytov model (13) is considered as a zero-mean random field [24], [27] characterized by its correlation function Rμ (ρ, ρ ) = μ(ρ)μ∗ (ρ ). In the first-order approximation (i.e., μ ≈ uf /u0 = χ + iϕ), Rμ (ρ, ρ ) can be next expressed via the corresponding correlation functions: Rχ (ρ, ρ ) = χ(ρ)χ(ρ of the log-amplitude fluctuations, Rϕ (ρ, ρ ) = ϕ(ρ)ϕ(ρ ) of phase fluctuations, and the cross correlations Rχϕ (ρ, ρ )) = χ(ρ)ϕ(ρ). The models of such correlation functions are well studied in the propagation theory, e.g., [24]–[27], where those are expressed over the electromagnetic characteristics of the medium that (for notational simplicity) we specify here by the general term m referred to as “medium model characteristics.” Dependent

86


Fig. 1. Linear array geometry with sensor possible mismatches. ×—nominal sensor positions. —actual sensor positions. m = 1, . . . , M .

on the particular medium propagation model, such characteristics m include [24]–[27]: the particular waveband, polarization, and antenna radiation patterns; the scattering particle size; the angular scattering characteristic of the particle; the medium angular scattering filter function; etc. It is clear that the resulting correlation function Rμ (ρ, ρ ) = Rμ (ρ, ρ : m) is implicitly model dependent, i.e., depends on the particular medium model characteristics m. In the RS experiments, because of random variations of the zero-mean multiplicative amplitude-phase distortions μ = μ(ρ, r) over different propagation paths (ρ, r), the corresponding integral-form representation of the additively noised and multiplicatively distorted data field can be presented by the following expansion of the stochastic EO (1): ˜ (r)) (p) + n(p) u(p) = (Se ˜ r)e(r)dr + n(p) = S(p,

R

=

S(p, r)e(r)dr + R

δS(p, r)e(r)dr + n(p) (14) R

where ˜ r) − S(p, r) = μ(ρ, r)S(p, r) δS(p, r) = S(p,

(15)

represents now stochastic perturbations of the SFO caused by the multiplicative amplitude-phase propagation distortions μ(ρ, r). For the purpose of our study of the inverse radar imaging problems, we assume that the propagation model uncertainties are associated with the a priori unknown medium perturbation statistics, in the sense that the observer has no prior knowledge of the particular values that take the medium model characteristics m in the particular operational scenario, thus cannot specify the particular model Rμ (p, p : m). In other words, the perturbation statistics are conditioned by the set of the medium model characteristics m that in a particular operational scenario are unknown to the observer. D. Composite Observation Uncertainties Fig. 1 shows an example of a linear array with sensor position mismatches {˜ ρm = ρm + δρm }. In the descriptive experiment design formalism, possible displacements of array elements with respect to the presumed nominal positions, as well as distorted antennas shapes are attributed to the unknown disturbances {δhm ; m = 1, . . . , M } in the decomposition functions ˜ m = hm + δhm } [19], [20], [36]. In imaging SAR applica{h tions, such disturbances account to the deviations of a carrier from the nominal trajectory [28]–[30], as shown in Fig. 2.

Fig. 2. SAR system geometry: nominal and in the presence of trajectory deviation.

These position mismatches and random SFO perturbations due to the propagation distortions result in the uncertain SFO ˜ in the vector-form EO (6). The elements of S ˜ are dematrix S fined now by the scalar products (7) with respect to the distorted ˜ m = hm + δhm }. Hence, the uncertain decomposition set {h SFO matrix admits the following representation [19], [23]: ˜ = S + Δ. S

(16)

The nominal M × K SFO matrix S in (16) is composed of the scalar products {Smk = [Sgk , hm ]U }M,K m,k=1 while all problem model uncertainties are attributed to the distortion term Δ which incorporates also the multiplicative signaldependent speckle noise that arises due to the presence of many independent scatterers within the resolution cell, whose scattered field contributions randomly interfere at the receiver [18], [31]. In practice, one cannot attribute the exact portion of the composite SFO perturbation term Δ to a particular source of disturbances, thus cannot separate in (16) the uncertainties caused by the turbulent medium propagation effects, speckle noise, or the observation mismatch errors as those are randomly mixed in the resulting Δ. These practical aspects motivate our adopting the robust statistical treatment of the irregular SFO perturbations Δ as a random zero-mean matrix with the bounded second-order moment [23], i.e.,

(17) Δ = 0;

Δ 2 = tr{ΔΔ+ } ≤ η where Δ 2 = tr{ΔΔ+ } defines the squared matrix norm, tr{·} is the trace operator, superscript + defines the Hermitian conjugate (conjugate transpose), and η is the bounding constant. Thus, the uncertainties regarding the propagation medium distortions and observation model mismatch are combined and formalized by introducing one robust statistical bounding constraint (17) that distinguishes our approach from other studies where these sources of uncertainties were treated in a deterministic fashion, e.g., [20], [21], [32], and [38]. Our proposition for the robust statistical treatment of the model uncertainties (17) is motivated by the practical RS considerations [18], [30], [33]. For example, in the SAR applications [2], [28], [30], the speckle noise, antenna vibrations, and carrier trajectory deviations are always modeled in a statistical way to avoid rather restrictive assumptions of considering and bounding the maximal possible values of some deterministic metric of multiplicative noise, random antenna vibration, and carrier trajectory deviation.


E. Discrete-Form Problem Model Making the use of the perturbed SFO model (16), we can express the degraded vector-form EO (6) as follows: ˜ + n = Se + Δe + n. u = Se

(18)

The correlation matrix Ru of the data vector u is specified via performing the statistical ensemble averaging + ˜ ˜+ p(e) S + nn+ (19) Ru = uu+ = See

over the pixel-formatted observation scene R r by processing the recorded data vector u (in the operational scenario with the single processed data realization) or J > 1 whatever available independent realizations {u(j) ; j = 1, . . . , J} of the data (in the scenario with multiple observations) collected with a particular system operating in the uncertain RS environment. Recall that in this paper, we intend to develop and follow the DEDR strategy. III. M OTIVATIONS FOR DEDR

p(Δ)

where the subscripts p(e) and p(Δ) indicate that the mathematical expectation operations are performed with respect to the independent randomness [factorized probability density functions (pdfs)] of e and Δ. The second term at the right-hand side of (19) defines the observation noise correlation matrix Rn = nn+ . In addition, by definition (3), the correlation matrix of vector e is Re = ee+ = D(b) = diag(b), i.e., a diagonal matrix with the SSP vector b at its principal diagonal composed of the elements bk = ek e∗k = |ek |2 ; k = 1, . . . , K, hence, vector b is referred to as a K-D vector-form representation of the SSP. Using these definitions for Re and Rn , we next express the model data correlation matrix (19) as follows: ˜ ˜+ Ru = SD(b) S

p(Δ)

+ Rn

= SD(b)S+ + ΔD(b)Δ+ p(Δ) + Rn

(20)

where ·p(Δ) defines the averaging performed over the randomness of Δ. In the uncertain observation scenario, the pdf p(Δ) is unknown to the observer. In the simplest case of a given (adopted by the observer) coarse approximation of the pdf model p(Δ) = p(Δ|m) conditioned by some accepted values of the model characteristics m, one can compute the corresponding model-dependent coarse approximation to ΔD(b)Δ+ p(Δ|m) , in which case the approximated model of the data correlation matrix becomes +

Ru = SD(b)S + Rcomp

(21)

where Rcomp = Rcomp (b; m) = ΔD(b)Δ+ p(Δ|m) + Rn represents the term that corresponds to the composite distortions in the observation data. Note that the latter depends on both the SSP vector b and the particular values of the model characteristics m. In the sequence, we consider the general case of uncertain scenarios, in which the data vector u is characterized by the correlation matrix (20) dependent on the model characteristics unknown to the observer. The hypothetical coarse model (21) will be also treated to explain the motivations for the DEDR method. The uncertain inverse problem of SSP reconstruction from the finite-dimensional data measurements is formulated now as ˆ and use it follows: to derive an estimator for the SSP vector b to reconstruct the SSP distribution ˆb(K) (r) =

K k=1

ˆbk gk (r)

(22)

87

A. ML Estimators for Certain Scenarios To motivate the DEDR strategy, in this section, we recall the ML method for statistically optimal estimation of the SSP vector b in the operational scenarios without model uncertain˜ = S, and Rn = N0 I. We ties, i.e., putting Δ = 0, hence, S commence with the simplest scenario, in which the only single recorded realization u = Se + n of the random complex data vector is available for further processing. For the adopted zeromean Gaussian u with the correlation matrix Ru = Ru (b) = SD(b)S+ + Rn , the log-likelihood function of vector b can be expressed in the form of a logarithm of the conditional pdf p(u|b)

ln p(u|b) = − ln det SD(b)S+ + Rn −1 − [u, SD(b)S+ + Rn u]

(23)

where the terms that do not contain b are ignored. The ML solution ˆ = arg max {ln p(u|b)} = arg min {− ln p(u|b)} b b

b

(24)

of the SSP estimation problem implies now the minimization of the objective function

L(b) = − ln p(u|b) = ln det SD(b)S+ + Rn −1 + [u, SD(b)S+ + Rn u].

(25)

Because of a nonlinear dependence of ln det{SD(b)S+ + + −1 in (25) on the desired Rn } and R−1 u = (SD(b)S + Rn ) SSP vector b, no unique regular method exists for deriving the solution to the problem (24) in a closed algorithmic form [3], [34], [35]. Nevertheless, the solution to this problem can be found in the form of a technically tractable adaptive algorithm applying the standard gradient method to minimize (25) with respect to b [5] that yields the ML estimator

˜ ML = FW uu+ F+ b W diag

(26)

where {·}diag defines a vector composed of the elements of the principal diagonal of the embraced matrix and −1 −1 + −1 FW = S+ R−1 S Rn n S+D

(27)

ˆ ML )) Wiener-type is the solution-dependent (i.e., D = diag(b solution operator. The derivation of this ML estimator is detailed in Appendix A, where (A8) corresponds to (26). Observe

88


that the outer product uu+ in the ML estimator (26) is recognized to be a rank-1 (i.e., ill-conditioned) square matrix Y = Y(rank−1) = uu+

(28) C. Partially Uncertain Scenarios

formed from single observed random realization of the data vector u. Technically, the ML optimal algorithm (26) implies, first, application of the Wiener-type adaptive spatial filter (ASF) (27) to the observed data realization that provides an optimal ˆ = FW u of the complex scattering vector e, linear estimate e and second, performing its squared element-by-element detection to form the ML estimates {ˆbMLk = |{FW u}kk |2 ; k = ˆ ML . 1, . . . , K} that compose the desired SSP vector estimate b The extension of this ML estimator for the case of arbitrary recorded observations is obviously performed via averaging the correlations over J actually registered independent data realizations {u(j) ; j = 1, . . . , J} that yields

ˆ ML = {R ˜ e }diag = FW YF+ b (29) W diag where Y represents now the empirical estimate of the data correlation matrix [10], [20] Y = Y(rank−J) = (1/J)

J

u(j) u+ (j)

(30)

j=1

which is the full-rank matrix for J ≥ M . The derivation of this ML estimator is detailed in Appendix A, where (A9) corresponds to (29). Hence, in both operational scenarios, the adaptive (solution dependent) Wiener-type spatial filter (27) is to be constructed and applied to the data to form the ML estimates of the desired SSP vector, in particular, employing the ML algorithm (26) in the scenario with the single recorded data realization and the (29) in the scenario with the multiple independent data recordings, respectively. B. Balance of Risk Components It is well known, e.g., [8], [10], and [19], that for the Wienertype solution operator (27), the minimal value of the Bayesian ˆ = FW u is risk r(F) = Fu − e 2 of the linear estimate e attained, i.e., FW = arg min{r(F)}. Such Bayesian risk funcF

tion r(F) admits the following representation [19]: r(F) = Fu

− e = tr (FS − I)D(FS − I)+ + tr FRn F+

optimal balance between the attained spatial resolution and ˆ = FW u. noise suppression in the linear estimate e

Consider now some observer-adopted approximation Rcomp to the second (distortion) term in the data correlation matrix (21). Two distinct hypothetical scenarios could be treated in such a context. In the first one, the observer completely neglects any electromagnetic model regarding the propagation medium characteristics, as well as any statistical description of the observation uncertainties, and roughly approximates the distortion term in the correlation matrix (21) by the maximum entropymotivated [5], [35] white noise model Rcomp = Ncomp I with the accepted composite noise variance Ncomp . In the second scenario, the observer adopts some coarse electromagnetic model for the solution-dependent component Rcomp = Rcomp (b; m) and tries to expand the ML method for that hypothetical partially uncertain scenario. The first scenario (i.e., Rcomp = Ncomp I) relates directly to the certain one considered in the previous Section III-A. Really, the correlation term Rcomp , in this case, has the same statistical meaning as the noise term Rn in the certain model of the data correlation matrix. Thus, the ML-optimal solutions presented above in the Section III-A. can be easily expanded for that case simply by replacing Rn = Rcomp in all the relevant formulas. Hence, the fixed model case obviously coincides with the certain statistical model scenarios investigated in the previous studies, e.g., [5], [6], and [19]. In the second hypothetical scenario, the observer adopts some model-dependent approximation Rcomp = Rcomp (b; m) to the uncertain correlation term in the model (21) of the data correlation matrix. Let us demonstrate that an attempt to expand directly the ML method for this case leads to the extremely intensive computational problem. The original ML objective function (25) for such partially uncertain correlation data model (21) must be replaced by L(b; m) = − ln p(u|b; m)

= ln det SD(b)S+ +Rcomp (b; m) −1 u . + u, SD(b)S+ +Rcomp (b; m)

(32)

2

(31)

thus, it is decomposable into the systematic error component rs (F) = tr{(FS − I)D(FS − I)+ } and the fluctuation or noise error component rn (F) = tr{FRn F+ }, respectively. The systematic component of risk rs (F) is recognized to be a D-weighted squared matrix norm of the discrepancy between the composition FS and the identity matrix I, i.e., it measures “how close” F is to the pseudoinverse of S; thus, the systematic error term provides a measure of the spatial resolution gained with the employed optimal filter FW . The fluctuation component of risk rn (F) represents the squared norm of noise in the estimate (i.e., the energy of noise in the resulting solution). Hence, the linear minimum risk estimator minimizes the sum of these two error measures rs (F) and rn (F) providing an

Because of a nonlinearity, no unique regular method exists for deriving the solution to the minimization problem with the objective function specified by (32) in a closed algorithmic form. For the certain scenarios (treated in the previous section), the elegant ML solution was constructed in the adaptive algorithmic form employing the independence of the noise correlation term Rn in the data correlation matrix Ru = SD(b)S+ + Rn on the desired solution vector b. In the expanded objective function (32), the composite correlation term Rcomp (b; m) depends on both the SSP b and the uncertainty model characteristics m. Thus, the expanded ML problem of minimization of the objective function (32) results in a computationally extremely intensive (NP-hard) problem [34]. This provides additional motivations to seek for the statistically robust strategies for solving the SSP reconstruction


inverse problems both in the uncertain and partially uncertain operational scenarios. IV. U NIFIED DEDR M ETHOD A. DEDR Strategy We propose the DEDR strategy as an extension of the ML method for the uncertain operational scenarios pursuing the same idea of optimal weighted balancing between the attained spatial resolution and noise suppression in the estimate of the desired SSP vector

ˆ = FYF+ (33) b diag specified now by the DEDR solution operator F = FDEDR . The DEDR estimator (33) employs the outer product Y formed by (28) in the case of the single processed data realization (e.g., a single-look SAR) and the empirical data correlation matrix (30) in the scenario with arbitrary recorded independent data realizations. For both scenarios, the DEDR strategy is aimed at the minimization of the extended risk function + ˜ ˜ p(Δ) }+tr{FRn F+ } with respect to tr{(FS−I)D(F S−I) the solution operator F where the averaging in the first term (systematic risk component) is performed now over the randomness of the uncertain SFO. Observe that such the extended ˜ − I)D(FS ˜ − I)+ p(Δ) } is systematic risk rs (F) = tr{(FS an averaged D-weighted squared norm of the discrepancy ˜ and the identity matrix I. It between the composition FS measures “how far” is the optimal solution operator F from the ˜ in the averaged operator pseudoinverse to the uncertain SFO S metrics induced by the solution-dependent weight matrix D = D(b). In general, the systematic risk can be redefined for any sensible choice of the metrics inducing diagonal weight matrix A composed of positive {Akk > 0; k = 1, . . . , K} providing additional “degrees of freedom” in specifying such the ˜ ˜ S− extended systematic risk rs (F) = tr{(FS−I)A(F + I) p(Δ) }. Putting the idea of optimization of the balance between the attained spatial resolution and noise suppression in the desired solution as a base for the DEDR strategy, we next propose to treat the statistical model uncertainties via employing the constrained regularization approach that incorporates the WCSP optimization into the objective risk minimization problem [20], [21]. Thus, the DEDR strategy for determining the solution operator is formalized by the following constrained optimization problem: {R(F)} F = arg min max 2 F

Δ ≤η

(34)

where ˜ − I)A(FS ˜ − I)+ p(Δ) } + αtr{FRn F+ } R(F) = tr{(FS (35) represents the descriptive objective risk function, in which the balance (regularization) parameter α and the weight matrix A constitute the controllable/adjustable “degrees of freedom”: α is viewed as a tolerance factor that balances the systematic and noise measures in the risk function (35), while A induces

89

the weighted metrics structure in the systematic risk function defined by the first term at the right-hand side of (35). The DEDR strategy (34), (35) is recognized to be a structural extension of the statistical ML method outlined in the previous section. It implies the minimization of the weighted sum of the systematic risk rs (F) [defined by the first term in (35)] and the fluctuation or noise risk rn (F) [defined by the second term in (35)] in the desired solution (34), while the unknown disturbances of the SFO are treated through the robust WCSP optimization [constrained max problem in (34)] via imposing the statistical constraint (17) onto the bounded average squared norm of the random SFO perturbations. The selections (adjustments) of the regularization parameter α and the operator metrics inducing weight matrix A provide the additional DEDR degrees of freedom assigning the weights to the particular solution operator components. Formally, these weights can be viewed as user-defined parameters that may incorporate any descriptive metrics properties of a solution [21], [32], [36]. The ML optimal assignment α = 1, A = D(b) specifies the extended Bayesian risk [23], while in the simplest case, the uniform metrics can be induced by setting A = I. Note that different feasible assignments of the DEDR degrees of freedom α and A affect also the overall solution operator norm, thus result in different corresponding scalings of the SSP estimates (33). To avoid possible scaling mismatches, the resulting solution operators (34) can be further calibrated to some prescribed fixed norm (e.g., related to the case of the ML-optimal assignments ˆ Formally, the operator norm calibration α = 1, A = D(b)). does not influence the SSP estimation pattern; it scales the resulting image that, in practice, is always additionally rescaled to the user defined brightness representation format. We defer the detailed considerations of six practically motivated examples of different sensible assignments for α and A to the next section. To proceed with the derivation of the DEDR method, we now decompose the risk (35) incorporating directly the robust WCSP uncertainty constraint (17) into the DEDR strategy. The first term in (35) defines the averaged squared operator distance that we decompose as ˜ − I)A(FS ˜ − I)+ p(Δ) } rs (F) = tr{(FS

˜ − I 2 + FΔ 2 = FS A A p(Δ)

(36)

where C 2A = tr{CAC+ } represents the A-weighted squared norm of a matrix C. The second term on the right-hand side of (36) can next be bounded applying the Loewner ordering [21] of the weight matrix A ≤ γI with the Loewner ordering factor γ = min{γ : A ≤ γI} > 0 that yields

FΔ 2A

p(Δ)

≤ γ FΔ 2 p(Δ) ≤ γ F 2 Δ 2 p(Δ) (37)

where the second inequality follows directly from the Cauchy–Schwarz inequality [3], and C 2 = C 2I = tr{CC+ } defines a conventional squared norm of a matrix C. Next, using the inequality (17) that bounds the squared norm

90


of the SFO perturbation term, we evaluate the maximum value that may take the last term in the inequality (37)

max γ F 2 Δ 2 p(Δ) = ε F 2 (38) 2 Δ ≤η

with the bounding factor ε = ηγ ≥ 0. Under these conditions, the DEDR strategy (34) is transformed into F = arg min {RΣ (F)} F

(39)

where

(40) RΣ (F) = tr (FS − I)A(FS − I)+ +αtr FRΣ F+ represents now the aggregated robust risk function, and RΣ defines the composite noise correlation matrix RΣ = RΣ (β) = (Rn + βI)

(41)

in which the second regularization parameter (i.e., the diagonal loading factor) is given by β = ε/α = γη/α ≥ 0.

(42)

The evaluation of this error metric is detailed in Appendix C. Observe that the derived DEDR-optimal solution (44) manifests robustness against three possible problem model uncertainty factors, namely, the following. 1) Robustness against the observation noise and SFO perturbation model uncertainties: This is performed due to the employed diagonal loading in RΣ defined by (41) with the diagonal loading factor determined by (42). 2) Robustness against the numerical problem format dimension mismatch (i.e., M = K, specifically for the overdetermined representation format K > M that results in the ill-conditioned system ambiguity operator ΨΣ = S+ R−1 Σ S): This is performed due to the regularizing term αA−1 incorporated into the optimal reconstruction operator K defined by (45). 3) Robustness against the low-rank estimates of the uncertain data correlation matrix Y given by (28) and (30): This is due to the numerical structure of the DEDR estimator (44) that does not involve inversion of the matrix Y (in contrast to the conventional MVDR method and its robust extensions that require the inversion Y−1 , e.g., [20], [21], and [38]). V. F AMILY OF THE DEDR-R ELATED E STIMATORS

B. Unified DEDR Estimator Routinely solving the minimization problem (39), we obtain the DEDR-optimal solution operator F = KS+ R−1 Σ that provides via (33) the desired SSP vector estimate

ˆ DEDR = {FYF+ }diag = KS+ R−1 YR−1 SK b Σ Σ diag

(43)

(44)

where the self-adjoint matrix −1 −1 K = K(α, β, A) = S+ R−1 Σ (β)S + αA

(45)

is referred to as the DEDR-optimal reconstruction operator dependent on three degrees of freedom (two regularization parameters α, β and the diagonal weight matrix A); R−1 Σ = −1 (β) = (R + βI) , and subscript DEDR stands for the R−1 n Σ general DEDR estimator (see Appendix B where this solution is detailed). In the practical RS scenarios [2], [5], [16], [24] (and specifically, in the SAR imaging applications [4], [22], [23], [37]), it is a common practice to accept the robust white noise model, i.e., Rn = N0 I, attributing the unknown correlated noise component as well as the multiplicative speckle noise to the composite uncertain noise term Δe in (18), in which case RΣ = NΣ I with the composite noise variance NΣ = N0 + β. For the general-form solution operator (43), the resulting DEDR error measure r(F) is evaluated as a composition of the systematic error tr{(FS − I)A(FS − I)+ } and the augmented fluctuation error measure tr{FRΣ F+ }, and attains its marginal value

r(F) = tr (FS − I)A(FS − I)+ + tr{FRΣ F+ }

(46) = α2 tr{KA−1 K} + tr KS+ R−1 Σ SK .

In this section, we are going to present the examples from a family of the DEDR-related algorithms for estimating the SSP vector derived from (44) that employ different feasible adjustments of two regularization parameters α, β and the weighting matrix A (recall that these constitute the processing-level degrees of freedom of the unified DEDR method), i.e., ˆ (p) = {F(p) YF(p)+ }diag ; p = 1, . . . , P b

(47)

where different DEDR-related solution operators {F(p) ; p = 1, . . . , P } specify the corresponding RS imaging techniques. In the case of possible solution-dependent assignments, the DEDR solution operator (44) depends nonlinearly on the desired SSP ˆ hence, the nonlinear equation (47) must be next resolved b; ˆ applying some contractive iterative scheme with respect to b [3], [34]. Here beneath, we consider six (P = 6) important practically motivated examples, i.e., p = 1, . . . , 6. A. Robust Spatial Filtering Algorithm Consider, first, the white zero-mean noise in observations, hence, Rn = N0 I, the preassigned SFO uncertainty bound η and no preference to any prior model information in the objective risk function (40), i.e., putting α = 1, A = γI, with the Loewner ordering factor γ matched to the prior average gray level of the SSP b0 [23], i.e., γ = b0 , β = b0 η and NΣ = N0 + β. In this case, the solution operator F is recognized to be the Tikhonov-type robust spatial filter (RSF) FRSF = F(1) = (S+ S + αRSF I)−1 S+

(48)

i.e., a composition of the conventional matched filter given by the adjoint SFO matrix S+ and the Tikhonov-type robust reconstructive filter [23], [34] defined by the operator (S+ S +


αRSF I)−1 , in which the regularization parameter αRSF = NΣ /b0 is adjusted to the inverse to the compose signal-to-noise ratio αRSF = NΣ /b0 = (N0 /b0 ) + η. The latter is recognized to be the conventional Tikhonov’s regularization parameter N0 /b0 [19] augmented by the prescribed SFO uncertainty bound η. B. Pseudoinverse Spatial Filtering Algorithm

α→0

ˆ can be referred to as the DEDR-related Such FRASF = F(4) (b) robust extension of the conventional Wiener-type spatial filter (27), in which the initial noise correlation matrix is substituted by the diagonal loaded (i.e., robustified) RΣ given by (41), ˆ (instead of the a priori ˆ = diag(b) and the actual estimate D unknown D = diag(b)) serves as a stabilizing term in the solution operator (51). E. Conventional MVDR Algorithm

Consider now the model from the previous example but the case of α → 0, i.e., an absolute priority of the systematic error measure over the fluctuation error measure in the minimization problem (39). In this case, the solution operator becomes FPIF = F(2) = lim (S+ S + ααRSF I)− S+ = S−

91

(49)

which is essentially the pseudoinverse S− of the nominal SFO (in the Moore–Penrose sense [3]).

The relationship between the conventional MVDR beamforming algorithm and the adaptive spatial filtering method in the certain operational scenario was established in the previous studies [19], [22], [23] that we outline here as follows. The ASF of [23] is a simplified version of the specified above RASF for the case of no assumed model uncertainties in the ˜ = S, η = 0) and the full-rank unstructured estimate data (i.e., S (30) of the empirical data correlation matrix Y that satisfies (asymptotically, for J K) the correlation data agreement model [19], [23] ˆ = SD(b)S ˆ + + Rn . Y ≈ Ru (b)

C. MSF Algorithm Consider, next, the model from the first example but the alternative assumption α S+ S , i.e., a dominating priority of the second error measure (suppression of noise) over the systematic error minimization in the optimization problem (39). Under these conditions, the F is reduced to the scaled conjugate transpose (adjoint) SFO matrix, i.e., the MSF operator in the signal processing terminology [2], [4] (3)

FMSF = F

+

= cMSF S

(50)

where the scaling constant cMSF provides balance of the matrix norms c2MSF tr{S+ SS+ S} = tr{FDEDR SS+ F+ DEDR }. This constant is irrelevant for practical MSF estimators as it affects only the overall image calibration via scaling the mean gray level b0 in the resulting SSP estimate, i.e., it does not affect the fine spatial pattern of the SSP over the representation format (pixel frame) [22], [23]. D. Robust Adaptive Spatial Filtering Algorithm Consider the case of zero-mean additive observation noise with an arbitrary given positive definite invertible correlation matrix Rn , uncertain SFO with the bounded second-order moment (17) specified by the predetermined uncertainty factor η, equal importance of the systematic and fluctuation error measures in the aggregated objective risk function (35), i.e., α = 1, and the ML optimal (solution-dependent) adjustment ˆ Under these ˆ = diag(b). of the weight matrix [23] A = D conditions, the solution operator (43) becomes the Wiener-type robust ASF (RASF) ˆ = (S+ R−1 S + D ˆ −1 )−1 S+ R−1 FRASF = F(4) (b) Σ Σ

(51)

that defines via (33) the corresponding RASF estimator

ˆ RASF = FRASF YF+ b RASF diag .

(52)

(53)

Under these specifications, the solution operator (51) becomes the nonrobust (against the SFO model uncertainties) ASF operator (5) ˆ ˆ −1 −1 + −1 (b) FASF = (S+ R−1 n S + D ) S Rn = F

(54)

that admits the alternative representation form [23] ˆ ˆ + Y−1 = F(6) (b). FASF = DS

(55)

In turn, for the same regular SFO S and the invertible (fullrank) correlation data matrix estimate Y, the conventional MVDR method yields the celebrated estimator [4], [20], [27] ˆbMVDR k = s+ Y−1 sk −1 ; k

k = 1, . . . , K.

(56)

In the MVDR terminology [4], [20], [30], sk represents the socalled steering vector in the kth look direction, which in our notational conventions is essentially the kth column vector of the regular SFO S referred to as the steering matrix. In [22] and [23], it was shown that this familiar MVDR algorithm (56) directly relates to the ASF estimator ˆ ASF = {F(6) YF(6)+ }diag b ˆ + Y−1 YY−1 SD} ˆ diag = {DS −1 + −1 = { diag {S Y S}diag }diag

(57)

which in the component form coincides with (56), i.e., is algebraically equivalent to the conventional MVDR algorithm, ˆ MVDR (for the certain scenario). ˆ ASF = b i.e., b F. DEDR-Robustified MVDR Algorithm The DEDR-related robustification of the ASF algorithm (57) for the uncertain operational scenario can now be performed via substituting the ASF solution operator F(6) by the corresponding RASF operator F(4) (51) from the specified above DEDR

92


family. This yields the DEDR-extended robust version of the MVDR beamforming algorithm that we refer to as the RMVDR ˆ RASF = {F(4) YF(4)+ }diag b

−1 = KS+ R−1 Σ YRΣ SK diag −1 = { diag {Z+ Y−1 Z}diag }diag ˆ RMVDR =b

(58)

with the DEDR-modified steering matrix Z. The latter is to be determined from the solution to the −1 −1 matrix equation {KS+ R−1 Σ YY YRΣ SK}diag = + −1 −1 {[diag({Z Y Z}diag )] }diag that directly follows from (58). Note that for the RS applications, e.g., an imaging SAR considered in the companion paper [1], the RASF form ˆ RASF = {F(4) YF(4)+ }diag , is of the algorithm (58), i.e., b adopted because it does not involve the inversion Y−1 of the empirical correlation matrix Y required for implementing the RMVDR method specified by second row of (58), i.e., ˆ RMVDR = {[diag({Z+ Y−1 Z}diag )]−1 }diag . In addition, b the RASF algorithm (52) avoids numerical derivation of the solution-dependent steering matrix Z employed in the RMVDR version (58) of the robustified MVDR. Resuming the undertaken analysis, we may deduce that by specifying the corresponding degrees of freedom in the DEDR algorithmic family (47), one can proceed from the unified general estimator (44) to a variety of different DEDR-related algorithms, e.g., from the simplest MSF to the sophisticated RMVDR and RASF techniques. VI. POCS R EGULARIZATION Nonlinearity of the adaptive DEDR-related estimators implies the nonlinear dependence of the corresponding solution ˆ on the desired estimate b; ˆ thus, once we operators {F(p) (b)} have chosen the particular solution-dependent estimator from the DEDR family (47), we next need an adaptive algorithm to ˆ in some iterative form. The crucial find the desired solution b issue is to construct a convergent iterative scheme. This requires specific convex projection-form regularization [3], [39]–[41]. The idea of such regularization is to incorporate into the fixedpoint iteration method [3], [40] POCS that force the resulting POCS-regularized iteration rule to converge to a fixed point in the intersection of the sets specified by the employed POCS operators. In the optimization community, this general result is known as fundamental theorem of POCS [3, Sec. 15.4.5]. Implementation of this very general iterative convex optimization tool in different scenarios will require corresponding specifications of the POCS operators and the design of the relevant POCS-regularized fixed-point iteration rules [3], [40], [41]. We defer the design of the practical schemes particularly adapted for enhancement of the fractional SAR imagery in (near) real computational time as well as the simulations and analysis of their performances for the second companion paper [1]. Here, we resume the theoretical study with the general DEDR-related POCS regularization considerations. The desired SSP estimate must satisfy certain intrinsic RS image properties such as positivity, continuity, and agreement with the data [3], [9], [42], [43]. Observe that all DEDR-related

Fig. 3.

Illustration of a fixed-point algorithm in 1-D.

Fig. 4.

Complete framework of the DEDR method.

estimators from the family (47), by construction, provide the resulting solutions in the set of nonnegative vectors that is a convex set. When implemented iteratively, the DEDR method must incorporate (in addition to the positivity constraint) also some properly constructed convergence enforcing POCSregularizing operator P at each iteration step {i = 0, 1, . . .} to avoid numerical instability of the corresponding fixed-point iterated procedures [40], [41]. Thus, to convert the general-form DEDR estimator (44) to an iterative algorithm, one has to define a sequence of estimates

ˆ [i] ) = P K[i] S+ R−1 YR−1 SK[i] ˆ [i+1] = T [i] (b (59) b Σ Σ diag i = 0, 1, . . ., in which the POCS-regularized fixed-point operator T[i] (·) is defined by the right-hand side of (59), and ˆ [i] ) = ((S+ R−1 S + αdiag−1 (b ˆ [i] ))−1 K[i] = K(b Σ

(60)

represents the reconstruction operator adaptively updated at each iteration step. A graphical way of illustrating the fixedpoint algorithm (59) is shown in Fig. 3. General framework of the POCS-regularized DEDR method (DEDR-POCS) is next shown in Fig. 4. In general, the incorporation of the POCS regularization enforces the convergence of (59) but computationally intensive solution-dependent matrix inversions required to compute the reconstruction operator (60) at each iteration step complicates severely the numerical implementation of the overall adaptive DEDR procedure. As a resuming example, consider a typical imaging SAR experiment with M collected rangeazimuth trajectory signal samples and range-azimuth scene dimension K = Kx × Ky . Assume that numerical arithmetic with individual element has complexity O(1). Then, the overall


computational complexity of the general-form fixed-point iteration process (59) is of the order ∼ O((K × M 2 ) + (K 2 + K 3 ) × I) for I executed fixed-point iterations (see companion paper [1, Sec. V-B] for the complexity evaluation details). Because of high problem dimensionalities, this is an extremely intensive computational procedure [40], [41]; hence, it is questionable to address the general-form fixed-point DEDR technique (59), (60) as a practical SSP estimator executable in (near) real processing time. In the companion paper [1], we develop the computationally efficient POCS-regularized DEDR-related iterative algorithms particularly adapted to the imaging SAR applications along with their simulations and performance analysis. It may, therefore, also be a good starting point for those more interested in numerical algorithms than in theoretical derivations. VII. C ONCLUDING R EMARKS In this paper, a DEDR approach for estimation of the SSP of the scattered wavefield power distribution in the uncertain remotely sensed environment has been proposed as required for an imaging radar and SAR. With the unified DEDR method, the SSP reconstruction is performed in a descriptively balanced fashion via weighted maximization of the spatial resolution over minimization of the resulting noise energy algorithmically coupled with the WCSP optimization-based regularization. The family of the DEDR-related SSP estimators was constructed for different typical model assumptions regarding the uncertain operational scenarios with the statistical observation model mismatches and ill-conditioned (low-rank) estimates of the data correlation matrices. Being nonlinear and solution dependent, the general-form unified DEDR method requires rather complex signal processing operations ruled by the fixed-point iterative implementation process. To reduce the computational load and derive the practically attractive application algorithms, we are now ready to proceed with the application part of this research in the companion paper [1]. A PPENDIX A ML E STIMATORS FOR C ERTAIN S CENARIO This section provides a derivation of the ML estimators of the SSP in the operational scenario without model uncertainties. Here, we use the following notational conventions: D(a) = diag(a) diagonal matrix, where a is a vector of the principal diagonal of matrix D; vector composed of the elements of a = {A}diag the principal diagonal of the embraced matrix A; kjth element of the embraced matrix A. {A}kj The objective minimization function is defined by (25) in this paper, i.e.,

−1 L(b) = ln det SD(b)S+ +Rn +[u, SD(b)S+ +Rn u] (A1) where u = Se + n defines the single observed zero-mean random realization of the data vector with a given noise correlation

93

matrix Rn . The inner product for the complex-valued vectors in (A1) is defined as T −1 ∗ u, R−1 (A2) u u = u Ru u where Ru represents the correlation matrix of the data vector u given by Ru = SD(b)S+ + Rn .

(A3)

The inverse of such composite matrix Ru is [6] + −1 −1 −1 R−1 u = Rn − Rn SKS Rn

(A4)

−1 −1 . K = K(b) = S+ R−1 n S + D (b)

(A5)

where

ˆ ML = With these specifications, the ML optimization problem b arg min{L(b)} is structurally similar to that considered in the b

previous study [5]. Thus, borrowing from [5, Appendix A], we first specify the ML-optimal solution operator for the case of the single processed realization of the observation data FW = DS+ R−1 u

(A6)

that admits also the second representation form [5] + −1 −1 −1 + −1 S Rn = FW . DS+ R−1 u = S Rn S + D

(A7)

Observe that this solution operator, FW is recognized to be the ˆ ML )) regularized Wienersolution-dependent (i.e., D = diag(b type inverse [23] of the regular SFO S that yields the MLoptimal SSP estimate

ˆ ML = FW uu+ F+ (A8) b W diag . The extension of this ML estimator (A8) for the case of arbitrary recorded data observations is performed via averaging the correlations over J actually registered independent realizations of {u(j) } [5] that yields

ˆ ML = {R ˆ e }diag = FW YF+ b W diag

(A9)

where Y represents now the empirical estimate of the data correlation matrix [10], [20] ˆ u = aver{u(j) u+ } = (1/J) Y=R (j) j∈J

J

u(j) u+ (j) .

(A10)

j=1

Observe that (A8) and (A9) correspond to the DEDR strategy for the simplified scenario of no model uncertainties in the SFO and in the data acquisition process. A PPENDIX B D ERIVATION OF THE G ENERAL -F ORM DEDR S OLUTION O PERATOR (43) The problem to be resolved in this section is to derive the solution operator that is optimal in the sense of the DEDR strategy (39). Such problem is a direct expansion of the one considered

94


in the previous study [19]. The expansion is performed via the generalization of the aggregated risk function

i.e., the desired general-form DEDR solution operator defined by (43) in this paper.

F = arg min {RΣ (F)} F

= arg min tr (FS−I)A(FS−I)+ +αtr FRΣ F+ .

A PPENDIX C E VALUATION OF THE DEDR E RROR M EASURE (46)

F

(B1)

In this section, we evaluate the performance measure

To determine the optimal solution operator, we have to differentiate this objective risk function RΣ (F) with respect to F, set the result to zero, and solve the corresponding variational equation. To proceed with the calculations, we first decompose the first addend in the aggregated risk function using the formula FSA(FS)+ = FSAS+ F+ , and rewrite (B1) as follows:

F = arg min tr FSAS+ F+ − trFSA + tr AS+ F+ F

+ tr{A} + αtr{FRΣ F+ } . (B2)

r(F) = tr{(FS − I)A(FS − I)+ } + tr{FRΣ F+ }

Next, we invoke the formulas for differentiating the traces of the composition of matrices with respect to a matrix [19] ∂tr{FCF+ } = 2FC ∂F

∂ {tr{FT} + tr{T+ F+ }} = 2T+ . ∂F (B3)

To apply these formulas for solving the minimization problem (B2), we associate C with SAS+ for the first addend from (B2) and with RΣ for the last addend from (B2), correspondingly, while T is associated with SA. In addition, in calculations, we take into account that the A is a self-adjoint real-valued diagonal weighting matrix (metrics inducing operator), i.e., A = A+ , hence T+ = AS+ . Following these specifications, we apply now formulas (B3) to (B2) to get the expression for the matrix derivative ∂RΣ (F)/∂F and then set the result to zero. This yields the following matrix-form variational equation: ∂RΣ (F )}/∂F = 2FSAS+ − 2AS+ + 2αFRΣ = 2(FS − I)AS+ + 2αFRΣ = 0.

(B4)

(B5)

that yields the desired solution operator in its initial form F = AS+ (SAS+ + αRΣ )−1 .

(B6)

Next, we make the use of another algebraically equivalent form of representation of the matrix composition F = AS+ (SAS+ + αRΣ )−1 −1 −1 + −1 = (S+ R−1 Σ S + αA ) S RΣ

(B7)

detailed, for example, in [5, Appendix B] that results in the desired solution operator F = FDEDR = KS+ R−1 Σ

(B8)

with −1 −1 K = K(α, β, A) = S+ R−1 Σ (β)S + αA

attained at the DEDR-optimal solution operator (B8). This performance metric is a composition of the systematic error rs (F) = tr{(FS − I)A(FS − I)+ } and the augmented fluctuation (or noise) error measure rn (F) = tr{FRΣ F+ }. First, we evaluate the systematic error term

rs (F) = tr (FS − I)A(FS − I)+

= tr (KΨΣ − I)A(KΨΣ − I)+

= tr K(ΨΣ − K−1 )A(ΨΣ − K−1 )K = tr{KαA−1 AαA−1 K} (C2) = α2 tr{KA−1 K} where ΨΣ = S+ R−1 Σ S

(C3)

is referred to as the robust matrix-form ambiguity operator. Second, we evaluate the fluctuation error measure defined by the second addend in (C1) using (B7) to express the solution operator F that yields + } rn (F) = tr{FR

Σ+F −1 = tr KS RΣ RΣ R−1 Σ SK = tr{KΨΣ K}.

(C4)

Next, we aggregate these two measures (C2) and (C4) and obtain r(F)= rs (F) + rn (F) = α2 tr{KA−1 K}+tr{KΨΣ K} (C5)

Rearranging (B4), we obtain F(SAS+ + αRΣ ) = AS+

(C1)

(B9)

that corresponds to (46) in this paper. Last, we specify this quality metric for the algebraically equivalent RASF and RMVDR estimators (58) that optimally balance the systematic and fluctuation error measures, i.e., for ˆ in which case the the ML-optimal assignments αA = diag(b), error measure (C5) attains its minimal marginal value r(FRASF ) = r(FMVDR ) = tr{KA−1 K} + tr{KΨΣ K} = tr{K}. (C6) ACKNOWLEDGMENT The author would like to thank the anonymous reviewers who pointed out a number of suggested improvements. R EFERENCES [1] Y. V. Shkvarko, “Unifying experiment design and convex regularization techniques for enhanced imaging with uncertain remote sensing data— Part II: Adaptive implementation and performance issues,” IEEE Trans. Geosci. Remote Sens., vol. 48, no. 1, pp. 83–98, Jan. 2010. [2] F. M. Henderson and A. V. Lewis, Eds., Principles and Applications of Imaging Radar, Manual of Remote Sensing, 3rd ed., vol. 3. New York: Wiley, 1998.


[3] H. H. Barrett and K. J. Myers, Foundations of Image Science. New York: Wiley, 2004. [4] A. Farina, Antenna-Based Signal Processing Techniques for Radar Systems. Norwood, MA: Artech House, 1991. [5] Y. V. Shkvarko, “Estimation of wavefield power distribution in the remotely sensed environment: Bayesian maximum entropy approach,” IEEE Trans. Signal Process., vol. 50, no. 9, pp. 2333–2346, Sep. 2002. [6] Y. V. Shkvarko, “Unifying regularization and Bayesian estimation methods for enhanced imaging with remotely sensed data—Part I: Theory,” IEEE Trans. Geosci. Remote Sens., vol. 42, no. 5, pp. 923–931, May 2004. [7] Y. V. Shkvarko, “Unifying regularization and Bayesian estimation methods for enhanced imaging with remotely sensed data—Part II: Implementation and performance issues,” IEEE Trans. Geosci. Remote Sens., vol. 42, no. 5, pp. 932–940, May 2004. [8] L. E. Franks, Signal Theory. Englewood Cliffs, NJ: Prentice-Hall, 1969. [9] J. R. Fienup, “Phase retrieval algorithms: A comparison,” Appl. Opt., vol. 21, no. 15, pp. 2758–2769, Aug. 1982. [10] Adaptive Radar Detection and Estimation, S. Haykin and A. Steinhardt, Eds. New York: Wiley, 1992. [11] D. R. Wehner, High-Resolution Radar, 2nd ed. Boston, MA: Artech House, 1994. [12] A. Levi and H. Stark, “Image restoration by the method of generalized projections with application to restoration from magnitude,” J. Opt. Soc. Amer. A, Opt. Image Sci., vol. 1, no. 9, pp. 932–943, Sep. 1984. [13] S. W. Perry, H.-S. Wong, and L. Guan, Adaptive Imaging: A Computational Intelligence Perspective. Boca Raton, FL: CRC Press, 2001. [14] Y. V. Shkvarko and V. L. Kolyadin, “Spatial spread object scattering function restoration using the restricted radio data set of multichannel measurements,” Sov. J. Commun. Technol. Electron., vol. 33, no. 6, pp. 1174– 1181, Jun. 1988. English translation of Radiotekhnika i Elektronika. [15] J. R. Fienup, “Synthetic-aperture radar autofocus by maximizing sharpness,” Opt. Lett., vol. 25, no. 4, pp. 221–223, Feb. 2000. [16] D. C. Bell and R. M. Narayanan, “Theoretical aspects of radar imaging using stochastic waveforms,” IEEE Trans. Signal Process., vol. 49, no. 2, pp. 394–400, Feb. 2001. [17] Y. V. Shkvarko, Y. S. Shmaliy, R. Jaime-Rivas, and M. Torres-Cisneros, “System fusion in passive sensing using a modified Hopfield network,” J. Franklin Inst., vol. 338, no. 4, pp. 405–427, Jul. 2001. [18] J. S. Lee, “Speckle suppression and analysis for synthetic aperture radar images,” Opt. Eng., vol. 25, no. 5, pp. 636–643, May 1986. [19] Y. V. Shkvarko, “From matched spatial filtering towards the fused statistical descriptive regularization method for enhanced radar imaging,” EURASIP J. Appl. Signal Process., vol. 2006, pp. 1–9, 2006. [20] A. B. Gershman, “Robustness issues in adaptive beamforming and highresolution direction finding,” in High-Resolution and Robust Signal Processing, Y. Hua, A. B. Gershman, and Q. Cheng, Eds. New York: Marcel Dekker, 2004, pp. 63–100. [21] Y. Shkvarko, H. Perez-Meana, and A. Castillo-Atoche, “Enhanced radar imaging in uncertain environment: A descriptive experiment design regularization approach,” Int. J. Navig. Obs., vol. 2008, pp. 1–11, 2008. [22] Y. V. Shkvarko, J. L. Leyva-Montiel, and I. E. Villalon-Turrubiates, “Unifying the experiment design and constrained regularization paradigms for reconstructive imaging with remote sensing data,” in Proc. IEEE SP Int. Conf. Image Process., Atlanta, GA, 2006, pp. 3241–3244. [23] Y. V. Shkvarko, “Finite array observations-adapted regularization unified with descriptive experiment design approach for high-resolution spatial power spectrum estimation with application to radar/SAR imaging,” in Proc. 15th IEEE SP Int. Conf. Digital Signal Process., Cardiff, U.K., 2007, pp. 79–82. [24] A. Ishimary, Wave Propagation and Scattering in Random Media. New York: IEEE Press, 1997. [25] R. L. Fante, “Electromagnetic beam propagation in turbulent media,” Proc. IEEE, vol. 63, no. 12, pp. 1669–1692, Dec. 1975. [26] R. L. Fante, “Turbulence-induced distortion of synthetic aperture radar images,” IEEE Trans. Geosci. Remote Sens., vol. 32, no. 4, pp. 958–961, Jul. 1994. [27] Z. W. Xu, J. Wu, and Z. S. Wu, “Potential effects of the ionosphere on space-based SAR imaging,” IEEE Trans. Antennas Propag., vol. 56, no. 7, pp. 1968–1975, Jul. 2008. [28] G. Franceschetti, A. Iodice, S. Perna, and D. Riccio, “SAR sensor trajectory deviations: Fourier domain formulation and extended scene simulation of raw signal,” IEEE Trans. Geosci. Remote Sens., vol. 44, no. 9, pp. 2323–2334, Sep. 2006.

95

[29] G. Franceschetti, A. Iodice, S. Perna, and D. Riccio, “Efficient simulation of airborne SAR raw data of extended scenes,” IEEE Trans. Geosci. Remote Sens., vol. 44, no. 10, pp. 2851–2860, Oct. 2006. [30] M. Ruegg, E. Meier, and D. Nuesch, “Vibration and rotation in millimeterwave SAR,” IEEE Trans. Geosci. Remote Sens., vol. 45, no. 2, pp. 293– 304, Feb. 2007. [31] M. S. Greco and F. Gini, “Statistical analysis of high-resolution SAR ground clutter data,” IEEE Trans. Geosci. Remote Sens., vol. 45, no. 3, pp. 566–575, Mar. 2007. [32] V. Ponomaryov, A. Rosales, F. Gallegos, and I. Loboda, “Adaptive vector directional filters to process multichannel images,” IEICE Trans. Commun., vol. E90-B, no. 2, pp. 429–430, Feb. 2007. [33] M. C. Wicks, M. Rangaswamy, R. Adve, and T. B. Hale, “Space–time adaptive processing: A knowledge-based perspective for airborne radar,” IEEE Signal Process. Mag., vol. 23, no. 1, pp. 51–65, Jan. 2006. [34] A. N. Tikhonov and V. Y. Arsenin, Solution of Ill-Posed Problems. Washington, DC: W. H. Winston, 1977. [35] S. M. Kay, Modern Spectral Estimation. Englewood Cliffs, NJ: Prentice-Hall, 1988. [36] M. Younis, C. Fisher, and W. Wiesbeck, “Digital beamforming in SAR systems,” IEEE Trans. Geosci. Remote Sens., vol. 41, no. 7, pp. 1735– 1739, Jul. 2003. [37] F. De Zan and A. Monti Guarnieri, “TOPSAR: Terrain observation by progressive scans,” IEEE Trans. Geosci. Remote Sens., vol. 44, no. 9, pp. 2352–2360, Sep. 2006. [38] A. Hassanien, S. Shahbazpanahi, and A. B. Gershman, “A general Capon estimator for localization of multiple spread sources,” IEEE Trans. Signal Process., vol. 52, no. 1, pp. 280–282, Jan. 2004. [39] Y. Shkvarko, J. A. Gutierrez-Rosas, and L. G. Guerrero Diaz de Leon, “Towards the virtual remote sensing laboratory: Simulation software for intelligent post-processing of large scale remote sensing imagery,” in Proc. IEEE IGARSS, Barcelona, Spain, Jul. 2007, pp. 1561–1564. [40] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C: The Art of Scientific Computing, 2nd ed. New York: Cambridge Univ. Press, 1992. [41] J. H. Mathews, Numerical Methods for Mathematics, Science, and Engineering, 2nd ed. Englewood Cliffs, NJ: Prentice-Hall, 1992. [42] O. Loffeld, H. Niez, S. Knedlik, and W. Yu, “Phase unwrapping for SAR interferometry—A data fusion approach by Kalman filtering,” IEEE Trans. Geosci. Remote Sens., vol. 46, no. 1, pp. 47–58, Jan. 2008. [43] Y. Gu, Y. Zhang, and J. Zhang, “Integration of spatial spectral information for resolution enhancement in hyperspectral images,” IEEE Trans. Geosci. Remote Sens., vol. 46, no. 5, pp. 1347–1358, May 2008. [44] Z. Li, S. Papson, and R. M. Narayanan, “Data-level fusion of multilook inverse synthetic aperture radar images,” IEEE Trans. Geosci. Remote Sens., vol. 46, no. 5, pp. 1394–1406, May 2008.

Yuriy V. Shkvarko (M’95–SM’04) received the Dip.Eng.(Hon.) degree in electrical engineering, the Ph.D. degree in radio engineering, and the Dr. Sci. degree in radio physics, radar, and navigation from the Kharkov Aviation Institute, Kharkov, Ukraine, in 1976, 1980, and 1990, respectively. From 1976 to 1991, he was with the Scientific Research Department, Kharkov Aviation Institute, as a Research Fellow, a Senior Fellow, and Chair of the Research Laboratory in information technologies for radar and navigation. From 1991 to 1999, he was a Full Professor with the Department of System Analysis and Control, Ukrainian National Polytechnic Institute, Kharkov. He immigrated to Mexico in 1999. From 1999 to 2001, he was an Invited Professor with the Guanajuato State University, Salamanca, Mexico. Since 2001, he has been with the Centro de Investigación y de Estudios Avanzados del Instituto Politécnico Nacional (Superior Education and Research Center of the National Polytechnic Institute), Guadalajara, Mexico, as a Full Titular Professor. He is the holder of 12 patents and has published two books and some 150 journal and conference papers. His research interests are in applications of signal processing to remote sensing, imaging radar, navigation, and communications, particularly in inverse problems, random fields estimation, adaptive spatial analysis, statistical sensor array and multichannel remote sensing data processing, and system fusion.