多角度遥感图像超分辨

AnOperationalSuperresolutionApproachforMulti-TemporalandMulti-Angle

RemotelySensedImagery

JianglinMa,JonathanCheung-WaiChan,andFrankCanters

Abstract—Inthispaperweproposeanoperationalsuperreso-lution(SR)approachformulti-temporalandmulti-angleremotesensingimagery.Themethodconsistsoftwostages:registrationandreconstruction.Intheregistrationstageahybridpatch-basedregistrationschemethatcanaccountforlocalgeometricdistor-tionandphotometricdisparityisproposed.Obstacleslikecloudsorcloudshadowsaredetectedaspartoftheregistrationprocess.ForthereconstructionstageaSRreconstructionmodelcomposedoftheL1normdatafidelityandtotalvariation(TV)regularizationisdefined,withitsreconstructionobjectfunctionbeingefficientlysolvedbythesteepestdescentmethod.OtherSRmethodscanbeeasilyincorporatedintheproposedframeworkaswell.Thepro-posedalgorithmsaretestedwithmulti-temporalandmulti-angleWorldView-2imagery.Experimentalresultsdemonstratetheef-fectivenessoftheproposedapproach.

IndexTerms—Deconvolution,registration,superresolution,WorldView-2.

I.INTRODUCTION

R

EMOTEsensingisplayinganincreasinglyimportantroleinmappingandmonitoringtheEarth.Increasingtheavailabilityofhighspatialresolutionremotesensingdataisanactivedriveformanyremotesensingapplicationssuchasurbanmapping,militarysurveillance,intelligencegathering,anddisastermonitoring.However,thereareseveralconstraintsthatmakehighspatialresolution(HR)remotesensingimagerywithhighspectraldefinitiondifficulttoobtain.Sensorsystemsoftenrequiremultipledetectorarrays,oneforeachspectralband;byloweringthespatialresolutionofeachdetector,morebandsmaybeputonthesensoratthesamecost[1].Hencethereisatradeoffbetweenspectralandspatialresolution.Inairborneremotesensing,securityconcernsmaypreventtheairplanefromflyingatlowheight,particularlyinurbanareas,thuslimitingthepossibilitiesforobtainingHRimagery.Otherconstraintsarerelatedtothemanufacturingtechnology.Whenthesensor’spixelsizeisreducedbeyondacertainlevel,shotnoisewillprevailanddegradetheimage[2].

ManuscriptreceivedJuly26,2011;revisedOctober21,2011;acceptedNovember14,2011.DateofpublicationJanuary18,2012;dateofcurrentversionFebruary29,2012.ThisworkwassupportedinpartbytheBelgianSciencePolicyOffice(BELSPO)undertheframeworkoftheSTEREOIIprogram—projectHABISTAT(contractSR/00/103).

TheauthorsarewiththeCartographyandGISResearchGroup(CGIS),De-partmentofGeography,VrijeUniversiteitBrussel,Brussels1050,Belgium(cor-respondingauthor,e-mail:Jianglin.Ma@vub.ac.be).

Colorversionsofoneormoreofthefiguresinthispaperareavailableonlineathttp://ieeexplore.ieee.org.

DigitalObjectIdentifier10.1109/JSTARS.2011.2182505

Aninterestingwaytoincreasespatialresolutionistocol-lectlowspatialresolution(LR)imagesofthesamesceneandfusethemintoaHRimageusinganimageprocessingtechniquecalledsuperresolution(SR).ThoughthedefinitionofSRvariesfromfieldtofield,multiframeSRisbroadlyunderstoodtobeatechniquethatcanrecoverhighspatialfrequencycomponentsfromLRimageryaliasedduetoundersampling[3].Amajorad-vantageofthisapproachisthatitcanmakefulluseoftheavail-ableLRimagingsystems,andmaythereforebeaneconomicallyinterestingalternativeforimageresolutionimprovement.

SRimageenhancementapproacheshavereceivedalotofinterestinthepasttwodecades,andavarietyofSRalgorithmshavebeenproposedintheliterature.TheideaofSRwasfirstintroducedin1984byTsaiandHuangforLANDSAT4images[4].Kimetal.generalizedthisworktonoisyandblurredimages,usingleastsquareminimization[5].Inbothofthesestudies,SRisperformedinthefrequencydomain.Althoughfrequency-basedSRmethodshavetheadvantageofsimplicityandlowcomputationalcomplexity,theyrequirethatthegeometricdisparitybetweenLRimagesistranslational.Inordertoovercomethisproblem,variousspatial-basedSRapproacheshavebeenproposedastheyprovidemoreflexi-bilityinmodelingtheimagedegradationprocess.ASRspatialdomainalgorithmwasfirstpresentedbyUrandGrossin1992[6].Utilizingthegeneralizedmultichannelsamplingtheoremtheyproposedanon-uniforminterpolationmethodformultiplespatiallyshiftedLRimages.Theinterpolationisfollowedbyadeblurringprocess.AdifferentmethodcalledIterativeBack-wardProjection(IBP),whichwasadoptedasabasicalgorithmincomputer-aidedtomography,wasfirstsuggestedbyIraniandPeleg[7].Iftheimagedegradationelementsincludingblurring,warping,sampling,andadditivenoiseareknownexactly,SRbecomesaninverseproblemsimilartoimagerestoration,whichhasbeenstudiedforamuchlongertimethanSR.ManyalgorithmscomingfromtheimagerestorationdomaincanthenbeusedforSR,andinsomewaySRcanbethoughtofassecond-generationimagerestoration.Toreducenoiseandsolvesingularcases,severalSRalgorithmsincorporatepriorknowl-edgeintothecomputationbyconstrainingthesolution.Forexample,Lietal.proposedthemaximumaposterioribasedmethod(MAP)usingadiscontinuitypreservinguniversalHiddenMarkovTree(HMT)modelforSRreconstruction[8].Foranoverviewofstate-of-the-artSRmethods,wereferto[2].NotallLRimages,however,canbeemployedforSRimageenhancement.First,imagesmustbeundersampled,thatistosay,aliasesshouldexistinLRimageryofthesamescene.Second,noneoftheLRimagesmaybeobtainedfromotherLRimages

1939-1404/$31.00©2012IEEE

MAetal.:ANOPERATIONALSUPERRESOLUTIONAPPROACHFORMULTI-TEMPORALANDMULTI-ANGLEREMOTELYSENSEDIMAGERY111

intheset.Mostdigitalimagessatisfythefirstrequirementastheopticsusuallyresolvesmuchfinerdetailthanwhatcanbecapturedbythesensorarray.Thesecondrequirementisusu-allyfulfilledthroughthegeometricdisparitybetweentheLRimages.

IncontrastwiththevastamountofSRapplicationsindo-mainslikevideosequencingandmedicalimaging,successfulSRapplicationsinremotesensingarecomparativelyrare[9].ThemaindifficultyinremotesensingliesinthefactthatitisnoteasytoobtainasetofLRimagesthatcanbeusedforSR.Ingeneral,LRremotesensingimagessuitedforSRcanbecol-lectedinthreeways:

1)Shift-controlledsituation:inthiscasemultipleimagesareacquiredsimultaneouslywithknownpixelshiftbetweeneachacquisition.AtypicalexampleistheSPOT5imagingsystem,whichcapturestwoimageswithahalf-pixelshiftfromadoubleCCDlineararray.UsingSRtechniques,SPOT5providespanchromaticimageryat2.5mspatialresolutioninsteadoftheoriginal5mresolution[10].2)Multi-temporalsituation:manycurrentsensorshaveaveryhightemporalresolution.Forexample,theModerate-Res-olutionImagingSpectroradiometer(MODIS)carriedontheTERRAsatellitecanacquireimageryofthesamesceneeveryoneortwodays.Therefore,itcanprovideasequenceofLRimagesaccuratelyandquickly.

3)Multi-anglesituation:todaymanysatellite-basedimagingsystemsareequippedwithmulti-anglecapabilities,in-cludingMultispectralThermalImager(MTI),IKONOS,Quickbird,WorldView-2,Multi-angleImagingSpectro-Radiometer(MISR),AlongTrackScanningRadiometers(ATSR-1,ATSR-2,AATSR),andCompactHighRes-olutionImagingSpectrometeronboardtheProjectforOn-boardAutonomy(CHRIS/Proba).

Themulti-anglesituationcanberegardedasaspecialcaseofthemulti-temporalsituationasimagescapturedfromdif-ferentanglesareusuallynotavailablesimultaneously.Inthispaperwewillfocusonthemulti-temporalandmulti-anglesitu-ationasthisismostcommoninremotesensing.AlthoughmanystudiesonSRhavebeencarriedoutusingmulti-temporalandmulti-angleremotesensingimagery,mostSRworkisbasedonsimulateddatainsteadofrealremotesensingdata.AnobviousadvantageofusingsimulateddataisthatSRimprovementcanbeeasilyquantifiednumerically.AssuchitisnowonderthatmostproposedSRalgorithmsaretestedonsimulateddataasthecomparisonofnewapproacheswithconventionalSRalgo-rithmsbecomesmucheasier[1],[2],[11],[14].

