DTIC ADA580551: Quantitative Image Recovery From Measured Blind Backscattered Data Using a Globally Convergent Inverse Method

3.8 MB·English

by  Defense Technical Information Center

REPORT DOCUMENTATION PAGE Form Approved OMB NO. 0704-0188
The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions,
searching existing data sources, gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments
regarding this burden estimate or any other aspect of this collection of information, including suggesstions for reducing this burden, to Washington
Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington VA, 22202-4302.
Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any oenalty for failing to comply with a collection of
information if it does not display a currently valid OMB control number.
PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.
1. REPORT DATE (DD-MM-YYYY) 2. REPORT TYPE 3. DATES COVERED (From - To)
New Reprint -
4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER
Quantitative image recovery from measured blind backscattered W911NF-11-1-0399
data using a globally convergent inverse method
5b. GRANT NUMBER
5c. PROGRAM ELEMENT NUMBER
611102
6. AUTHORS 5d. PROJECT NUMBER
Andrey V. Kuzhuget, Larisa Beilina, Michael V. Klibanov, Anders
Sullivan, Lam Nguyen, Michael A. Fiddy 5e. TASK NUMBER
5f. WORK UNIT NUMBER
7. PERFORMING ORGANIZATION NAMES AND ADDRESSES 8. PERFORMING ORGANIZATION REPORT
NUMBER
University of North Carolina - Charlotte
Office of Sponsored Programs
9201 University City Blvd.
Charlotte, NC 28223 -0001
9. SPONSORING/MONITORING AGENCY NAME(S) AND 10. SPONSOR/MONITOR'S ACRONYM(S)
ADDRESS(ES) ARO
U.S. Army Research Office 11. SPONSOR/MONITOR'S REPORT
P.O. Box 12211 NUMBER(S)
Research Triangle Park, NC 27709-2211 60035-MA.6
12. DISTRIBUTION AVAILIBILITY STATEMENT
Approved for public release; distribution is unlimited.
13. SUPPLEMENTARY NOTES
The views, opinions and/or findings contained in this report are those of the author(s) and should not contrued as an official Department
of the Army position, policy or decision, unless so designated by other documentation.
14. ABSTRACT
The goal of this paper is to introduce the application of a globally convergent inverse scattering algorithm to
estimate dielectric constants of targets using time resolved backscattering data collected by a US Army Research
Laboratory (ARL) forward looking radar. The processing of the data was conducted blind, i.e. without any prior
knowledge of the targets. The problem is solved by formulating the scattering problem as a coefficient inverse
problem for a hyperbolic partial differential equation. The main new feature of this algorithm is its rigorously
15. SUBJECT TERMS
joint paper with ARL engineers Anders Sullivan and Lam Nguyen, remote sensing, inverse scattering, quantitative imaging,
experimental data
16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF 15. NUMBER 19a. NAME OF RESPONSIBLE PERSON
a. REPORT b. ABSTRACT c. THIS PAGE ABSTRACT OF PAGES Michael Klibanov
UU UU UU UU 19b. TELEPHONE NUMBER
704-687-2645
Standard Form 298 (Rev 8/98)
Prescribed by ANSI Std. Z39.18
Report Title
Quantitative image recovery from measured blind backscattered data using a globally convergent inverse method
ABSTRACT
The goal of this paper is to introduce the application of a globally convergent inverse scattering algorithm to estimate
dielectric constants of targets using time resolved backscattering data collected by a US Army Research Laboratory
(ARL) forward looking radar. The processing of the data was conducted blind, i.e. without any prior knowledge of
the targets. The problem is solved by formulating the scattering problem as a coefficient inverse problem for a
hyperbolic partial differential equation. The main new feature of this algorithm is its rigorously established global
convergence property. Calculated values of dielectric constants are in a good agreement with material properties,
which were revealed a posteriori.
REPORT DOCUMENTATION PAGE (SF298)
(Continuation Sheet)
Continuation for Block 13
ARO Report Number 60035.6-MA
Quantitative image recovery from measured blin ...
Block 13: Supplementary Note
© 2012 . Published in IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, Vol. Ed. 0 (2012), (Ed. ). DoD
Components reserve a royalty-free, nonexclusive and irrevocable right to reproduce, publish, or otherwise use the work for
Federal purposes, and to authroize others to do so (DODGARS §32.36). The views, opinions and/or findings contained in this
report are those of the author(s) and should not be construed as an official Department of the Army position, policy or decision,
unless so designated by other documentation.
Approved for public release; distribution is unlimited.
IEEETRANSACTIONSONGEOSCIENCEANDREMOTESENSING 1
Quantitative Image Recovery From Measured Blind
1
Backscattered Data Using a Globally
2
Convergent Inverse Method
3
4 AndreyV.Kuzhuget,LarisaBeilina,MichaelV.Klibanov,AndersSullivan,LamNguyen,andMichaelA.Fiddy
5 Abstract—Thegoalofthispaperistointroducetheapplication from only a single location of either a point source or from 34
6of a globally convergent inverse scattering algorithm to estimate a single direction of an incident plane wave. In particular, in35
7dielectricconstantsoftargetsusingtime-resolvedbackscattering [14], that method was extended from the 3-D case to the 1-D36
8data collected by a U.S. Army Research Laboratory forward-
case. Thus, that 1-D version of [14] is used here to work with37
9lookingradar.Theprocessingofthedatawasconductedblind,i.e.,
10withoutanypriorknowledgeofthetargets.Theproblemissolved theexperimentaldata.Theilluminatingfieldispulsedintime, 38
11by formulating the scattering problem as a coefficient inverse andthetEimehistoryofthebackscatteringfromtheilluminated39
12problemforahyperbolicpartialdifferentialequation.Themain targetvolumeconstitutesthemeasureddatathatareprocessed 40
13new feature of this algorithm is its rigorously established global by this algorithm. The authors are unaware of other groups 41
14convergenceproperty.Calculatedvaluesofdielectricconstantsare
working on MCIPs using data acquired from a single source42
15inagoodagreementwithmaterialproperties,whichwererevealed
16aposteriori. location. However, the single measurement case is clearly the43
most practical one, particularly for military applications. In-44
17 IndexTerms—Experimentaldata,inversescattering,quantita-
Edeed, because of many dangers on the battlefield, the number 45
18tiveimaging,remotesensing.
ofmeasurementsshouldbeminimized.f 46
19 I. INTRODUCTION Thealgorithmintheaforementionedcitedpublicationscom- 47
putes values for the spatial distribution of the dielectric con- 48
20A FUNDAMENTAL problem in remote sensing is the stants of objects within the targeot volume. It is important to49
21 processingofscatteredfielddatafromstronglyscattering stress that this algorithm requires neither no prior knowledge 50
22penetrabletargets.MultiplescatteringrendersthEisproblemex- ofwhatmightexistinthetargetvolumenorapriorknowledge 51
23tremelydifficulttosolve,itbeingillconditionedwithadditional ofagoodfirstguessaboutthesolution.Thereisarigorousguar-52
24questionsofuniquenessand,themostdifficult,nonlinearityto anteethatthisalgorithmgloballyconverges(seemathematical 53
o
25contend with. In practice, limited noisy data typically require detailsin[7],[14],[16],and[17]).Becauseoftheglobalcon-54
26thatsomephysicalmodelsbeassumed,fromwhichonehopes vergence property, estimates of spatially distributed dielectric 55
27to extract meaningful and preferably quantitative information constants are reliable and systematically improve with more56
I
28about the target in question. A number of recent publications measuredandlessnoisydata.Thetheoryoftheaforementioned 57
29by Beilina and Klibanov [3]–[8] and by Klibanov et al. [12], cited publicrations rigorously guarantees that this numerical 58
30[14]–[16]haveledtoanewapproachtoaddressthisimportant method delivers a good approximation to the exact solution59
31topic. This numerical method was originally developed for of an MCIP without any a priori information about a small60
32some multidimensional coefficient inverse problems (MCIPs) neighbPorhood of the exact solution as long as iterations start61
33for a hyperbolic partial differential equation (PDE) using data from the so-called “first tail function” V0(x), which can be62
easily computed using available boundary measurements (see63
Manuscript received March 24, 2012; revised July 22, 2012; accepted (2.27)–(2.29) in Section II-C). In addition, it is in this sense64
July27,2012.ThisworkwassupportedinpartbytheU.S.ArmyResearch that we use the term “global convergence" of the algorithm.65
LaboratoryandtheU.S.ArmyResearchOfficeunderGrantW911NF-11-1- The common perception of the term "global convergence" is66
0399;bytheSwedishResearchCouncil(VR);bytheSwedishFoundationfor
StrategicResearch(SSF)inGothenburgMathematicalModellingCentre;and thatonecanstartfromany pointand stillgetthesolution,but 67
bytheSwedishInstitute,VisbyProgram. we stress that we actually start not from any point but rather68
AQ1 AL..BVe.iKliunzahiusgweittihstwheithDMepoarrgtmanenSttaonflMeya,tNheemwaYtiocrakl,SNciYen1c0e0s,36ChUaSlmAe.rsUni- from the function V0(x), which can be easily computed from69
versityofTechnology,42196Gothenburg,Sweden,andalsowithGothenburg theboundarydata(see(2.27)–(2.29)inSectionII-C). 70
University,40530Gothenburg,Sweden(e-mail:[email protected]). It is well known that least squares functionals for MCIPs71
M. V. Klibanov is with the Department of Mathematics and Statistics,
suffer from multiple local minima and ravines. Hence, local 72
UniversityofNorthCarolinaatCharlotte,Charlotte,NC28223USA(e-mail:
[email protected]). convergence of numerical methods to incorrect estimates will73
A. Sullivan and L. Nguyen are with the U.S. Army Research Labora- occurunlessaninitialguessthatisclosetothetruesolutionis74
tory,Adelphi,MD20783USA(e-mail:[email protected];lam.h.
used. Such a guess is rarely available in most applications. In75
[email protected]).
M.A.FiddyiswiththeOptoelectronicsCenter,UniversityofNorthCarolina contrast,ouralgorithmdoesnotusealeastsquaresfunctional, 76
atCharlotte,Charlotte,NC28223USA(e-mail:mafi[email protected]). andhence,itisfreefromtheproblemoflocalminima.Instead, 77
Colorversionsofoneormoreofthefiguresinthispaperareavailableonline
thisalgorithmreliesonthestructureofthedifferentialoperator78
athttp://ieeexplore.ieee.org.
DigitalObjectIdentifier10.1109/TGRS.2012.2211885 ofthewave-likePDE. 79
0196-2892/$31.00©2012IEEE
2 IEEETRANSACTIONSONGEOSCIENCEANDREMOTESENSING
80 Prior to the work reported here, a major focus by the We assume that the constitutive parameter of interest, i.e.,135
81U.S. Army Research Laboratory (ARL) had been on the de- mapping the target volume, is a relative permittivity εr(x). In136
82velopment of image processing techniques [19] that would otherwords,weignoremagneticeffectsinthispaper.Wealso137
83improve radar images, which is through postprocessing tech- assume for convenience that εr(x)=1 outside of the target 138
84niquesratherthanthroughtheapplicationofinversescattering volume, which is x∈(0,1) in our case. We assume that the139
85methods. By incorporating more physics of the target-wave source x0 <0 lies outside of the target volume. We can write 140
86electromagnetic response into the data processing, one can theforwardscatteringproblemas 141
87greatlyimprovetargetdetectionandidentification.Presentdata
88processing provides an electromagnetic field brightness or an εr(x)utt =uxx, x∈R (2.1)
89intensity map of the target volume, which need not relate u(x,0)=0, ut(x,0)=δ(x−x0). (2.2)
90in a simple fashion to the scattering structures themselves.
Thesubscriptsin(2.1)indicatethenumberofpartialderivatives 142
91Our method estimates dielectric constants of targets, which
with respect to the variable indicated. The coefficient inverse143
92obviously adds an important new dimension to the interpre-
93tation of data acquired by the radar system since this allows problem (CIP) is to recover εr(x), assuming that the initial144
illuminating pulse is known and that we measure the function145
94specific bounds on the dielectric properties of a feature in
g(t),i.e., 146
95the target volume, which can help identify its likely material
E
96properties.Sincenopriorknowledgeisrequired,themeasured
u(0,t)=g(t) (2.3)
97datawereprocessedbyKuzhuget,Beilina,Klibanov,andFiddy
98in the most challenging scenario, i.e., without any knowledge for sufficiently large times t that all multiple scattering events 147
99of the actual target structures and their dielectric properties. within the target volume, which can produce a measurable 148
100Once this had been done, Sullivan and Nguyen compared a signal at the detector, do so. Practically, we gate the radiation 149
101posterioritheimageestimateswiththeactuallyknownmaterEial source in time; and since the Laplace transform (LT), i.e., 150
102characteristics. w(x,s), is used to solve this CIP, the decay e−st,s>0 of151
f
103 We draw attention to the fact that this algorithm has been the LT kernel further limits the duration of the measured time 152
104used with forward-scattered data from experiments. These history.Itisworthpointingoutthat,moretypically,scattering 153
105results were previously reported, which are also in a blind datawouldbemeasuredatdifferoentscatteringanglesforfixed 154
106experiment (see [12, Tables 5 and 6] and [7, Tables 5.5 and frequency illumination at various incident angles. One can 155
1075.6]). In this case, the images in [12] were fuErther improved easily appreciate that this leads to the acquisition of Fourier156
108andpresentedinafollow-uppublication[6]usingtheadaptivity information about the target or the secondary source function, 157
109techniqueof[1],[2],[4],[5],and[7]. dependingupontheextentofthemultiplescattering;andonce158
110 In Section II, we outline the basic steps in the underlying one has sufficient daota, a reasonable estimate of the target159
111theory upon which the new algorithm is based. In Section III, properties becomes possible. By taking measurements in the 160
112we formulate the global convergence theorem. In Section IV, time domain, one can see that this is essentially simultane-161
113we outline results obtained usinIg time-resolved backscatter ously gathering information in a transform space from many162
114electric field measurements collected in the field. Measure- illuminatiorn frequencies. The Laplace and Fourier transforms 163
115mentswerecarriedoutbyaforward-lookingradarsystembuilt providecomplimentaryrepresentationsofthetargetintermsof 164
116andoperatedbytheARL.Thedatawerenoisyandlimited,and momentsormodes,respectively. 165
117the target volumes included miscellaneous sources of clutter. ThPeLTis 166
118The purpose of this particular radar system is to detect and
(cid:1)∞
119possiblyidentifyshallowexplosive-liketargets.
w(x,s)= u(x,t)e−stdt:=Lu, s≥s=const.>0 (2.4)
0
120 II. THEORETICALBACKGROUND and we assume that the so-called pseudofrequency s≥167
121A. IntegralDifferentialEquation s(εr(x)):=sissufficientlylarge.Thisgives[7] 168
122 Since we were given only one time-resolved experimental w −s2ε (x)w = −δ(x−x ), x∈R (2.5)
xx r 0
123curvepertarget,wehadnochoicebuttousea1-Dmathemati- lim w(x,s)=0. (2.6)
124calmodel,althoughtherealityis3-D(seeSectionIIIforsome |x|→∞
125details about the data collection). In addition, since only one
Let 169
126componentoftheelectricwavefieldwasbothtransmittedand
127measured,wemodelthescatteringprocesswithonewave-like w(0,s)=ϕ(s)=Lg (2.7)
128PDEratherthanusingcompleteMaxwellequations.Westress
129that the method is designed for use with 3-D problems, and betheLTofthemeasuredfunctiong(t)in(2.3).Sinceεr(x)=170
130onewouldnormallycollectdatawithcopolarizationandcross 1for x<0,then, using (2.5)and (2.6),one can prove that,in171
131polarizationinordertocaptureallofthepertinentinformation addition to the function w(0,s) in (2.7), the function wx(0,s) 172
132about the target. Here, we simply wish to show the steps isalsoknownas(see[17]) 173
133employed by the method and demonstrate their quantitative
134reconstructionaccuracygivennoisymeasureddata. wx(0,s)=sϕ(s)−exp(sx0). (2.8)
KUZHUGETetal.:IMAGERECOVERYFROMBLINDBACKSCATTEREDDATAUSINGANINVERSEMETHOD 3
174 Letw0(x,s)bethesolutionoftheproblemin(2.5)and(2.6)
175forthecaseoftheuniformbackgroundεr(x)≡1.Then
exp(−s|x−x |)
w (x,s)= 0 . (2.9)
0 2s
176When implementing the algorithm, given the assumption of a
177uniformnormalizedεr(x)=1outsideofthetargetvolume,we
178considerthefunction
(cid:2) (cid:3)
1 w
r(x,s)= ln (x,s) . (2.10)
s2 w
0
179Sincethesourcex0 <0,thenthefunctionr(x,s)isthesolution
180ofthefollowingequationintheinterval(0,1):
r +s2r2−2sr =ε (x)−1, x∈(0,1). (2.11)
xx x x r
181Inaddition,by(2.7)and(2.8)
E
Fig.1. Flowchartofthealgorithm.
r(0,s)=ϕ (s), r (0,s)=ϕ (s) (2.12)
0 x 1
ϕ (s)=lnϕ(s)−ln(2s) + x0 wherefunctionsψ0(s)=ϕ(cid:9)0(s)andψ1(s)=ϕ(cid:9)1(s)arederived 194
0 s2 s from (2.13). The condition qx(1,s)=0 can be easily derived 195
ϕ (s)=2 − esx0 . (2.13) from(2.6)sinceεr(x)=1outsideoftheinterval(0,1). 196
1 s s2ϕ(s) In (2.17) and (2.18), both functions q(x,s) and V(x) are197
Eunknown. The reason why we can approximate both of them 198
118823froTmhe(i2d.1ea1)novwiaisdtioffeerleimntiinaatitoenthewiutnhknreoswpnecctoetofficpiesenutdεorf(rxe)- isthatwefindupdatesforq(x,s)viainfneriterationsexploring 199
(2.17)and(2.18)insideoftheinterval(0,1).Atthesametime, 200
184quencys.Differentiating(2.11)withrespecttosanddenoting
weupdatethetailfunctionV(x)viaouteriterationsexploring 201
185q(x,s)=∂sr(x,s),weobtain theentirereallineR.Inshort,giveonanapproximationforV(x),202
qxx+2s2qxrx+2srx2−2sqx−2rx=0, x∈E(0,1). (2.14) the algorithm updates q and then updated for εr(x). Next, the203
forward problem in (2.5) and (2.6) is solved for the function204
186Wenowneedtoexpressin(2.14)thefunctionrviathefunction w(x,s)fors=s.Next,thetailfunctionV(x)isupdatedusing205
187q.Wehave (2.16).Thismightseemoreminiscentofthestepsinalgorithms206
(cid:1)s suchasthemodifiedgradientinversescatteringtechnique[20];207
butweemphasizethat,unlikeourcase,suchmethodshaveno208
r(x,s)=− q(x,τ)dτ +V(x,s) (2.15)
I globalconvergenceproperties. 209
s
B. IterativerProcess 210
188whereV(x):=V(x,s)isreferredtoasthetailfunction,which
189is small in practice for large positive s. Here, the truncation We now outline the formulation of our algorithm and the 211
190pseudofrequency s serves as a regularization parameter. The iterative process (see details in [7], [14], [16], and [17]; see212
P
191exactformulaforV(x)is Fig. 1). Unlike computationally simulated data in [14], we213 AQ2
(cid:2) (cid:3)
do not use prior knowledge of the function q(1,s) on the214
1 w(x,s)
V(x,s):=V(x)=r(x,s)= ln . (2.16) transmittededgesincethisfunctionisunknowntous.Wehave 215
s2 w (x,s)
0 observedinourcomputationalexperimentsthattheknowledge 216
192Substituting(2.15)in(2.14),weobtainthefollowingnonlinear of q(1,s) only affects the accuracy of the calculation of the217
193integraldifferentialequation: locationofthetarget,butitdoesnotaffecttheaccuracyofthe218
q −2s2q (cid:1)sq (x,τ)dτ +2s(cid:1)sq (x,τ)dτ2 coonmlypiunttehdattacrognettr/absatc(ksegeroSuencdtiocnonIItrIa).sSt.inHceerεe,r(wxe)=are1ifnotrexres≥ted1222109
xx x x x and x0 <0, then one can easily derive from equations (2.5),221
s s (2.9),and(2.10)that∂xq(1,s)=0. 222
(cid:1)s Considerapartitionoftheinterval[s,s]intoN smallsubin-223
−2sqx+2 qx(x,τ)dτ tervalswiththesmallgridstepsizeh>0andassumethatthe224
s functionq(x,s)ispiecewiseconstantwithrespecttos.Thus 225
(cid:1)s
+2s2q V −4sV q (x,τ)dτ s=sN <sN−1 <···<s0 =s, si−1−si =h
x x x x q(x,s)=qn(x), fors∈(sn,sn−1]. (2.19)
s
+2s(Vx)2−2Vx =0, (2.17) Foreachsubinterval(sn,sn−1]weobtainadifferentialequation 226
x∈(0,1); s∈[s,s] forthefunctionqn(x).Weassign,forconvenienceofnotations, 227
q(0,s)=ψ0(s), qx(0,s)=ψ1(s) q0 :≡0.Followingtheaforementionedideaofacombinationof 228
qx(1,s)=0, s∈[s,s] (2.18) innerandouteriterations,weperformforeachninneriterations 229
4 IEEETRANSACTIONSONGEOSCIENCEANDREMOTESENSING
230withrespecttothetailfunction.Thisway,weobtainfunctions C. ComputingFunctionsqn,k(x)andV0(x) 262
231qn,k andVn,k.Theequationforthepair(qn,k,Vn,k)is Atfirstglance,itseemsthat,foragiventailfunctionVn,k(x),263
(cid:10)n−1 the function qn,k(x) can be computed as the solution of a 264
qn(cid:9)(cid:9),k−A1,nh qj(cid:9)−A1,nVn(cid:9),k−2A2,nqn(cid:9),k conventional boundary value problem for the ordinary differ-265
j=0    ential equation in (2.20) with any two out of three boundary 266
(cid:10)n−1 2 (cid:10)n−1 (cid:10)n−1 conditions in (2.21). However, attempts to do so led to poor267
=−A h2 q(cid:9) +2h q(cid:9)+2A V(cid:9) h q(cid:9) quality images (see [14, Remark 3.1]). At the same time, the268
2,n j j 2,n n,k j
QRMhasresultedinaccuratesolutionsbothin[14]andinTest269
(cid:13) j=(cid:14)0 j=0 j=0
−A V(cid:9) 2+2A V(cid:9) , x∈(0,1) (2.20) 1 for synthetic data (see succeeding discussion). The QRM is 270
2,n n,k (cid:9) 2,n n,k (cid:9) welldesignedtocomputeleastsquaressolutionsofPDEswith 271
q (0)=ψ , q (0)=ψ , q (1)=0 (2.21)
n,k s(cid:1)n0−,n1 n,k 1,n s(cid:1)nn−,k1 overdeterminedboundaryconditions,suchas,e.g.,theproblem 272
1 1 in (2.20) and (2.21). We refer to [18] for the originating work 273
ψ = ψ (s)ds, ψ = ψ (s)ds.
0,n h 0 1,n h 1 about the QRM and to [7], [9], [13], [15], and [16] for some274
sn sn follow-uppublications. 275
232Here,A1,n andA2,n arecertainnumbers,whoseexactexpres- LetL(qn,k)(x)andPn,k(x)beleft-andright-handsidesof 276
233sionsaregivenin[3]and[7]. (2.20),Erespectively. In our numerical studies, L(qn,k)(x) and 277
234 The choice of the first tail function V0(x) is described in Pn,k(x) are written in the form of finite differences. Let α∈278
235Section II-C. Let n≥1. Suppose that, for j =0,...n−1, (0,1)betheregularizationparameter.TheQRMminimizesthe 279
236functions qj(x) and Vj(x) are already constructed. We now followingTikhonovregularizationfunctional: 280
237need to construct functions qn,k and Vn,k for k =1,...,m. J (q )=(cid:12)L (q )−P (cid:12)2 +α(cid:12)q (cid:12)2 (2.25)
238WesetVn,1(x):=Vn−1(x).Next,usingthequasi-reversibilEity α n,k n,k n,k n,k L2(0,1) n,k H2(0,1)
239method (QRM) (see Section II-C), we approximately solve subject to boundary conditions in (2.21). Here, again norms 281
240(2.20) for k =1 with overdetermined boundary conditions in in L2(0,1) and in the Sobolev spacefH2(0,1) are understood 282
241(2.21)andfindthefunctionqn,1.Next,wefindtheapproxima- in the discrete sense. The functional Jα(qn,k) in (2.25) is283
242tionεr(n,1) fortheunknowncoefficientεr(x)viathefollowing quadratic.UsingthisfactandtheotoolofCarlemanestimates,it 284
243twoformulas: canbeshownthatJα(qn,k)hasauniqueglobalminimumand 285
E nolocalminima[14],[15],[17].Wefindthatglobalminimum286
(cid:10)n−1
r (x)= −hq −h q +V , x∈[0,1](2.22) viatheconjugategradientmethod,minimizingwithrespectto 287
n,1 n,1 j n,1
ε(n,1)(x)=1+r(cid:9)(cid:9) (x)+j=s20(cid:15)r(cid:9) (x)(cid:16)2 t1h0e0vgarliudespooifnttsheinfuthneoctiionnterqvna,lk(a0t, g1r)i.dTphoeinsttse.pWsiezehainvethuesesd-228898
r −2snr,1(cid:9) (x), n nx,1∈[0,1]. (2.23) direction was h=0.5. The s-interval was [s,s]=[3,12]. For 290
n n,1 each n=1,...,N, we take functions qn,k for k =1,...m,291
244Next, we solve the forward probIlem in (2.5) and (2.6) with and we typically choose m=10. The reason for the choice 292
224456εthri(sxw)a:=y.Aεr(fnte,1r)t(hxi)s,,wesu:p=dasteatnhdefitanildvtihaethfeufnocrtmiounlawinn,1(2(x.1,6s)), oeafcmho=ft1h0ernis,ttahialstsntaubmileirziecaaltkex≈per1i0e.nAceshtoasthsehroewgunlathriazta,tifoonr229934
247inwhichw(x,s):=wn,1(x,s).Thisway,weobtainanewtail parameter α, we have found, when testing synthetic data, that 295
248Vn,2(x).Similarly,wecontinueiteratingwithrespecttotailsm α=P0.04istheoptimalone,andwetakeitinourcomputations. 296
249times.Next,weset We note that we determined the regularization parameter 297
when testing simulated data. These data were for the target 298
qn(x):=qn,m(x), Vn(x):=Vn,m(x), εr(n)(x):=εr(n,m)(x) depicted in Fig. 7(a), for which we varied the regularization299
parameter between 0.03 and 0.05. The resulting images for300
250replacenwithn+1andrepeatthisprocess.Wecontinuethis these data showed only an insignificant change. We also var-301
251proces(cid:17)suntil[15] (cid:17) ied the regularization parameter between 0.03 and 0.05 for 302
either(cid:17)(cid:17)ε(n)−ε(n−1)(cid:17)(cid:17) ≤10−5 the experimental data. Again, we only observed insignificant 303
r r L2(0,1) changes, which lead us to select the average value of 0.04. 304
or (cid:12)∇Jα(qn,k)(cid:12)L2(0,1) ≥105 (2.24) Althoughtheregularizationtheorystatesthattheregularization 305
parametershoulddependonthenoiselevelinthedata[23],we306
252wherethefunctionalJα(qn,k)isdefinedinSectionII-C.Here, do not actually know the noise level for our data. Further, for307
253the norm in the space L2(0,1) is understood in the discrete nonlinear problems (as we have), this dependence is claimed308
254sense. In the case when the second inequality in (2.24) is byregularizationtheoryonlyforthelimitingcaseofarelatively 309
255satisfied,westopatthepreviousiteration,takingεr(n,k−1)(x)as smalllevelofnoise,whichisnotourcase.Inourcomputations 310
256oursolution.Ifneitheroftwoconditionsin(2.24)isnotreached usingmeasureddata,oneworkswithsomelevelofnoise,which311
257at n:=N, then we repeat the aforementioned sweep over the isnotlikelytobesmallandisunknown.Therefore,inpractice,312
258interval [s,s], taking the pair (qN(x),VN(x)) as the new pair when applying this algorithm to experimental data, we were 313
259(q0(x),V0(x)).Usually,atleastoneoftheconditionsin(2.24) guided by results from simulations to choose a value for the314
260is reached either on the third or on the fourth sweep, and the regularizationparameter.Ifwehadpriorknowledgeaboutsome315
AQ3 261processstopsthen. objectsinthetargetvolume,thenwewouldchoosetheoptimal 316
KUZHUGETetal.:IMAGERECOVERYFROMBLINDBACKSCATTEREDDATAUSINGANINVERSEMETHOD 5
317regularizationparameterforthatobject.Becauseweprocessed
318the data without any prior knowledge whatsoever about the
319objects, we chose the regularization parameter based on the
320simulateddataprocessing,andfortunately,ouranswersforfive
321outoffivetargetswerewellwithintabulatedlimits.
322 We now describe an important step in choosing the first
323tail function V0(x). To choose it, we consider the asymptotic
324behavior of the function V(x,s) in (2.16) with respect to the
325truncationpseudofrequencys→∞.Thisbehavioris[14],[17]
(cid:2) (cid:3)
p (x) 1
V(x,s)= 0 +O , s→∞.
s s2
326WetruncatethetermO(1/s2),whichissomewhatsimilarwith
327thedefiningofgeometricalopticsasahigh-frequencyapprox-
328imation of the solution of the Helmholtz equation. Hence, we
329take Fig. 2. (a) Schematic diagram of the forward-looking radar system illumi-
natingadEielectrictarget.(b)Typicalmeasuredtimehistoryofthebackscatter
p (x) field.(c)Compositeofunprocessedreturnshighlightingthedielectric target
V(x,s)≈ 0 . (indicatedbytheredtriangle).(d)Downrangecutofthepermittivityprofile,
s
whichthenewalgorithmwillgenerate.
330Sinceq =∂srandV(x,s)=r(x,s),then thenthefollowingestimateisvalid: 354
p (x) (cid:17) (cid:17)
q(x,s)=− 0s2 . (2.26E) (cid:17)(cid:17)εr(n,k)−ε∗r(cid:17)(cid:17) ≤(cid:18)hω (2.31)
L2(0,1)
331Hence,settingin(2.17)s:=sandusing(2.26),weobtainthe f
332followingapproximateequationforthefunctionp0(x): wherethenumberω ∈(0,1)isindependentofn,k,(cid:18)h,εr(n,k),355
d2 andε∗r. o 356
dx2p0(x)=0, x∈(0,1). (2.27) Therefore, Theorem 1 guarantees that, if the total number 357
E Q of computed functions εr(n,k) is fixed and error parameters 358
333Boundary conditions for p0(x) can be easily derived from σ,h are sufficiently small, then obtained iterative solutions359
334(2.18)and(2.26)as εr(n,k)(x)aresufficientlyclosetotheexactsolutionε∗r,;andthis360
p (0)=−s2ψ (s), p(cid:9)(0)=−s2ψ (s), p(cid:9)(1)=0. (2.28) closenessisdefinedbyotheerrorparameters.Therefore,thetotal 361
0 0 0 1 0
numberofiterationsQcanbeconsideredastheregularization 362
335We find an approximate solution p0,appr(x) of the problem in parameterofourprocess,whichistheadditionalregularization 363
336(2.27) and (2.28) via the QRM, simIilarly with the aforemen- parameter to the number s. The combination of inequalities 364
AQ4 337tionedequation.Next,wesetforthefirsttailfunction,i.e., in (2.30) and (2.31) has a direct analog in the inequality in365
r
p (x) [11, Lemma 6.2, p. 156] for classical Landweber iterations,366
V (x):= 0,appr . (2.29)
0 s which are defined for a substantially different ill-posed prob-367
lem.AstothetotalnumberofiterationsQbeingaregulariza- 368
338 A simplified formal statement of the global convergence P
tion parameter here, there is no surprise in this. Indeed, it is369
339theorem is as follows (see [7, Th. 6.1] for more details and
statedon[11,p.157]thatthenumberofiterationscanserveas 370
340[7,Th.6.7]forthe3-Dcase).
341 Theorem 1: Let the function ε∗(x) be the exact solution of aregularizationparameterforanill-posedproblem. 371
r
342ourCIPforthenoiselessdatag∗(t)in(2.3).Fixthetruncation
343pseudofrequency s>1. Let the first tail function V0(x) be III. IMAGINGRESULTS 372
344definedvia(2.27)–(2.29).Letσ ∈(0,1)betheleveloftheerror
Theschematicofthedatacollectionbytheforward-looking 373
345intheboundarydata,i.e.,
radar is shown in Fig. 2(a). Time-resolved electromagnetic 374
|ψ (s)−ψ∗(s)|≤σ, |ψ (s)−ψ∗(s)|≤σ, fors∈[s,s] pulses are emitted by two sources installed on the radar. Only375
0 0 1 1
one component of the electric field is both transmitted and376
346wherefunctionsψ0(s)andψ1(s)dependonthefunctiong(t)in measured in the backscatter direction. The data are collected377
347(2.3)via(2.7),(2.13)and(2.18);andfunctionsψ0∗(s)andψ1∗(s) by sixteen detectors with the step size in time of 0.133 ns.378
348depend on the noiseless data g∗(t) in the same way. Le√t h∈ Data from shallow targets placed both below and above the 379
349(0,1)bethegridstepsizeinthes-directionin(2.19);let α= groundwereprovided.Theonlypieceofinformationprovided 380
(cid:18)
350σ and h=max(σ,h). Let Q be the total number of functions bytheARLteam(SullivanandNguyen)toKuzhuget,Beilina,381
351εr(n,k) computed in the aforementioned algorithm. Then, there Klibanov,andFiddywaswhetherthetargetwaslocatedabove 382
352existsaconstantD =D(x0,d,s)>1suchthat,ifthenumbers the ground or was buried. The depth of the upper surface of a383
353σandharesosmall,that buried target was a few centimeters. GPS was used to provide 384
(cid:18) 1 thedistancebetweentheradarandapointontheground,which385
h< D2Q+2 (2.30) is located above that target to within a few centimeters error.386
6 IEEETRANSACTIONSONGEOSCIENCEANDREMOTESENSING
Fig.3. (a)Scatteredfieldfromametallictarget.(b)Scatteredfieldfromahigh-permittivitytargetwiththesameshape(εr(target)=10).Notethesimilarity
betweenthebackscatterelectricfieldsincases(a)and(b).
E
387Thetime-resolvedvoltagesinducedbythebackreflectedsignals illustrates this. We use the upper bound εr(target)=30 for425
388wereintegratedovertheradartotargetdistancesrangingfrom8 the metallic targets because our calculations show that LT in 426
389to20m,andtheywerealsoaveragedwithrespecttobothsource (2.7), from the response function g(t), almost coincides for427
390positions and with respect to the output of the 16 detectors. εr(target)≥30. 428
391Since we can assume here that the radar/target distance was In both cases of a metal structure and a high-permittivity 429
E
392known,thenitwasalsoapproximatelyknownwhichpartofthe structure, one can expect enhanced backscatter if the incident 430
393measured time-resolved signal would correspond to scattering pulseincludesfrequenciesthatcorrespfondtoanormalmodeof431
394events from that target (see Fig. 2). Despite the presence of thetarget.Hence,weassign 432
395clutter, a single time-dependent curve is extracted from the o
396measuredreturntimehistories,asillustratedinFig.2(b).This 10≤εr(metallictarget)≤30. (3.2)
397is the form of the data that have been procesEsed in each of
398the five measured data sets generated by the ARL. A typical We call (3.2) the appearing dielectric constant of metallic tar-433
399plotofreturnswithoutapplyingtheinversealgorithmisshown gets.Inotherwords,weconsiderin(3.2)thatregionsappearing434
400in Fig. 2(c), where the triangle denotes a possible target of tohaveahighdielectroicconstantcouldalsobemetallictargets.435
401interestamongtheclutterfromthebackscattergeneratedfrom Toappreciatethekindofbackscatterdataandimagerecov-436
402the volume of the region illuminated by the radar in Fig. 2(a). ery expected from a simple dielectric block, a 1-D example 437
403WeprocessasetofaveragedtimeIhistorieslikethoseshownin illustrated in Fig. 3 was investigated. Computations in this 438 AQ5
404Fig.2(b)tocreateadown-rangecutofthepermittivityprofile, example were performed using the software package WavES 439
405asindicatedinFig.2(d). [24]. The prermittivity profile, i.e., εr(target)=4, is shown in440
406 Ourobjectivewastocalculateratios Fig. 4(a); and the computed function u(0,t)=g(t) for 0<441
t<3 is shown in Fig. 4(b) [see (2.3) for g(t)]. We assume442
R= εr(target) (3.1) tempPoral units here for which at t=3, a distance of x=3443
εr(background) units is traversed; the source is at x0 =−1, and the block’s 444
frontfaceisatx=0.2.Sincetheblockis0.2unitswide,g(t) 445
407ofdielectricconstants.Iftheεr(background)isknown,thenit representsthebackscatterreturnfromthefrontandbackfaceof446
408istrivialtodeduceεr(target).Clearly,foratargetlocatedabove theblock.Thereasonwhy,inFig.4(b),g(t)=0fort<1and 447
409the ground, εr(background)=1. In general, we would expect g(t)=1/2for1≤t≤1.4isthatthesolutionoftheproblemin448
441110tthiaellytarvgaerytivnogluεmr(ex)to.Acownetaiginhtmedanavyeirnahgoemofogdeienleeicttireiscwcoitnhstsapnats- (w2h.1e)reanHd((z2).2is)tfhoerεHre(axv)is≡ide1fiusnuc0ti(oxn,,t)i.e=.,H(t−|x−x0|)/2, 444590
412of these constituent materials will be found over the volume (cid:19)
413spatial resolution cell that corresponds to the particular data 0,z <0
H(z)=
414acquisition configuration. In the examples presented here, we 1,z ≥0.
415showresultsobtainedfromjustonetime-historycurveforeach
416target,correspondingtoonepolarizationcomponentofthein- Hence, u(0,t)=g(t)=H(t−1)/2 for 1≤t≤1.4; and at451
417cidentelectromagneticfieldandbackscatterdatameasuredand t=1.4, the return wave from the block hits the observation 452
418averaged over all 16 receiver locations. Clearly, this severely point{x=0}forthefirsttime. 453
419limitsthetransverseresolutionbutimprovesthesignal-to-noise The measured data are also challenging to process since454
420ratio for 1-D imaging in the propagation direction. The model they arise from oblique illumination, and the exact location 455
421is further simplified by using the 1-D CIP employing only and the amplitudes of the incident pulses were not known. In456
422one hyperbolic PDE. Consequently, the interpretation of the addition, a comparison of Fig. 4(b) with Fig. 5(b), (d), and (f) 457
423backscattering radiation will assign a high-permittivity value showsthatthemeasureddataarehighlyoscillatory,whichare458
424to metal structures. A comparison between Fig. 3(a) and (b) unlike their simulated counterparts. Consequently, we applied459
KUZHUGETetal.:IMAGERECOVERYFROMBLINDBACKSCATTEREDDATAUSINGANINVERSEMETHOD 7
E
Fig. 4. (a) Function εr(target)=4; note that εr(background)=1. (b) u(0,t)=g(t) for 0<t<3.0. The source is located at x0=−1, and the first
backscatterreturnisthereforeshownatapproximatelyt=2.4with“ringing”determinedbyinterferenceofmultiplyscatteredwavesbetweenthetwoboundaries
oftheblock.ComputationswereperformedusingthesoftwarepackageWavES[24].
E
f
o
E
o
I
r
P
Fig.5. Threetargetsandtheirassociatedmeasureddata.Thegroundisdrysandwith3≤εr ≤5[21],[22].Theinformationshownin(a),(c),and(e)were
onlyprovidedaftercomputationsweremade.(a)Depictsabushthatwaslocatedonaroad,whichgeneratedbackgroundclutter.(b)Scaledexperimentaldatafor
(a),wherethehorizontalaxisrepresentstimeinnanosecondshavingatimestepof0.133ns;andtheverticalaxisistheamplitudeofthemeasuredvoltageatthe
detector.(c)Woodenstake.(d)Scaledexperimentaldatafor(c).(e)Metalboxburiedindrysand.(f)Scaledexperimentaldatafor(e).Themismatchbetween
experimentalandsimulateddata[seeFig.4(b)]isevident.