Whilesomestudiesarebasedonrealdatatheyusuallydealwithsmallimagepatches,ratherthanwithlargerimageex-tracts.Forinstance,in[9],theSRmethodwastestedonLandsatimageextractsof256256pixels,whilein[8],MAP-basedSRwithaHuberpriorwastestedonmulti-temporalMODISdatausingimageextractsof50by50pixels.Whendealingwithsmallimagepatches,thegeometricdisparitybetweenLRimagescanbeeffectivelyexpressedasarigidtransformmodel(translational,affineorprojective).Moststate-of-the-artSRmethodsmakeuseofrigidtransformmodelsbecauseitenablessimulationofLRimagesfromapreliminaryreconstructedHRimage.StudieswhereSRisappliedtolargerimageextractsusuallyrelyonnon-uniformbasedSR,whichdecomposesSR

intothreeindependentprocedures:registration,interpolationanddeconvolution.Inthiscase,intheregistrationphase,anon-rigidtransformmodelisusedtoaccountforthecomplexgeometricdisparitybetweentheLRimages.AfterLRimagepixelpositionsareidentifiedintheHRgrid,non-uniformin-terpolationisperformed.Forexample,MerinoandNunez[12]andChanetal.[13]appliedVariable-PixelLinearReconstruc-tion(VPLR)andDelaunaytriangulation-basednon-uniforminterpolationmethodsrespectivelyonLANDSATETM+andCHRIS/Probaimagery.However,implementingHRimagereconstructionasasequenceofthreeindependentproceduresissub-optimalastherestorationdoesnottakeintoaccounttheerrorsthatoccurintheinterpolationstep[2].Nexttotheneedtodevelopapproachesthatenableapplicationofstate-of-the-artSRmethodsonfull-sizeimagery,workingwithrealimagedataalsorequireseffectivesolutionstodealwithobstaclesthatmayseverelydegradeSRresultssuchascloudsorcloudshadows[9].

Inthispaperweproposeanoperationally-orientedSRap-proachformulti-temporalandmulti-angleremotesensingim-agery,wheretheissuesoflocaldistortionandpresenceofcor-ruptedpixelsaretackledbyapatch-basedSRschemeandageneralizedSRreconstructionmodelrespectively.Thebasicideabehindthepatch-basedSRschemeisthatSRisperformedonthebasisofsmallimagepatches,forwhichthetransformmodelcanbeexpressedasaffine.Usingapatch-basedapproach,thetransformmodelislocallyrigidbutgloballynon-rigid.ThegeneralizedSRreconstructionmodelisbasedonminimizingaHRobjectfunctionincludinganL1normdatafidelitytermandatotalvariation(TV)regularizationterm.Byintroducinganobstaclemaskduringtheimageregistrationprocedure,andtheL1normdatafidelitytermduringtheimagereconstruc-tionprocedure,onlynon-corruptedLRpixelsareusedforesti-matingtheregistrationtransformmodelandforreconstructingtheimage.Thereconstructionobjectfunctioncorrespondingtothemodelcanbeefficientlysolvedbythesteepestdescentmethod.TheoperationalSRapproachdescribedinthispaperalsoprovidesaframeworkforapplyingotherstate-of-the-artSRmethodstomulti-temporalandmulti-angleremotesensingim-agery.Indeed,aswillbeshown,otherSRmethodscanbeeasilyincorporatedintheapproachthatisproposed.

Theorganizationofthispaperisasfollows.SectionIIde-scribestheproposedmethodindetail,whiletheapplicationoftheproposedmethodonmulti-angleWorldView-2imageryisreportedinSectionIII.SectionIVconcludesthepaperwithrec-ommendationsforfuturework.

II.SUPERRESOLUTIONMETHODOLOGY

A.ObservationModel

ReconstructingaHRimagefromLRdatarequiresthedefini-tionofanobservationmodeldescribingtheimagingprocess.ThemostpopularmodelwasintroducedbyEladandFeuer[14],wherethepotentialHRimageundergoeswarping,blur-ringanddownsamplingbeforeLRimagesaregenerated.Capel[15],GevrekciandGunturk[16]generalizedthismodelbyin-cludingphotometricmodeling,asexternalilluminationcondi-tionsmayvaryfordifferentinputLRimages.AlthoughCapel’smodelhasbeenusedinremotesensing[9],aparticulardifficulty

112IEEEJOURNALOFSELECTEDTOPICSINAPPLIEDEARTHOBSERVATIONSANDREMOTESENSING,VOL.5,NO.1,FEBRUARY2012

Fig.1.ObservationmodelthatrelatestheLRimagestothepotentialHRimage.

intheapplicationofthismodelliesinthefactthatobstaclessuchascloudsandcloudshadowsmayoccurinremotesensingim-agery,leadingtocorruptedpixels.Also,forthemulti-temporalandmulti-angleSRcase,changesmayoccurovertime,leadingtoinconsistenciesinSRreconstruction.Wethereforeproposeanewimageobservationmodel:

(1)

whereistheimagesetcomposedofLRimageswitheachLRimageofsizebeingrepresentedinlexicographicalorder

anddenotes

theobstaclemaskandthenoiseforLRimagesrespectively;representstheHRimagetobeobtainedinlexicographicalorder

,withtheHRimagesize

andthemagnificationfactor,usuallylargerthan1;istheelementbyelementmultiplicationoperator,constraining

therelationshipbetween

andtothenon-obstacleregion,andandareoperatormatricesthatrelateto:

(2)

where

andareofthesamesize,standingforthewarpmatrixandblurmatrixrespectively,and

isthedecimationmatrixofsize.Thewarp

matrix

representsthegeometricrelationshipbetweenLRimagesandtheunderlyingHRimage,whichisusuallyobtainedviaimageregistration.Theblurmatrixiscloselyrelatedtothepointspreadfunction(PSF)oritsequivalentinthefrequencydomain,theopticaltransferfunction(OTF).UsuallythePSFisassumedtobespatiallyinvariantandismodeledbylow-passkernelssuchasGaussianfunctions.Thedecimation

matrix

isdeterminedbyboththePSFandthemagnificationfactorordownsamplingrate.Itisdownsamplingthatdistin-guishesSRfromimagerestoration.andarethecontrastandbrightnessmatrixrespectively:

(3)

(4)

with

istheunitvectorofsize.

andareaffineandcanbeestimatedbyphotometricreg-istration.Suchasimpleaffinemodelisoftenfoundusefulforimagescapturedatdifferenttimes[9]orfromdifferentangles[15].istheobstaclemaskmatrixidentifyingtheLRimagepixelsthatcanbeusedforSR.canbeobtainedeitherbyreg-istration(theoutliersintheregistrationprocedure)orbychangedetectionafterregistrationhasbeendone.isthenoisematrix,whichisintroducedtoaccountforelectronicnoise,shotnoiseofphotondetectorsaswellaserrorsintroducedbyinaccurategeo-metricregistration.Dependingontheassumptionaboutnoise,thenoisemodelcanbeeitherGaussianorPoisson.

TheblockdiagramshowninFig.1depictsthesystemtobemodeled,andlike(1),describestherelationshipbetweentheobservedLRimagesandthepotentialHRimage:aHR

imagegoesthroughwarping

,blurring,down-sam-pling

,photometricvariation,andiscontami-natedbyobjectobstaclesandnoise,resultingintheLRimagesthatareobserved.AccuratelyestimatingeverydegradationcomponentintheobservationmodelisapremiseforSR.Ifthedifferentdegradationcomponentsareknown,ob-tainingtheHRimagebecomesanoptimizationproblem,whichcanbesolvedbyavarietyofoptimizationalgorithms.Basedontheproposedobservationmodel,thecorrespondingSRal-gorithmprocedureisstraightforwardasFig.2illustrates.Eachstepintheprocedurewillbedescribedinmoredetailbelow.B.RegistrationandObstacleDetection

AccurategeometricregistrationiscriticalforsuccessfulSR.Multi-temporalandmulti-angleimagesposeamajordifficultyforSR:localgeometricdistortioncausedbytopographiceffectsand/orplatforminstabilitymaybeimportant,whichmakesrigidtransformationmodelssuchasaffineorprojectivetransformsunreliable.Recently,weproposedtheuseofanon-rigidthinplatespline(TPS)transformmodelformulti-angleSR[17]anddemonstratedthatanon-uniformSRapproachbasedonanaccuratelyestimatedTPSmodelcanperformbetterthanotherconventionalSRalgorithms[13].However,astheTPSmodelisnotinvertible,itrestrictsSRalgorithmsthatcanbeemployedincombinationwiththisregistrationframeworktothenon-uni-forminterpolationapproach.Inordertoobtainhigh-qualitySR

MAetal.:ANOPERATIONALSUPERRESOLUTIONAPPROACHFORMULTI-TEMPORALANDMULTI-ANGLEREMOTELYSENSEDIMAGERY113

Fig.2.FlowchartshowingtheSRprocessingchain:inthemethodproposed

steps

areperformedsimultaneously.results,however,moststate-of-the-artSRalgorithmsrequire

theregistrationtransformmodeltobeinverseconsistent,meaningthattheforwardandreversemappingmatricesshouldbemutuallyinversetoeachother.Aseeminglyobvioussolutionwouldbedensetransformoropticalflow,whereeverypixelisdefinedbyonedisplacementvector.However,densetransformmayhavetoomanydegreesoffreedom,resultinginlocalirregularitiesbecauseofthepresenceofcorruptedpixels.TheresultingregistrationerrorsthereforemayrenderSRinfeasible[18].Moreover,itscomputationalburdenisnotablyheavy.Inourproposedpatch-basedregistrationscheme,insteadofassigningeachpixelonedisplacementvectorweassigneachimagepatchaparameterizedlocallyaffinetransformationmodel.Whileestimatingtheaffinetransformationforanimagepatchpair,thephotometrictransformandtheobstaclemaskcanbeestimatedsimultaneously.Thepatch-basedregistrationschemeisshowninFig.3andiscomposedoffoursteps:

1)Firstlyathird-orderpolynomialtransformmodelissetuptodefinetheoverallgeometricrelationshipbetweentheinputimage(imagetoberegistered)andthereferenceimage.

2)Thenthereferenceimageisdividedintooverlappingpatches.Foreachpatchitsgeometricallycorrespondingimagepatchintheinputimagecanbeeasilyfoundviatheabove-estimatedthird-orderpolynomialtransformmodel.Foreachpatchpair,anaffinegeometricandphotometricmodelexplicitlyaccountingforobstaclesintheimageryisestimated.

3)Eachpatchintheinputimageisresampledwiththeesti-matedaffinetransformmodel.

4)Inordertoobtainspatiallyconsistentresults,overlappingpartsofpatchesarefusedbysimpleaveraging,andafullimageisobtained.

HereweshouldmentionthatintheSRprocedure(Fig.2)theoutputoftheregistrationisthelocalaffinetransformmodelratherthantheregisteredimage.Therefore,intheSRprocedurestep(3)intheregistrationprocessissubstitutedwithHRrecon-

Fig.3.Flowchartshowingtheregistrationprocessingchain.

structionandisfollowedbyfusingandmosaicingofHRimagepatches,aswillbeexplainedinmoredetaillater.

Thebasicideabehindthepatch-basedregistrationschemeinourapproachistorefinetheimageregistrationprocessbyusingahierarchicaltransformationschemethatcombinesdif-ferentmethods.Asimilarideaisalsofoundin[17]and[19].However,theschemedefinitionandthemethodsapplieddifferfromtheapproachesdevelopedinthesepapers.Thepatch-basedSRreconstructionschemeismotivatednotonlybytheregistra-tionbutalsobytheimagerestorationprocess.Manyimagede-convolutionalgorithmsassumestationarityofsignalandnoiseprocessesandspaceinvarianceoftheblur.Whiletheseassump-tionsrarelyholdforanentireimage,theycanbeassumedtoholdforalocalimagepatch.Therefore,itseemsnaturaltodividethelargerimageintooverlappingimagepatchesandprocesseachpatchbasedonitsowncharacteristics,whichisexactlythephilosophybehindourpatch-basedSRreconstructionscheme.Inthecontextofthespatially-varyingPSFimagerestorationproblem,thepatch-basedschemeisalsoreferredtoasthesec-tioningmethod[20].

1)GlobalThird-OrderPolynomialModelEstimation:Inordertoestimateanaccuratethird-orderpolynomialmodelplentyofhigh-qualitycontrolpoints(CPs)arerequired.Inearlierworkonmulti-angleimagerestoration[17],afterana-lyzingtheadvantagesanddisadvantagesofregion-basedandfeature-basedmethods,weadoptedthescale-invariantfeaturetransform(SIFT)localfeaturedescriptortofindcandidateCPs,asSIFT-baseddescriptorsarepreferablewhenchangesoccurinscale,illumination,andtoacertainextentinthe3Dcameraposition.ToensurethequalityofCPs,inourapproacharig-orousscreeningprocedureisappliedincluding:1)anambiguitytestbasedonSIFTfeatures,whicheliminatesambiguousCPs;2)theapplicationofaspatialconstraintbasedonthe-esti-matorsampleconsensusmethod,whichdiscardsCPsthatdonotconformtoaunanimousprojectivetransformmodel;and3)aniterativemethodforCPsrefiningusingthestatisticalcharacteristicsofthedistributionofthird-levelpolynomialtransformrandomerrors.

114IEEEJOURNALOFSELECTEDTOPICSINAPPLIEDEARTHOBSERVATIONSANDREMOTESENSING,VOL.5,NO.1,FEBRUARY2012

Theimagepatchdivisionschemewillbepresentedlater.Inthenextsection,weassumetheimagepatchpairsareavailable.2)GeometricandPhotometricRegistration,andObstacleDetectionforLocalImagePatchPairs:Althoughitisimpos-sibletorepresentthegeometricdisparitybetweenmulti-tem-poralandmulti-angleimageswithanaffinetransform,forsmallimagepatchesaffinetransformationissufficient.Toreachsub-pixelaccuracywithaffinetransformmodelestimation,therearetwogroupsofmethods,asfollows:

1)Oneapproachisbasedoninterpolation;inthiscaseithasbeenshownthatthemainfactorsaffectingregistrationac-curacyaretheinterpolationfunction,samplingfrequency,numberofbitsperpixel,andfrequencycontentoftheimage[21].

2)Theotherapproachisbasedonthedifferentialpropertyoftheimagesequenceundertheintensityconservationas-sumption;forthisapproach,whichhasbeenfoundmostpopularintheliterature,thereisnoneedforinterpolation[7],[9],[13].

Theintensityconservationassumption,however,doesnotholdformulti-temporalandmulti-angleremotesensingimagepatchpairsasluminancevariationsareoftenpresentbetweenmulti-angleremotesensingimages.Eventhoughanaffinepho-tometricmodelcanbeappliedtocompensateforintensitydis-parity,corruptedpixelsdueto1)saltandpeppernoise;2)pres-enceofcloudsorcloudshadowsand3)surfacetypeswhosebidirectionalreflectancedistributionfunctionchangesdramati-callyoftenfailtheaffinephotometricmodel.

Inordertosolvetheabove-mentionedproblemsweadoptarobustdifferentialregistrationapproachfirstproposedformed-icalimageprocessing[22].Wefoundthisapproachespeciallysuitedforourcasethoughourpatchdivisionschemeandthesmoothnessconstraintappliedforobtainingagloballyconsis-tentresultisdefineddifferently.Tocompensateforintensityvariations,themethodincorporatesanexplicitphotometricchange(contrastandbrightness)intothegeometricmodel,performinggeometricandphotometricregistrationsimultane-ously.Toaccountforoutliereffects,themodelalsoincludesanobstaclemask,separatinguncorruptedimagepixelsfromcorruptedones.

Supposetheinputimagepatch

canbewellregisteredtothereferenceimagepatch,andthenon-obstacleregionand

obstacleregionaredenotedas

andrespectively.Fortheimageregionthegeometricandphotometricrelationship

between

andcanbeexpressedas(5)

wherearetheaffinegeometric

parameterswhilearetheaffinephotometricparam-eters.Inordertogetanaccurateestimationofparameters

onlyuncorruptedpixelsbelongingto

theimageregionshouldbeemployed.

Nowweassumethatthepixels(eitherin

orin)fromthedifferenceimagepatchobtainedbysubtractingthe

registeredpatchfromthereferencepatcharespatiallyindepen-dentandidenticallydistributed.Forpixelsin

weassumeaGaussiandistributionwithzeromeanandstandarddeviation,whichisacommonassumptioninregistration[15],[17]:

(6)

where

isthegrayvalueofthepixellocatedatinthedifferenceimagepatch:

(7)

Forpixelsin

weassumeauniformdistribution:

(8)

Thedifferenceimagepatchisanindicatorofregistrationquality.Supposenowtheinputimagepatchiswellregisteredtothereferenceimagepatchviatheregistrationmodel.Thenthelikelihoodofthedifferenceimagepatchshouldreachthemaximum:

(9)

In(9)weassumethatthepossibilitythatapixelbelongstoandisequal,andandareconstantsthatcanbede-terminedbyand.Accordingtothemaximumaposteriori

criterion,when

reachesthemaximumitsfirstorderderivativewithregardtoshouldbeequaltozero:

(10)

wheredeter-mineshowapixelaffectstheestimationofregistrationparame-ters.Ifthepixelhappenstobecorrupted

tendstobelarge,leadingtoasmallweightingfactor,henceitsinfluenceonregistrationparameterestimationisplayeddown.

canbefurthersimplifiedwiththefirstorderTaylor

seriesexpansionof:

(11)

whereandarethe-spatial,-spatialand

temporalderivativeof

respectively.Bysubstituting(11)in(10)wegetaquadraticerrorfunction,whichisnowlinear

MAetal.:ANOPERATIONALSUPERRESOLUTIONAPPROACHFORMULTI-TEMPORALANDMULTI-ANGLEREMOTELYSENSEDIMAGERY115

foritsunknown.Therefore,weobtainananalyticexpression

for:

(12)

where

and

.From(12)itseemsasifwe

canobtaintheregistrationparametersdirectlywithoutitera-tion.However,theestimatorforrequiresanestimatorfor,whichalsorequiresanestimatorfor.Therefore,theexpec-tation-maximization(EM)algorithmisusedforsolving,andthewholealgorithmprocedureisdefinedasfollows:

1)Initialization:assignthefixedconstantsand,andpro-videinitialestimatefor.

2)E-Step:compute.Thisstepalsoinvolvestheprocessofclassifyingpixelseitherintotheregionorintotheregion.Pixelsbelongingtohaveaclosetozeroweight.3)M-Step:calculate.Thiscanbedonewiththeclassicaldifferentialmethodproposedin[7],whichreliesonaNewton-Raphsonstyleiterativeschemetoapproximatethesolution,andaGaussianpyramidschemetoincreasetheconvergencespeed.

4)Thelasttwostepsarerepeateduntilthedifferencebetween

successiveestimatesof

fallsbelowasmallthresholdvalue.

3)Patch-BasedRegistrationandSuperresolution:Inourex-perimentstheimagepatchsizewassetat9696asacompro-misebetweenaccuratemodeladaptationandestimation.Onceimagepatchesareaccuratelyregistered,itisimportanttoelimi-natetheblockingeffecttypicalforpatch-basedapproachesandmakesurethatadjacentimagepatchesmatch.Therearetwowaystoeliminatetheblockingeffect[23]:oneistoapplyanav-eragingkernelwithanappropriatesizetosmooththeresultingimage,whichwillinevitablydecreasetheimagequality;theotheristodividetheimageintooverlappingblocksasshowninFig.4,whichisthemethodadoptedinthispaper.InFig.4thesmalltileisofsize48by48,andeachimagepatchiscom-posedoffourtiles.Usingthemovingwindowapproach,intheexamplethereare49imagepatchesalltogetherinsteadof16,withthepatchcenterindexedasillustratedinFig.4.Foranyimagetile(exceptthoseneartheedge)spatialconsistencyisobtainedbyaveragingthedifferences(geometricdisplacement,photometricdisparityandreconstructeddiscrepancy)forthefourfloatingpatchesthetilebelongsto.Byadoptingthisslidinghalf-windowmethod,wecanmakesurethatthetransformationbetweenimagesislocallyaffinebutgloballysmooth.

Thepatch-basedschemehasthefollowingmeritsforSR:1)Althoughtheglobaltransformisnotinvertible,foreachpatchthetransformisinvertible,thereforestate-of-the-artSRalgorithmscanbeimplementedinthisscheme.

2)Remotesensingimagesoftenhavealargesize;therefore

theconstructionofthe

matrixin(2)mayposecompu-tationalproblems.Constructingthematrixatpatchlevelisefficientandallowsapplicationofthemethodonfull-Fig.4.Anillustrationofoverlappingimagepatches.Thesmalltileisofsize4848,andtheimagepatchiscomposedof4tiles.Thesolidlinerepresentsthepositionsofnon-overlappingimagepatcheswhilethedashedlinesdenotetheslidingpositions.Thecentersofoverlappingimagepatchesareindexedfrom1and49.Onlytheimagetilesshowningreencanbeprocessedusingtheapproachproposed.

sizeimagery.Thepatch-basedschemealsoenablespar-allelprocessingofSR.

3)Althoughpatch-basedimageregistrationhasbeencriti-cizedbecauseimagedetailsmaybeoverlysmoothed,andapolyaffineschemehasbeenproposedtosolvethisproblem[24],wefoundthatsimplefusionofSRpatchesbyaver-agingdoesnotdeterioratetheSRresult.Inthespatially-varyingPSFimagedeconvolutionliterature,thissimplefu-sionmethodiscommonpractice[20].C.PSFEstimationandConstructionoftheHMatrixFrom(2)wecanseethatthePSFthatdefinestheblurmatrixisappliedontheunderlyingHRimageinsteadofontheLRimage.Therearetwowaystoderiveit,asfollows:

1)Oneisreferredtoastheblind-SRmethod,whichassumesthattheHRimagecanbemodeledbymeansofaninho-mogeneousimagemodelsuchasthePotts-Markovmodel[25].However,thecomplexityofaremotesensingsceneoftendoesnotallowtheuseofthisapproach.Moreover,thesimultaneousestimationofthePSFandtheHRimageusuallyrequiresiterationand,thereforeiscomputationallyheavy.

2)TheotherpossibilityistofirstestimatethePSFoftheLRimageandthenlinkingittothePSFoftheunderlyingHRimage.ComparedtoestimatingthePSFblindly,on-orbitdeterminationofthePSFofLRimages,wherethePSFiscalculatedbywayofspecificobjectsintheobservedimages,ismorereliable.ThelocalaffineassumptionalsomakesiteasiertodefinetherelationshipbetweenthePSFoftheLRimageandthatoftheHRimage.Therefore,inthisstudythisstrategyisadopted.

116IEEEJOURNALOFSELECTEDTOPICSINAPPLIEDEARTHOBSERVATIONSANDREMOTESENSING,VOL.5,NO.1,FEBRUARY2012

1)PSFEstimationofLRImages:IntheSRliteraturethesystemPSFisusuallyrepresentedbya2D-Gaussianshapefunc-tion,whichcanbedeterminedbyitscovariancematrix[7],[9],[13],[14].Ofallavailableon-orbitPSFestimationmethodsweadoptedtheknife-edgeinputmethod,asthismethodhasbeenwidelyusedforIKONOS,CBERS-2,SPOT,OrbView-3andBiLSATPSFestimation.Anewimplementationoftheknife-edgeinputmethodhasbeenproposedrecentlyanditseffec-tivenesshasbeendemonstratedwithdifferentsimulatedim-agery[26].Forapplicationonrealdata,arobustoptimizationschemewasadoptedtoreducetheinfluenceofbadmeasure-ments.Bywayofthismethod,ananisotropic2-DGaussianPSFisobtained:

(13)

whereisthealong-trackdirection,andandarethe-co-ordinateand-coordinatestandarddeviationmeasuredintermsofLRimagepixels.

2)PSFEstimationoftheUnderlyingHRImage:WhenthePSFisprojectedfromtheLRimagetotheHRimage,itsshapeundergoesthesamewarpingasspecifiedinthewarpmatrix.Sinceinourpatch-basedSRschemethelocaldefor-mationisassumedaffine,itisreasonabletoassumethattheanisotropicGaussianPSFbecomesamoregeneralanisotropic2-DGaussiankernel,whosecovariancecanberepresentedby

(14)

where

isthecovariancematrixofthePSFoftheLRimage,whichisdeterminedbyandin(13):

(15)

andistheJacobianoftheaffinetransformmatrixwhichrelatestheLRimagewiththepotentialHRimagegeometrically:

(16)

withandthe-coordinateand-coordinateup-sampling

scalerespectively,and

thesameaffinegeo-metrictransformmodelparametersasin(5).

3)Constructionofthe

Matrix:InmostSRpaperscon-structionofthematrixanditsadjointoperatorisunnec-essaryasitcanbeimplementedbydirectimageprocessingop-eratorssuchasimageconvolution(forand),downsam-pling,andupsamplingwithzeropadding.However,thiskindofimplementationstrategyisconsistentwiththeob-servationmodelonlyinthecaseoftranslationaldisplacement(undoubtedlyimproperforouraffinecase)andintegerup-sam-pling(increasingpixelnumbersineachdirectionbyafactor2,3,4,)[25].Theintegerup-samplingschemeisnotonlyinconvenientinpracticebutalsoinefficientfromatheoreticalpointofview.Indeed,ithasbeenfoundthatbetterSRperfor-manceisobtainedbyemployingnon-integermagnificationfac-torsratherthanintegerones[27].

Inordertoovercometheaboveimplementationconstraints,

weadoptedthestrategyofdirectlyconstructingthe

matrix,whichismadepossibleduetotheuseofapatch-basedSRscheme.DirectconstructionofthematrixalsomakesiteasiertotakeadvantageofavastrangeofoptimizationalgorithmsavailableintheliteratureforthefinalHRimagereconstruction.Fromtheobservationmodel,wecanseethateachLRimagepixeliscreatedbymultiplyingablurkernel(anisotropic2-DGaussian)withitsunderlyingHRimagepixelsandtakingthesumofthepixelintensityvalues.ThecenteroftheblurkernelisdeterminedbythelocationoftheLRimagepixelwhenitisprojectedontotheHRgrid.NowsupposeisthepixelindexinlexicographicalorderforthethLRimage

denotesitspositionintheLRgrid,while

representsitspositionintheHRgrid.Ifisdefinedas

thepixelindexinlexicographicalorderfortheHRimage,and

standsforitspositionintheHRimagegrid,then

theelementofcanberepresentedas[28]

(17)

(18)

whereisequaltoin(14).Sincetheblurkernelisspatiallyconstraineditisreasonabletoassumethatpixelsfarfromthekernelcenterhavezeroweight,henceleadingtoaverysparsematrix.D.HRImageReconstruction

Aftergeometricregistration,photometricregistrationand

on-orbitPSFestimation,the

matrixcanbereconstructed,whichisapremiseforHRimagereconstruction.However,isofteneitherunder-determinedorill-conditioned,whichmeansthatthereisaninfinitenumberofHRimagesthatconformtotheobservationmodelorthatsmallamountsofnoiseintheLRimageswillresultinlargeperturbationsintheunderlyingHRimage.Hence,HRreconstructionbecomesahighlyill-conditionedinverseproblem.

The‘goldenrule’forsolvingthisproblemistosearchforanapproximatesolutionsatisfyingadditionalconstraintsrelatedtothephysicsoftheproblem[29].Theadditionalconstraintsareoftenreferredtoasaprioriorpriorinformationastheyarenotderivedfromtheobservationmodel.Arealisticassump-tionaboutthepotentialHRimageisthattheimageconsistsofsmoothregions,separatedbycrispedges.Thisisthebasicideaofedge-preservationregularizationmethodssuchastheTVmethod.TVregularizationhasbeenappliedinSRHRimagere-construction[30],[31],[32],however,thewayhowwedefinedTVisdifferentaswewilldetailinthefollowingparagraphs.1)TVPriorInformation:TheTVregularizationtermcanbeapproximatedby

(19)

MAetal.:ANOPERATIONALSUPERRESOLUTIONAPPROACHFORMULTI-TEMPORALANDMULTI-ANGLEREMOTELYSENSEDIMAGERY117

whereistheL1norm,andisthegradientoperatorwithitsdiscreteexpressionwrittenas

(20)

andarethederivativeoperatorinthedirectionwhere

isanddirectionrespectivelyand

thedifferencebetweenthecolumn-stackvectorformofimagewhichisshiftedhorizontallybypixelsandverticallybypixelswithrespecttotheunshifted.

2)ObjectFunctionforHRImage:Applyingregularizationtechniques,theHRimagecanbesolvedbyminimizinganobject

:function

(21)

Inthefirsttermiscalledthedatafidelityterm,whichis

isbasedontheuseoftheL1norm.Thesecondtermin

theTVregularizationtermof(20),andistheregularizationparameter,whichcontrolsthetrade-offbetweenfidelitytothedataandsmoothnessofthesolution.Althoughtherearetech-niquessuchasgeneralizedcross-validation(GCV)availableforautomaticallyestimating,thesetechniquesarecomputation-allyintensiveandinmanycasestheresultsarenotsatisfactory.ThereforeinourstudywemanuallyselectedvaluesforalltheSRexperiments.

In(21)weconsidertheobstaclemaskin(1)implicitly.AnalternativewouldbetoconsiderexplicitlywiththedatafidelitybasedontheL2norm,leadingtothefollowingobjectfunction:

(22)

Wepreferobjectfunction(21)to(22)though,as(22)stronglyreliesontheaccuracyofthedetectedobstaclemask.TheL1normismorerobusttoinaccuratelyestimatedregistrationandcorruptedpixels[30].

Thereareseveralpossibilitiestominimize(21),yetofallmethodsthesteepestdescentmethodisthemoststraightfor-ward,andhasbeenacceptedasthestandardSRobjectfunc-tionminimizationmethodintheliterature[8],[9],[30].Be-isnotdifferentiableatzeroasmallconstvaluecause

(inthisexperiment)isaddedto(20)sothat

(24)

Fig.5.WorldView-2Band2imageforRiodeJaneiro(Brazil)capturedat

44.7intheforwarddirectioninJanuary2010.TheregionindicatedbythesquareshowsthetestareadefinedforSRreconstruction.

WiththesteepestdescentmethodthereconstructedHRimagecanthenbeobtainediteratively:

(25)

whereisascalardefiningthestepsizeinthedirectionofthegradientandisthedescentgradientoftheTVregular-izationterm(24),whichcanbewrittenasshownin(26),atthe

usedtoinitializethemethodisobtainedbottomofthepage.

byaresamplingschemethatisappliedtoalltheregisteredinputimages,formingameanimage(initialHRimage).Capel[15]arguedthatthisoverlysmoothapproximationtotheHRimageisextremelyrobusttonoise,thusprovidingagoodstartingpoint.

III.APPLICATION

A.ExperimentDatasets

Theapproachproposedwastestedonmulti-temporalandmulti-angleWorldView-2imagery.

WorldView-2waslaunchedinlate2009,andprovidescom-merciallyavailablepanchromaticimageryat0.5mresolutionand8-bandmultispectralimageryat1.8mresolution(

(26)

118IEEEJOURNALOFSELECTEDTOPICSINAPPLIEDEARTHOBSERVATIONSANDREMOTESENSING,VOL.5,NO.1,FEBRUARY2012

Fig.6.Multi-angleWorldView-2imageryforthestudyarea(RiodeJaneiro,Brazil):(a)44.7intheforwarddirection;(b)56.0intheforwarddirection;(c)81.4intheforwarddirection;(d)59.8inthebackwarddirection;(e)44.6inthebackwarddirection.

,and

,withcenterwavelengthsat425,480,5,

605,660,725,835,and950nmrespectively).Withinathreeminutetimeframe,fiveimageswithsatelliteelevationanglesof44.7,56.0,and81.4intheforwarddirection,and59.8and44.6inthebackwarddirection,wereacquiredforRiodeJaneiro(Brazil)inJanuary2010andthenpreprocessedasOrthoReady2Aproduct.ThepreprocesseddatawasprovidedbyDigitalGlobeforthepurposeofthe2011GRSSDataFusionContest.FromWorldView-2Band2(thebluebandwithcenterwavelengthat480nm)asub-regioncomposedof500420pixelswithlimiteddifferenceinelevationwasselectedasstudyarea(Fig.5).Thereasonwhyonlyasub-regionoftheimagewasselectedforthisexperimentisbecausetheelevationvaria-tioninsomepartsoftheimageissostrongthatitwouldmakesub-pixelimageregistrationimpossible.TheSRstudyareain-cludespartofSantosDumontAirport(coordinates:22.9076S,43.1626W).ThefiveangularimagesareshowninFig.6.Theimagecapturedatanangleof44.7intheforwarddirectionwaschosenasthereferenceLRimage.InordertomakeiteasiertosettheSRparameters,reflectancevaluesweremaximallyrescaledonthe[0,1]range.B.Results

Withtheexperimentswewanttoshowthat1)withourpro-posedSRmethodhighgeometricandphotometricregistrationaccuracycanbeobtained;2)byutilizingmulti-temporalandmulti-angleremotesensingimages,ourapproachallowsrecon-structionofhighqualityHRimageryincludingmorespatialde-tailthantheoriginaldata;3)otherstate-of-the-artSRrecon-structionalgorithmscanbeeasilyincorporatedintheproposedSRprocessingchain.

1)RegistrationAccuracy:IntheSRprocedureitisimpor-tantthatahigh-accuracyimageregistrationisattained.AsSRisperformedonnormalizedreflectancedatawithvaluesrangingbetween0and1,valuesinthedifferenceimagearebetweenand1.Afterinitialexperiments,wesetthedistributionparam-etersandin(6)and(8)as0.03and1.5respectively.

Visualinspectionoftheregistrationresultshowsagoodfitwiththereferenceimageacrossthewholeimagearea.Toevaluatetheregistrationaccuracy,weadoptedthemethodproposedin[33].Oncetheimagehadbeenregistered,wecomputedthenormalizedcrosscorrelation(NCC)overlocal

TABLEI

REGISTRATIONEVALUATIONFORWORLDVIEW-2INTERMSOFSTD

windows(3232)foreachpixelintheimage.Discardingallpoorlycorrelatedareas,wecalculatedthestandarddeviation(STD)ofthelocaldisplacements.Wealsocomparedourregis-trationresultwiththeresultobtainedbyapplyingtheapproachproposedin[17],wheretheregistrationtransformmodelwasnon-rigidTPSandCPsofhighqualitywerecollectedusingSIFTandNCC.TableIshowsthattheregistrationmethodproposedattainsacomparableorhigherregistrationaccuracythannon-rigidTPSforthisdataset.Apossibleexplanationfortheimprovedperformanceisthattheproposedmethodisaregion-basedregistrationmethodwhile[17]isinessenceafeature-basedregistrationapproach.Region-basedregistrationmethodsareconsideredtobeabletoreachhigherregistrationaccuraciesthanfeature-basedregistrationapproaches.

2)PhotometricRegistration:Theeffectivenessofthepho-tometricregistrationcanbeverifiedbymeansoftheimagehis-togramasFig.8shows.Inthegraphtheinputimagereferstotheimagetakenat59.8inthebackwarddirection,whiletherefer-enceimageistakenat44.7intheforwarddirection.Ascanbeseentheinputimageprofileistotallydifferentfromtherefer-enceimageprofile.Theoutputimageisthegeometricallyandphotometricallycorrectedversionoftheinputimage,obtainedwiththeproposedregistrationmethod.Itisclearthatafterpho-tometriccorrection,thehistogramoftheoutputimagecloselymatchesthatofthereferenceimage.Whilephotometricregis-trationcanbeperformedindifferentways,e.g.,viahistogrammatchingaftergeometricregistration,asdescribedin[16],ourapproachismoreefficientbecauseitperformsgeometricreg-istration,photometricregistrationandoutlierdetectionatthesametime.

3)ObstacleDetection:Theobstaclemapdenotespixelswhosephotometricrelationshipcannotberepresentedwithanaffinemodel.Sincethereferenceimageisfixedinourcase,weobtainfourobstaclemapsfortheWorldView-2dataset.Amajorityvotingcriterionwasusedtodefinetheobstaclemapforthereferenceimage.Foreachpixelinthereferenceimage,

denotethepixelifthemajorityoftheobstaclemaps

MAetal.:ANOPERATIONALSUPERRESOLUTIONAPPROACHFORMULTI-TEMPORALANDMULTI-ANGLEREMOTELYSENSEDIMAGERY119

Fig.7.Obstacledetection:obstaclemapfortheWorldView-2imageacquiredatanangleof44.7intheforwarddirection(referenceimage)(a),detailoftheobstaclemapforthereferenceimagepatch(b),referenceimagepatch(c),reg-isteredimagepatchfor56.0intheforwarddirection(d),81.4intheforwarddirection(e),59.8inthebackwarddirection(f),and44.6inthebackwarddirection(g).

ascorrupted,thenthepixelisconsideredascorruptedinthereferenceobstaclemap.

Fig.7showstheobstaclemapfortheWorldView-2dataset.Sincethestudyareacoveredbytheimagecontainsartificialobjectsthatcanbeeasilyrecognized,aswellasanobjectthatchangedpositionwithinthetimeframeofmulti-angleimageacquisition(airplane),hencewecanverifythedetectedout-lierseasilybycomparingtheobstaclemapwiththeregisteredmulti-angleimages.Fig.7(a)showstheobstaclemapfortheWorldView-2imagetakenunderanangleof44.7inthefor-warddirection(thereferenceimage),whileFig.7(b)showsadetailforasmallsubregionoutlinedbytheyellowsquareinFig.7(a).InordertoverifytheobstaclesdetectedinFig.7(b)theregisteredmulti-angleimagepatchescorrespondingtotheobstaclesubregionareshowninFig.7(c)–(g).Closeexamina-tionindicatesthatthepixelsintheobstaclemaskhavehighergreydisparitybetweenthereferenceimagepatchandtheimage

Fig.8.Photometricregistration:theinputimagereferstotheimagetakenat59.8inthebackwarddirection,thereferenceimageistakenat44.7intheforwarddirection,theoutputimageisthegeometricallyandphotometricallycorrectedversionoftheoriginalimage.

Fig.9.SRresultforWorldView2multi-angleimagery(Band2),withsubre-gionsfordetailedanalysisoutlinedinyellow,green,red,blueandcyan.

patchesfromtherestviewangles.Ascanbeobserved,mostofthecorruptedpixelsdetectedinthissmallsubregioncanbere-latedtoanairplaneshiftingpositionduringthetimeofimageacquisition.

InourSRapproach,theobstaclemasklabelspixelsthatcannotbeusedforestimatingtheregistrationmodel.Theob-staclemask,however,canalsobeemployedduringHRimage

120IEEEJOURNALOFSELECTEDTOPICSINAPPLIEDEARTHOBSERVATIONSANDREMOTESENSING,VOL.5,NO.1,FEBRUARY2012

Fig.10.Detail1ofWorldView-2SRresult:(a)referenceLRimage;(b)bicubicimage;(c)POCS;(d)L1normdeconvolution;(g)proposedmethod;(h)panchromaticimage(spatialresolution0.5m).

Huberregularizer;(e)VPLR;(f)non-uniform

Fig.11.Detail2ofWorldView-2SRresult:(a)referenceLRimage;(b)bicubicimage;(c)POCS;(d)L1normdeconvolution;(g)proposedmethod;(h)panchromaticimage(spatialresolution0.5m).

Huberregularizer;(e)VPLR;(f)non-uniform

reconstructionasitprovidesinformationonthereliabilityofpixelsintheLRimages.

4)ComparisonofSRResults:Toshowthatourprocessingchainisoperationallyoriented,itisimportanttodemonstratethatdifferentstate-of-the-artSRreconstructionmethodscanbeincorporatedintheproposedframework.WewillalsocompareourTV-basedSRreconstructionmethodwithotherSRalgo-rithms.Tothisend,threepopularSRreconstructionalgorithmswereimplemented:theprojectionsontoconvexsets(POCS)method[11],theL1normandHuberregularizer[9]andtheVPLRalgorithm[12].TheseSRalgorithmssharethesamereg-istrationprocedure,theonlydifferencebeingtheHRimagere-constructionpart:

1)ForthePOCSalgorithm,whenthedifferencebetweensim-ulatedLRimages(fromthereconstructedHRimage)andobservedLRimagesbecomessmall,theiterationwillstopandthepotentialHRimageisconstructed.

2)TheL1normandHuberregularizermethodisverysimilartoourmethodasitsharesthesamedatafidelitytermtype.Italsoadoptsthesteepestdescentmethodtominimizetheobjectfunction.TheonlydifferenceliesintheregularizertermasitemploysaHuberpriorratherthanaTVprior.3)TheVPLRalgorithmisanon-uniforminterpolationbasedSRreconstructionmethod.UnliketheotherSRreconstruc-tionmethodsitdoesnotdealwiththeblurringintroducedbythePSFintheimagereconstructionphase.However,itisoneofthefewSRalgorithmsthatcanincorporateanob-staclemaskduringthereconstructionprocess.

Finally,wealsocomparedtheSRmethodproposedinthisstudywiththenon-uniforminterpolationanddeconvolution

MAetal.:ANOPERATIONALSUPERRESOLUTIONAPPROACHFORMULTI-TEMPORALANDMULTI-ANGLEREMOTELYSENSEDIMAGERY121

Fig.12.Detail3ofWorldView-2SRresult:(a)referenceLRimage;(b)bicubicimage;(c)POCS;(d)L1normdeconvolution;(g)proposedmethod;(h)panchromaticimage(spatialresolution0.5m).

Huberregularizer;(e)VPLR;(f)non-uniform

Fig.13.Detail4ofWorldView-2SRresult:(a)referenceLRimage;(b)bicubicimage;(c)POCS;(d)L1normdeconvolution;(g)proposedmethod;(h)panchromaticimage(spatialresolution0.5m).

Huberregularizer;(e)VPLR;(f)non-uniform

methodwhichwepresentedin[13],whereDelaunaytriangula-tion-basednon-uniforminterpolationisfollowedbyaWienerfiltertoremovetheblurringafterLRimagepixels’positionshavebeendeterminedwithaTPSmodel.

FortheSRexperimentthemagnificationfactorwassetto

,whichleadstoquadrupledpixelnumbers.

ThestandarddeviationsoftheellipticalPSFforeachobservedimagerangefrom0.45pixelsto0.55pixelsfortheWorldView-2dataset.ThemaximumnumberifiterationsfortheproposedSRalgorithmwassetat100,thethresholdforthesuccessivediffer-.ThesteepdescentenceimageL2normwassetat

stepsizein(25)issetas0.001,andtheregularizationparam-eterwaschosenequalto2bytrial-and-error.ForotherSRal-gorithmstheparametersettingswerebasedontheinstructionsdetailedinthecorrespondingpapers[9],[11]–[13].

AsWorldView-2suppliesapanchromaticimagewhosespa-tialresolutionisfinerthanthatofthemulti-spectralimages,thepanchromaticimagecanbeusedasgroundtruthforvisualas-sessmentoftheSRalgorithmsasMerinoandNunezdidin[12].Fig.9showstheSRresultfortheWorldView-2studyareaasawhole.Amoredetailedvisualcomparisonforseveralsub-regions,highlightedinFig.9,isshowninFigs.10to14.ForeachsubregiontheLRreferenceimage,asingle-framebicubicSRenhancedimage,SRresultsfromdifferentmethodsandthepanchromaticimageareshown.ComparedtotheLRimageandthebicubicSRimage,improvementswithmulti-frameSRare

122IEEEJOURNALOFSELECTEDTOPICSINAPPLIEDEARTHOBSERVATIONSANDREMOTESENSING,VOL.5,NO.1,FEBRUARY2012

Fig.14.Detail5ofWorldView-2SRresult:(a)referenceLRimage;(b)bicubicimage;(c)POCS;(d)L1normdeconvolution;(g)proposedmethod;(h)panchromaticimage(spatialresolution0.5m)

Huberregularizer;(e)VPLR;(f)non-uniform

obviouswhencomparingthedifferentresultsobtainedwiththepanchromaticdata.Forexample,inFig.10theninespotsaroundthecircularfeatureontheairporttarmacaredifficulttodetectintheLRimagepatchandinthebicubicimage,whileintheSRimagepatchestheninespotsareclearlypresent.ThecircularshapeinFig.11isnotrecognizableintheLRimagepatchwhileitsshapeiswellpronouncedintheSRimages.CharacteristictexturesandlinearstructuresarealsoeasiertoidentifyintheSRimagesthanintheLRimageasFig.12shows.InFig.13wecanseethatviaSRtechniquessomecharacteristicairplanefeatureslostintheLRimagearesuccessfullyretrieved.

Ingeneralallmulti-frameSRalgorithmstestedhelpincreasethelevelofspatialdetailintherepresentationofsurfaceob-jects.However,wealsonoticethatsomeSRalgorithmscreateartifactsthatdonotcorrespondwiththegroundtruthasshowninthepanchromaticimage.IfwecomparethePOCSSRre-sultinFig.11(c)withthegroundtruthaswellaswithotherSRenhancedimages,wenoticethatthenoisearoundtherecon-structedcircularobjectisexaggerated.ThismightbeexplainedbythefactthatofalltheSRreconstructionmethodsappliedinthisstudyPOCSusestheleastpriorinformationaboutthepo-tentialHRimage.POCSSRisalsoverysensitivetocorruptedpixelsasthismethoddoesnotconsidertheeffectofobstacles.TakingthefiveLRimagepatchesinFig.7(c)–(g)asanex-ample,SRalgorithmssuchasVPLR,theL1normandHuberregularizermethod,andourproposedmethodareallabletoselectpixelsintelligentlyandavoidtheuseofcorruptedLRpixelsbyemployingtheobstaclemask(Fig.7(b))eitherex-plicitlyorimplicitly.POCSSR,however,treatseachLRpixelequally.Hencecorruptedpixelswillhaveanegativeeffectonthefinalreconstructedimage.Assuch,thecentralairplaneintheupperpartofthereferenceimagepatch(Fig.7(c)),canstillbetracedinthereconstructedPOCSSRresult,asFig.14(c)illustrates.Sensitivitytocorruptedpixelsisalsoaproblemforthenon-uniforminterpolationanddeconvolutionSRmethod,asFig.14(f)shows,althoughcomparedtothePOCSmethodtheimpactisless.

VPLRdoesnotconsidertheeffectofPSFduringtheSRreconstructionprocedure.Therefore,ifwecomparetheVPLRresultwiththeresultobtainedwithotherSRmethodswecaneasilynoticethatitsresolutionenhancementcapa-bilityiscomparativelyweakwithlessretrieveddetails(seeFigs.10(e)to14(e)).Theadvantageofthismethodisitsinsensitivitytonoisecomparedtoothermethods.Sincenodeconvolutionisinvolved,itwillnotintroduceanyartifactsduringtheHRimagereconstructionprocedure.Itisalsotoleranttocorruptedpixelsasitemploystheobstaclemask.However,theaccuracyofthedetectedobstaclemaskwillbecriticaltoitsperformance.Also,inthecasewhereonlyafewusefulLRpixelsareavailablearoundthetargetpixelintheHRgrid,itsHRimagereconstructioncapabilitywillbegreatlyreducedsincenoaprioriinformationaboutthepotentialHRimageisutilized.This,combinedwiththereducedresolutionenhancementcapabilityduetoitsinabilitytoreversetheeffectofPSF,mayexplainwhyinFig.14(e)VPLRSRisincapableofproperlyreconstructingtheshapeoftheairplane.

TheL1normandHuberregularizermethodandthemethodproposedinthisstudyproduceverysimilarresults.BothmethodsarelesssensitivetonoisebecauseoftheHuberregu-larizerandtheTVregularizer,andarerobusttocorruptedpixelsaswellduetotheL1normbaseddatafidelityterm.BasedonvisualcomparisonitishardtotellwhichalgorithmissuperiorasboththeHuberregularizerandtheTVregularizerpreservetheedgeinformationinthereconstructedimage.However,comparedtotheL1normandHuberregularizermethodourapproachhaslessparameterstoset.IngeneralwecansaythatoverallinthisstudybothmethodsproducethebestSRresults,asshowninFigs.10to14.

MAetal.:ANOPERATIONALSUPERRESOLUTIONAPPROACHFORMULTI-TEMPORALANDMULTI-ANGLEREMOTELYSENSEDIMAGERY123

IV.CONCLUSION

ThispaperproposesanoperationalSRapproachformulti-temporalandmulti-angleremotesensingimagery.Thecontri-butionsofthispaperarefourfold:

1)Anewobservationmodelforremotesensingimagery,usefulasabasisforoperationalSRimagereconstructionisproposed.

2)Basedonthenewobservationmodelanoperational,patch-basedSRapproachispresentedtacklingtwoofthemaindifficultiesinSR:localdistortioncausedbysmalltopographiceffectsand/orplatforminstability,andpresenceofcorruptedpixelscausedbyobstaclesintheimagery,makingsimpleapplicationofstate-of-the-artSRalgorithmsimpossible.

3)Thepatch-basedSRschemeproposedissuitableforpro-cessinglargesizeremotesensingimages,andtheimple-mentationofthe

matrixisconsistentwiththeimagingmodel.Inordertoavoidblockeffects,asimplebuteffec-tivefusionschemehasbeenproposed.

4)TheTV-basedSRreconstructionmethodisintroducedwithitsobjectfunctionefficientlybeingsolvedbythesteepestdescentmethod.

Therearesomelimitationstothemethodproposedthough.Althoughtheobjectiveofourstudywastoproposeagener-allyapplicableSRapproachforremotesensingimagery,themethodwillnotworkwellforareaswithmoderatetostrongvariationsinelevation.Insuchcasesitwillbemoreappropriatetoapplyamethodthatstartsfromortho-rectifiedimagery.Also,theselectionofaproperimagepatchsizeforSRreconstructionneedsfurtherinvestigation.Theknife-edgeinputPSFestima-tionmethod,whichwasadoptedinthispapertoestimatethePSFofLRimages,alsohasitslimitations.HowtheuseofthisPSFestimationmethodaffectsthequalityofthereconstructedimagerequiresfurtherwork.AnotherinterestingissuefromanapplicationpointofviewiswhetherSRimagereconstructionmayleadtoimprovementsinimageclassificationperformance.Thisiscurrentlybeinginvestigated.

ACKNOWLEDGMENT

TheauthorswouldliketothankDr.LyndseyC.PickupfromtheVirtualRealityResearchGroupofthePhysiologyDepartmentofOxfordUniversityandDr.PascalGetreuerfromtheCentredeMathématiquesetdeLeursApplicationsoftheEcoleNormaleSupérieuredeCachanforvaluablediscussionsontheimplementationofthealgorithmandforprovidingusefulsoftware.

REFERENCES

[1]T.Akgun,Y.Altunbasak,andR.M.Mersereau,“Super-resolutionre-constructionofhyperspectralimages,”IEEETrans.ImageProcess.,vol.14,no.11,pp.1860–1875,Nov.2005.

[2]S.C.Park,M.K.Park,andM.G.Kang,“Super-resolutionimagere-construction:Atechnicaloverview,”IEEESignalProcess.Mag.,vol.20,no.3,pp.21–36,2003.

[3]P.Vandewalle,S.Susstrunk,andM.Vetterli,“Afrequencydomainap-proachtoregistrationofaliasedimageswithapplicationtosuper-res-olution,”EURASIPJ.Appl.SignalProcess.,vol.2006,2006,ArticleID71459.

[4]R.Y.TsaiandT.S.Huang,“Multiframeimagerestorationandreg-istration,”Adv.ComputerVisionImageProcess.,vol.1,pp.317–339,1984.

[5]S.P.Kim,N.K.Bose,andH.M.Valenzuela,“Recursivereconstruc-tionofhighresolutionimagefromnoisyundersampledmultiframes,”IEEETrans.Acoust.,Speech,SignalProcess.,vol.38,no.6,pp.1013–1027,June1990.

[6]H.UrandD.Gross,“Improvedresolutionfromsub-pixelshiftedpic-tures,”CVGIP:GraphicalModelsandImageProcessing,vol.,pp.181–186,Mar.1992.

[7]M.IraniandS.Peleg,“Improvingresolutionbyimageregistration,”in

CVGIP:GraphicalModelsandImageProc.,May1991,vol.53,pp.231–239.

[8]F.Li,X.Jia,D.Fraser,andA.Lambert,“Superresolutionforremote

sensingimagesbasedonauniversalHiddenMarkovTreemodel,”IEEETrans.Geosci.RemoteSens.,vol.48,pp.1270–1278,Mar.2010.[9]H.Shen,M.Ng,P.Li,andL.Zhang,“Super-resolutionreconstruction

algorithmtoMODISremotesensingimages,”Comput.J.,vol.52,no.1,pp.90–100,2009.

[10]C.LatryandB.Rouge,“Superresolution:Quincunxsamplingandfu-sionprocessing,”inProc.Int.GeoscienceandRemoteSensingSymp.2003(IGARSS’03),Jul.2003,vol.1,pp.315–317.

[11]A.Galbraith,J.Theiler,K.Thome,andR.Ziolkowski,“Resolution

enhancementofmulti-lookimageryforthemultispectralthermalimager,”IEEETrans.Geosci.RemoteSens.,vol.43,pp.19–1977,2005.

[12]M.T.MerinoandJ.Nunez,“Super-resolutionofremotelysensedim-ageswithvariable-pixellinearreconstruction,”IEEETrans.Geosci.RemoteSens.,vol.45,no.5,pp.1446–1457,2007.

[13]J.C.-W.Chan,J.Ma,P.Kempeneers,andF.Canters,“Superresolution

enhancementofhyperspectralCHRIS/Probaimageswithathin-platesplinenonrigidtransformmodel,”IEEETrans.Geosci.RemoteSens.,vol.48,pp.2569–2579,May2010.

[14]M.EladandA.Feuer,“Restorationofasinglesuperresolutionimage

fromseveralblurred,noisy,andundersampledmeasuredimages,”IEEETrans.ImageProcess.,vol.6,no.12,pp.16–1658,Dec.1997.[15]D.P.Capel,“Imagemosaicingandsuper-resolution,”Ph.D.disserta-tion,Univ.ofOxford,Oxford,U.K.,2001.

[16]M.GevrekciandB.K.Gunturk,“Superresolutionunderphotometric

diversityofimages,”EURASIPJ.Adv.SignalProcess.,vol.2007,2007,ArticleID36076.

[17]J.Ma,J.C.-W.Chan,andF.Canters,“Fullyautomaticsub-pixelimage

registrationofmulti-angleCHRIS/Probadata,”IEEETrans.Geosci.RemoteSens.,vol.48,no.7,pp.2829–2839,July2010.

[18]W.Y.ZhaoandS.Sawhney,“Issuper-resolutionwithopticalflow

feasible?,”inProc.ECCV,LNCS2350,2002,pp.599–613.

[19]G.Yang,C.V.Stewart,M.Sofka,andC.-L.Tsai,“Registrationof

challengingimagepairs:Initialization,estimation,anddecision,”IEEETrans.PatternAnal.Mach.Intell.,vol.29,no.11,pp.1973–19,Nov.2007.

[20]H.TrusselandB.Hunt,“Sectionedmethodsforimagerestoration,”

IEEETrans.Acoust.,SpeechSignalProcess.,vol.ASSP-26,no.2,pp.157–1,1978.

[21]Q.TianandM.N.Huhns,“Algorithmsforsubpixelregistration,”

Comput.Vision,Graphics,ImageProcess.,vol.35,pp.220–233,Aug.1986.

[22]S.PeriaswamyandH.Farid,“Medicalimageregistrationwithpartial

data,”Med.ImageAnal.,vol.10,no.3,pp.452–4,Jun.2006.

[23]X.Feng,“Analysisandapproachestoimagelocalorientationestima-tion,”M.S.thesis,Dept.Comput.Eng.,Univ.California,SantaCruz,Mar.2002.

[24]M.A.Martín-Fernández,M.Martín-Fernández,andC.Al-berola-Lopez,“Alog-euclideanpolyaffineregistrationforarticulatedstructuresinmedicalimages,”inProc.MICCAI,2009,pp.156–1.[25]A.Mohammad-Djafari,“Super-resolution:Ashortreview,Anew

methodbasedonhiddenMarkovmodelingofHRimageandfuturechallenges,”Comput.J.,vol.52,no.1,pp.126–141,2009.

[26]F.Chen,J.Ma,J.C.-W.Chan,andD.Yan,“Quantitativemeasurement

ofthehomogeneityandcontrastofstepedgesintheestimationofthepointspreadfunctionofasatelliteimage,”Int.J.RemoteSens.,vol.32,no.22,pp.7179–7201,2011.

[27]S.BakerandT.Kanade,“Limitsonsuper-resolutionandhowtobreak

them,”IEEETrans.PatternAnal.Mach.Intell.,vol.24,no.9,pp.1167–1183,Sep.2002.

[28]L.C.Pickup,“Machinelearninginmulti-frameimagesuper-resolu-tion,”Ph.D.dissertation,Univ.Oxford,Oxford,U.K.,Feb.2008.

124IEEEJOURNALOFSELECTEDTOPICSINAPPLIEDEARTHOBSERVATIONSANDREMOTESENSING,VOL.5,NO.1,FEBRUARY2012

[29]M.BerteroandP.Boccacci,IntroductiontoInverseProblemsin

Imaging.Bristol,U.K.:IOPPublishing,1998.

[30]S.Farsiu,D.Robinson,M.Elad,andP.Milanfar,“Fastandrobust

multi-framesuperresolution,”IEEETrans.ImageProcess.,vol.13,no.10,pp.1327–1344,Oct.2004.

[31]M.K.Ng,H.Shen,L.Zhang,andE.Lam,“Atotalvariationbased

super-resolutionreconstructionalgorithmfordigitalvideo,”EURASIPJ.Adv.SignalProcess.,vol.2007,2007,ArticleID74585.

[32]M.K.Ng,H.Shen,S.Chaudhuri,andA.C.Yau,“Zoom-basedsuper-resolutionreconstructionapproachusingpriortotalvariation,”Opt.Eng.,vol.46,no.12,2007.

[33]S.Leprince,S.Barbot,F.Ayoub,andJ.P.Avouac,“Automaticandpre-ciseorthorectification,coregistration,andsubpixelcorrelationofsatel-liteimages,applicationtogrounddeformationmeasurements,”IEEETrans.Geosci.RemoteSens.,vol.45,no.6,pp.1529–1558,Jun.2007.

JianglinMareceivedtheB.S.degreefromXidianUniversity,Xi’an,China,in2003andtheM.S.de-greefromtheInstituteofRemoteSensingApplica-tions,ChineseAcademyofSciences,Beijing,China,in2006.Currently,whilepursuingthePh.D.degree,heisanAssociateResearcherwiththeCartographyandGISResearchGroup,DepartmentofGeography,VrijeUniversiteitBrussel,Brussels,Belgium.

Hisresearchinterestsareimageprocessinganditsapplicationinremotesensing.

JonathanCheung-WaiChanreceivedtheB.S.andM.S.degreesingeographyfromtheChineseUni-versityof,Shatin,,andthePh.D.degreefromtheCenterofUrbanPlanningandEnvironmentalManagement,University,,in1999.

Between1998and2001,hewasaresearchscientistattheGeographyDepartment,UniversityofMaryland,CollegePark,UnitedStates.Since2001,hehasbeenworkingatVrijeUniversiteitBrussel(VUB),Brussels,Belgium.Hismainre-searchinterestsaremachinelearningalgorithms,landcoverclassificationusinghyperspectraldata,spatialenhancementofhyperspectralimagesusingsuperresolutionrestoration,andenvironmentalmonitoringusingtemporalremotesensingimages.HeisnowwiththeDepartmentofElectronicsandInformaticsatVUB.

FrankCantersreceivedtheDiplomadegreesingeographyandappliedcomputerscience,andthePh.D.degreeinsciencesfromtheVrijeUniversiteitBrussel,Belgium.

Currently,heisAssociateProfessorintheDe-partmentofGeographyoftheVrijeUniversiteitBrussel,whereheisheadingtheCartographyandGISResearchGroupandisteachingcartography,remotesensingandgeo-informationscience.Since2001heisalsoaVisitingProfessorintheGeographyDepartmentofGhentUniversity.Hismainresearch

interestsareurbanremotesensing,multi-sensor/multi-resolutionimageanal-ysis,andmodelingofuncertaintyinspatialdata.