Statistical Track-Before-Detect Methods Applied to Faint Optical Observations of Resident Space ...

14 Pages·2015·2.87 MB·English

StatisticalTrack-Before-DetectMethodsAppliedto
FaintOpticalObservationsofResidentSpaceObjects
KoheiFujimoto
UtahStateUniversity,UnitedStates
MasahikoUetsuhara
Independentresearchengineer,Japan
ToshifumiYanagisawa
JapanAerospaceExplorationAgency,Japan
ABSTRACT
In this paper, we apply a statistically rigorous track-before-detect (TBD) method, the Bernoulli particle
filter,toactualimageryofresidentspaceobjects(RSOs). Robustmethodsinthisrealmwillleadtobetter
space domain awareness (SDA) while reducing the cost of sensors and optics. We focus on estimating
sensor-level kinematics of RSOs for low signal-to-noise ratio (SNR) short-arc observations. For four
minutearcsof29images, targetswithSNRaslowas2.6983aretrackedto3arcsecaccuracyinimage
planepositionand0.07arcsec/secinvelocity.
1. INTRODUCTION
Automated detection and tracking of faint objects in optical, or bearing-only, sensor imagery is a topic of immense
interest in space surveillance [7,32,40,48,56,62]. Robust methods in this realm will lead to better space domain
awareness (SDA) while reducing the cost of sensors and optics [5]. They are especially relevant in the search for
high area-to-mass ratio (HAMR) objects, as their apparent brightness can change significantly over time [45]. A
track-before-detect(TBD)approachhasbeenshowntobesuitableforfaint, lowsignal-to-noiseratio(SNR)images
of resident space objects (RSOs). TBD does not rely upon the extraction of feature points within the image based
onsomethresholdingcriteria, butratherdirectlytakesasinputtheintensityinformationfromtheimagefile[8,44].
Notonlyisalloftheavailableinformationfromtheimageused, TBDavoidsthecomputationalintractabilityofthe
conventionalfeature-basedlinedetection(i.e.,“stringofpearls”)approachtotrackdetectionforlowSNRdata.
ImplementationofTBDrootedinfinitesetstatistics(FISST)theoryhasbeenproposedrecentlybyVo,etal.[36,37,58]
Compared to other TBD methods applied so far to SDA, such as the stacking method or multi-pass multi-period
denoising, the FISST approach is statistically rigorous and has been shown to be more computationally e�cient,
thus paving the path toward on-line processing. In this paper, we apply a Bernoulli filter to actual CCD imagery of
RSOs [42]. The Bernoulli filter can explicitly account for the birth and death of a target in some measurement arc.
TBDisachievedviaasequentialMonteCarlo(SMC)implementation. Forfourminutearcsof29images,targetswith
SNRaslowas2.6983aretrackedto3arcsecaccuracyinimageplanepositionand0.07arcsec/secinvelocity.
Although the advent of fast-cadence scientific CMOS sensors have made the automation of faint object detection a
realistic goal, it is nonetheless a di�cult goal, as measurements arcs in space surveillance are often both short and
sparse [16,50]. FISST methodologies have been applied to the general problem of SDA by many authors, but they
generallyfocusontrackingscenarioswithlongarcsorassumethatlinedetectionistractable[3,4,10,11,19,20,24–
26,28–30,33]. Wewillinsteadfocusthisworkonestimatingsensor-levelkinematicsofRSOsforlowSNRtoo-short
arcobservations. Oncesaidestimateismadeavailable,trackassociationandsimultaneousinitialorbitdetermination
maybeachievedviaanynumberofproposedsolutionstothetoo-shortarcproblem,suchasthoseincorporatingthe
admissibleregion[12,17,39,43,47,51,54,61]. ThebenefitofcombiningFISST-basedTBDwithtoo-shortarcasso-
ciationgoesbothways; i.e., theformerprovidesconsistentstatisticsregardingbearing-onlymeasurements, whereas
thelattermakesbetteruseoftheprecisedynamicalmodelsnominallyapplicabletoRSOsinorbitdetermination.
Theoutlineofthispaperisasfollows. WefirstmotivatetheuseofTBDforspacesurveillancewithopticalsensors
(Section2). Then,wediscussthecurrentliteratureonFISSTapplicationstoSDAandTBD(Section3). Afterabrief
introductionoftheBernoulliparticlefilter(Section4),the“fineprint”ofTBDencounteredduringtheimplementation
aredetailed,withafocusonimprovingfilterperformancethroughtunableparametersandthemeasurementlikelihood
function(Section5). Finally,filterresultsareanalyzedforthreetracksofvaryingSNR(Section6).
2. SENSOR-LEVELTRACKINGINSDA
Traditionally, thedetectionofastronomicalobjectsintelescopeimageryhasbeenachievedthroughadetect-before-
track(DBT)approach,inwhichagroupofpixelsaboveauser-definedintensitythresholdareassignedtoacorrespond-
ing hypothesized object [46]. In more modern algorithms, one may additionally require that the pixels follow some
point spread function (PSF) [64] or assign hypothesized objects based on the signal’s local energy instead of a uni-
formthreshold[32]. Ifsomesubsetofthesehypothesizedobjectsfollowanadmissiblekinematicmotion–nominally,
linearconstantvelocitymotionforRSOs–thentheobjectsaresaidtoformatrackortracklet. Althoughthis“string
of pearls” approach has been very successful thus far and certainly has its benefits, such as drastically compressing
the information from each image to a small set of coordinates, it also could become a computational bottleneck in
space surveillance should either the total number of images or hypothesized objects in a track drastically increase.
Quantitatively,therun-timecomplexityofDBTisorderO(nm) ifnhypothesizedobjectsareassignedineachimage
foratrackofm images. Anincreaseinbothm andnarerealisticconcernsinSDAasthesensitivityand,especially
withtheadventofscientificCMOSsensors,thetemporalresolutionofimagesensorsimproverapidlyinthecoming
years. Furthermore, there is an inherent loss of information by treating detection as a boolean decision, rather than
aprobabilisticone,aswellasbyseparatingitfromsubsequentSDAtaskssuchastracking,characterization,andthe
deliveryofactionableinformation.
Analternativeapproachistrack-before-detect(TBD),inwhichintensityinformationfromentireimagesareprocessed
insteadofsolelytheircharacteristicpoints.TherearethreemajoradvantagestoTBD.First,therun-timecomplexityis
orderO(m);i.e.,itisnolongerdependentonthenumberofobjectsintheimageandonlyscaleslinearlywithrespect
tothenumberofimagesinthetrack. Itispossibletoreapthebenefitsofthislineargrowthevenwithsmallvaluesof
m. Mostnotably,thedimensionalityofthekinematicsisfour(twopositionstatesandtwovelocitystatesintheimage
plane)insteadofthesixinthethree-dimensionalorbitproblem,makinghighlyrobustimplementationsofTBD,like
SMC,aviableoption. Infact,althoughonlyproposedasfutureworkinthispaper,weenvisionreal-timeprocessing
of fast-cadenceimagery ofRSOs byemploying general-purposecomputing ongraphics processingunits (GPGPU),
towhichTBDishighlyamenable.
Second,thestateuncertaintyofRSOsisrigorouslycharacterizedatthetimeofmeasurement.Thisinformationmaybe
carriedonthroughouttheSDAworkflowvianumeroustechniquesstudiedinrecentyears[9,18,21,22,27,31,38,57],
aswellasinformthedevelopmentofbetterphysicalmodelsformeasurementerror. Suchconsistentquantificationof
RSOstateuncertaintystartingatthesensor-levelisintegralto“improv(ing)spacesafety.”[15]
Finally, as TBD methods infer the existence of objects based on multiple images rather than assign hypothesized
objects per image, they have been shown to excel when SNR is low [8,44,58]. A study of the Pareto frontier for
a Raven-class telescope suggested that “the SNR threshold required by the detection algorithm largely influences
theoveralldetectioncapabilityofthesystem.”[5]ImproveddetectionatlowSNRwillalsoenableminimalistpassive
sensorstomakesubstantialcontributionstoSDA[48].Therefore,TBDwillgreatlyenhancetheutilityofbothexisting
andfutureSDAopticalsensorswithoutsignificant,ifany,hardwareinvestment.
InitialworkinTBDforopticalobservationsofRSOsemergedfromthestackingmethodproposedbyYanagisawa,et
al., where low SNR RSOs are surfaced through a brute-force search for its motion in the image plane [63]. Images
inatrackaresuccessivelydisplacedbysomeconstantamount. IfanRSOintheimageplanemovesatthesamerate
as the displacements, then the median intensity value in its vicinity should be higher than the background. In order
to reduce the number of displacement hypotheses that must be tested per track, the method has been implemented
on an FPGA [62] as well as combined with a priori information from RSO population models [55]. More recently,
a Bayesian estimation approach was taken in [56], where the measurement likelihood function of RSOs was first
determinedbyageneticalgorithm. RSOsweresubsequentlytrackedwithaparticlefilter. Independently,thebanked
matchedfilter(MF)wasstudiedbyMurphy,etal.[40]. Here,beginningwithpartialaprioristateinformation,such
asashort-arcmeasurement,onemaygenerateabankoftemplatesregardingwhatonecanexpecttoimageatalater
time. When a convolution with the actual observations is computed, templates that resemble the true state result in
anincreaseinSNR;asubsequenthypothesistestgivestheprobabilityofdetection. Asthesetemplatesareinherently
assigned to only a localized area in the sensor frame, the banked MF is computationally e�cient. Furthermore, the
aforementionedprobabilityofdetectioncandirectlybeinsertedasaPDFintoaBayesianestimator.Anotherapproach
proposed by Dao, et al. is based on the concept of random sampling and consensus (RANSAC) [7]. Detection is
achieved by searching for a statistically significant series of characteristic points, chosen according to user-defined
criteriaonthePSFandintensity,whichalignlinearlyinadatacube.
The goal in this work is to combine the stacking method’s a priori state free generality with the statistical rigor of
filtering approaches. We apply a new family of estimation algorithms called random finite set (RFS) filters. In the
nextseveralsections,weoutlinetheirbenefitsandhowtheyhavebeenappliedtoSDAthusfar. Wethenfocusonone
particularflavor: theBernoullifilter.
3. FISSTASAPPLIEDTOSDAANDTBD
One important assumption in Bayesian estimation, which has been the workhorse of orbit determination since the
beginning of spaceflight [53], is that the measurements being processed are of and only of the object being tracked.
Therefore,dataassociationisusuallytreatedbyaseparateheuristicpreprocessor.Inmanyscenariospertinenttospace
surveillance,however,thisassumptionisnolongervalid;notallRSOsarebrightenoughtoproduceaconsistentsignal
fortheentiretyofatrack,forinstance. Inaddition,dataassociationbyitselfcanbecomecomputationallyintractable
unlessassociationhypothesesareproperlygated[1,52].
RFS filters, which operate on unordered finite sets, are one proposed solution to these problems. RFSs are fully
specifiedbythedistributionofthenumberofelementsintheset(cardinality)andthejointdistributionoftheelements
conditioned upon the cardinality. As such, RFS filters can explicitly account for multiple tracking targets or, more
importantly for TBD, the complete lack thereof. Furthermore, one may retain notions such as probability density
functions(PDFs),theirstatisticalmoments,andBayesruleforRFSsviatheFISSTframework;e.g.,theFISSTPDF
ofanRFSX= x1,...,xn comprisedofnelementsis[42]
{ }
f(X) =n! ⇢(n) pn(x1,...,xn), (1)
· ·
where ⇢isthecardinalitydistributionand pn isthejointdistributionoftheelements. Wereferreaderstoashortlist
ofexistingliteratureforfurtherdetails[34–36,42].
PerhapsthemostdirectapplicationofFISST-basedmethodstoSDAistoassignanRSO’sstatetoeachsetelement.
As the exact multi-Bayesian RFS filter is often intractable, many approximations exist and have been applied to
SDA[11,24],includingbutnotlimitedtotheGM-PHD[10],GM-CPHD[19,29],ISP[4],andGLMBfilters[30]. If
welooktosurveypapers,comparisonsbetweendi↵erentRFSfiltersandcompetingmulti-targetmulti-sensortracking
methodologiesarepresentedin[3,20,33],andchallengesspecifictotheSDAproblemareoutlinedin[28]. Finally,
solving the joint sensor management and allocation problem using FISST is discussed in [26,65,66], and the joint
classificationproblemin[25].
TheBernoullifilterisaRFSfilterforwhichthecardinalitydistributionisgivenasaBernoullidistribution;thatis,we
specifyandsubsequentlyestimatetheprobability q thatatargetexists[42]. Itthusallowsforanexplicitdistinction
in the measurement likelihood function based on whether or not an RSO is in the field of view. Furthermore, it
maybeextendedtoamulti-targetproblembydefiningamulti-Bernoulliorlabeledmulti-BernoulliRFS[35,59,60].
ThefamilyofBernoullifiltershassuccessfullybeenappliedtoTBD[23,41,58]; thisproventrackrecordisanother
motivationforustoapplytheBernoullifiltertoactualCCDimageryofRSOs.
Finally, wenoteherethat, inthispaper, wechoosetofocusonsensor-leveltrackingasweareremindedofthelack
of astrometric information one gains from short-arc optical observations of RSOs [16]. Even as sensor technology
advances, survey designs that emphasize coverage will continue to produce tracks that cover only a minute fraction
of an object’s orbit. On the other hand, it has been shown that additional information such as the estimated rate-of-
changeofthesensorbearing[12,17,54,61],integralinvariantsoftheorbitalmotion[39],orspecialsolutionstothe
Lambertproblem[43,47,51]allowsonetogreatlyconstrainthemeasurementlikelihoodfunction. Assuch,insteadof
a“top-down”approachwheretheimmediatepurposeofRFSfilteringistoholisticallytrackandcharacterizetheRSO
population,wetakea“bottom-up”approachwhereweinitiallyaimtodistributemeasurementinformationtoalimited
yet su�cient set of states. As time passes and more measurement or dynamical information becomes available, the
problemcanbeconstrainedinhigherdimensionssoastoaddmorestatesandparameterstoourknowledgeofRSOs.
We present RFS filters, then, as a means to e↵ectively conduct this initial distribution of information – e↵ective
both in terms of computational e�ciency and low SNR performance – as well as gain probabilistic metrics on said
e↵ectiveness.
4. BERNOULLIPARTICLEFILTERS
WeintroducethepredictionandupdateequationsfortheBernoullifilter;amorethoroughtreatmentisgivenin[42].
Assumptionsmadeinthederivationarelistedbelow
• ThemeasurementsareassumedtobestatisticallyindependentofthestateRFS
• ThestateRFSdynamicsisaMarkovprocesswithaknowntransitionaldensity
• TheinitialstatePDF,initialtargetexistenceprobability,andtargetbirth/survivalprobabilitiesareknown
Giventhatmeasurementstakenattimet0throughtk 1havebeenprocessedtoobtainstatex0,thepredictionandupdate
�
stepstoobtaintheprobabilityofexistenceofatargetobjectqandaPDFs(x)forthestatexis
Prediction
qk k 1 = pb(1 qk 1k 1)+(1 ps)qk 1k 1 (2)
| � � � | � � � | �
sk k 1(x) = pb(1�qk�1|k�1)bk|k�1(x) + psqk�1|k�1 ⇡k|k�1(x|x0)sk�1|k�1(x0)dx0, (3)
| � qk k 1 R qk k 1
| � | �
where pb is the probability of a target being born duringtk 1 totk, ps is the probability of a target surviving
�
during the same timeframe, b(x) is the PDF of the target birth process, and ⇡(xx0) is the Markov process
|
transitionaldensity.
Update
lk(zk xk) = n g1(i) fzk(i)|xkg (4)
| g(i) z(i)
Yi=1 0 k
f g
Ik = lk(zk x)sk k 1(x)dx (5)
| | �
Z
q I
k k 1 k
qk k = | � (6)
| 1 qk k 1+qk k 1Ik
� | � | �
sk k(x) = in=1g1(i) fzk(i)|xgsk|k�1(x) , (7)
| Qin=1g1(i) zk(i)|x sk|k�1(x)dx
f g
wherethenmeasurementelementszk(1),...,zkR(n)Q,nominallyindividualpixelintensities,attimetk areexpressed
asavectorzk, g1(i) isthelikelihoodfunctionofelement zk(i) whenatargetexistsatxk,andg0(i) isthatwhena
targetdoesnotexist.
Wenotethattheintegralsintheequationsabovegenerallycannotbecomputedanalytically. Themostrobustquadra-
tureisanSMC,inwhichPDFsareexpressedasanensembleofweightedsamplepoints. Supposethattheaposteriori
statePDFattimetk 1isapproximatedas
�
N
sˆk�1|k�1(x) =Xi=1wk(i�)1�xk(i�)1(x), (8)
wherewk(i)1istheweightandxk(i)1thestateofsamplei,respectively,and�x(i) (x)istheDiracdeltafunctioncentered
atx(i) . W�e’dliketoexpressthe�predictedstatePDFas k�1
k 1
�
N
sˆk|k�1(x) =Xi=1wk(i|)k�1�xk(i|)k�1(x). (9)
Ofthe N samples, B samplesareusedtorepresenttheobjectbirthdensity; i.e.,thefirsttermineq.(3). Ifweindex
thesesamplesasi = N B+1,...,N,
�
psqk�1|k�1 ⇡k|k�1fxk(i|)k�1|xk(i�)1gw(i) i =1,...,N B
wk(i|)k�1 = 8>>>>>>>>< pbq(1k|�k�q1k�1|⇢kk�1fx)k(⇡i|)kk�|k1�|1xfk(ix�)k(1i,|)kz�k1g|xk(ik�)�11g 1, i = N B+�1,...,N (10)
>>>>>>>> qk|k�1 ⇢kfxk(i|)k�1|xk(i�)1,zkg B �
wsahmeprele⇢skufnxdk(ei|)rkg�o1|txhk(ei�)s1a,zmkegdi:synthaemiimcspaosrttahneceressat.mSpilminiglaprlryo,pfoosratlhdeiustprdibautetiostnepo,fxk|k�1,assumingthatthebirthdensity
N
Ik = lk zk|xk(i)k 1 wk(i)k 1, (11)
Xi=1 f | � g | �
andqk k maybecomputedviaeq.(6). TheintegralinthestatePDFupdateis,inessence,anormalizationfactor
|
N
sˆk|k(x) =Xi=1wk(i|)k�xk(i)(x) (12)
n g(j) z(j) x w(i)
w(i) = ( j=1 1 f k | g) k|k�1 , (13)
k|k NQ n g(j) z(j) x w(i)
i=1 j=1 1 k | k k 1
( f g) | �
where N B samples of x(i) are chosen fromPx(1) Q,...x(N) with a probability corresponding to their individual
� k k k 1 k k 1
weights,andtheremainingBsamplesaredrawnfro| m� thebirt|h�density.StandardmethodsinSMCtoavoiddegeneracy
areapplicabletotheSMCimplementationoftheBernoullifilteraswell,suchasregularization.Again,wereferreaders
toashortlistofexistingliteratureonSMCforfurtherdetails[2,13,14].
5. PROCESSINGTELESCOPEDATAUSINGBERNOULLIFILTERS
To test the performance of the Bernoulli filter in processing RSOs, we process three optical tracks of varying SNR
takenduringthenightofOctober21,2011attheNationalCentralUniversityLulinObservatoryinTaiwan[56]. We
shall refer to them as Tracks 1, 2, and 3. Images in the tracks have been preprocessed to remove limb-darkening,
bias, and bright stars, then cropped such that they contain exactly one RSO target. Finally, the intensity values are
normalizedsuchthatthemedianofthemaximumintensityofeachimagecorrespondstoanintensityof1. Thevisual
magnitudes of the targets are 16.70 in Track 1, 17.68 in Track 2, and 18.22 in Track 3. In terms of SNR as defined
through the coe�cient of variation, they are 9.0095, 4.4797, and 2.6983, respectively. Relevant specifications and
parametersofthemeasurementdataarelistedinTable1,andrepresentativeimagesareinFigures1and2. Additional
informationregardingthedatasetcanbefoundin[56].
Thetargetstatewe’dliketoestimatearethekinematicpositionandvelocity[x,y,x˙,y˙]intheimageplane. Weassume
thatRSOsappeartomoveinalinearconstantvelocitymotionwithaconstantacceleration,unbiasedprocessnoiseof
Table1: Relevantinstrumentspecificationsanddataparameters. (Originallyin[56])
Spec Value
Detectortype CCD(whitelight)
(#ofImages)/track 29
Size 301px 301px
⇥
Resolution 3arcsec/px
FOV 0.25deg
Exposuretime 5.9sec
Readouttime 2.66sec
Mounttracking O↵(Earth-fixed)
Figure1: Thefirstimageframefromeachtrack. Fromlefttoright,Track1,2,and3.
Figure2: RSOscroppedoutfromTracks1(toprow),2(middlerow),and3(bottomrow). Timeelapsesfromleftto
right. (Originallyin[56])
0.1px/s2 standarddeviation. Weassumethatthereisnoaprioriinformationregardingtheinitialstateofthetarget
noritsexistenceprobability;thus,wesettheinitialstatePDFtobeuniformwithasupportof[0,300]pxinposition
and [ 2,2] px/sec in velocity. We similarly assume no a priori information for the target birth process. The target
�
birthprobabilityissetto0.05,andthesurvivalprobabilityto1.
Webaseourmeasurementlikelihoodfunctionontheoneproposedin[58]
g1(i) zk(i)|xk =N zk(i);hk(i)(xk),�02 (14)
f g f g
g(i) z(i) = z(i);µ ,�2 , (15)
0 k N k 0 0
f g f g
where (x;µ,�2) isanormaldistributionoverstate x withmean µandvariance�2, µ isthenoisebias,and� is
0 0
thenoiNsestandarddeviation. h(i) isthepointspreadfunction
k
hk(i)(xk) = 2⇡I�k 2exp2�(ix �x)22�+2(iy �y)23, (16)
h 6 h 7
6 7
6 7
where Ik is the source intensity, �h is the blurring facto64r, and (ix,iy) are the i75mage plane coordinates of pixel (i).
Ik determines the sensitivity of the conversion from pixel intensity values in the vicinity of a particle to the particle
weight. Thelikelihoodratiois
lk(zk x) =exp (µ0�hk(i)(xk))(µ0�2zk +hk(i)(xk)) . (17)
| 2 2�2 3
66 0 77
6 7
6 7
Although functions better suited for raw CC6D intensity output certainly exist – e.g7., those with a strictly positive
4 5
support – the one proposed is fast to compute. Further cautionary notes regarding the likelihood function are listed
here.
Parameterselection Nominally,onewouldchooseparameters�0, µ0, �h,and Ik baseduponstatisticsofthemea-
surementdata. Suchvalues,however,makeeq.(4)muchlargerthancanbehandledevenwithextendedpreci-
sionarithmetic. Thus,inthispaper,weestimateonly µ ,whichareenumeratedinTable2,bytakingthemean
0
intensityovereachimage. Therestaretreatedastuningparameters: �h = 1foralltracks,whereas�0 and Ik
aresetasinTable2. ToelevatetheBernoullifiltertospacesurveillanceoperations,methodstoautomatically
choosetheseparametersshouldbestudied.
Clutter ItisapparentfromFigures1and2thatnotallofthesignalinthetracksoriginatefromRSOs. Thelikelihood
function, then, mustaccountforspuriousinputssuchasbrightstarsthatareonlypartiallyremoved, dimstars
thatarenotremoved,andfalsedetectionsduetocosmicrays. AsthetargetRSOsandexposuretimeschosenin
thecurrentdatasetaresuchthatRSOsdonotappearasstreaks,weamendeq.(14)suchthat
g1(i) fzk(i)|xkg = 8>>< Ng0(i)fzk(zik()i;)hk(.i)(xk),�02g ((optth<erwpti,smei)nandpk < pk,min) (18)
>> f g
For each particle, a Student’s t te:st for normality is conducted on the intensity values for the fourth to ninth
closestpixelsforallhorizontalscanlinesintheimage. Ifthe p-valueislessthansomeuser-definedminimum
pt,min,thenthesignalisdeemedtobesospreadoutthatitisduetoastreakingstarratherthananRSO.Similarly,
anormalityhypothesistestbasedonkurtosisisconductedontheintensityvaluesforthenineclosestpixelsfor
allhorizontalscanlinesintheimage[6]. Ifthep-valueislessthansomeuser-definedminimumpk,min,thenthe
signalisdeemedtobesoleptokurtic,i.e.,“peaked,”thatitisduetoacosmicrayratherthananRSO. pk,min is
setto5 10�5 foralltracks,whereas pt,min issetasinTable2. Thehigher pt,min isset,themoreaggressively
⇥
thefilterwillignoresignalthatarespreadout.
Muchmoreworkisrequiredtodeviseastatisticallyrigorousyetcomputationallytractablelikelihoodfunctionsuitable
forRSOsinallorbitalregimes, andthuswillbethefocusoffuturedevelopment. Weareencouragedbytheresults
in[40],whichsuggestthatevenasmallamountofaprioriinformationcanstronglyconstrainboththesearchvolume
andshapeofRSOimagery.
Table2: Filterparameterschosenforeachtrack.
Track1 Track2 Track3
µ (Estimated) 0.079066 0.10309 0.10722
0
Ik 0.375 0.50309 0.60722
pt,min 0.01 0.1 0.2
N 10,000 20,000 20,000
TheparticularflavorofSMCimplementedforthispaperisapost-regularizedparticlefilterusingtheoptimalband-
width for a Gaussian kernel [13]. All computation was conducted in GNU Octave on a quad-core Intel Core i5
(Haswell,3.40GHz)workstationwith16GBofRAM.ThenumberofparticlesintheSMCwaslargelyconstrained
byRAM;Octavebecameunstableifmorethan20,000particleswereprocessed. Thenumberofparticles N usedto
processeachtrackisgiveninTable2. Computationtimewasapproximately2frames/minuteforevery104 particles
intheSMC.
6. RESULTS
Figures3–6areresultsfromsuccessfulBernoullifilterruns. Foreachtrack,Figures3–5showtheevolutionofthe
SMCparticlecloudalongwiththeestimatedprobabilityoftargetexistence,andFigure6containthetimehistoriesof
theestimatedmeanaccuracyandcovariance.NotethatforTrack3,thefilterdidnotconvergereliablywiththenumber
of particles that fit in memory. A “hot start” was thus implemented; the support of the initial state PDF in position
spacewasreducedtoa100px 100pxregionaboutthetrueposition.
⇥
Allrunsshownconvergetowithin1pxofthetruepositionand0.2px/frameinvelocityinboththexandydirections,
whichcorrespondto3arcsecand0.07arcsec/secerrorinangle/angle-rate. Theestimatedcovariancesareconsistent
withtheactualestimationerrors.Asnotedpreviously,thesestatisticsmaydirectlybehandeddowntheSDAworkflow,
beitinitialorbitdeterminationorconjunctionassessmentalgorithms. Additionally,theproposedlikelihoodfunction
ise↵ectiveindistributingSMCparticlesonlyinthevicinityofsignalfromRSOs.
One“lessonlearned”thatemergedfromtheseresultsisthateveninsolvingfortheprojectedmotionoftheRSOinthe
imageplane,thenumberofparticlesrequiredforthefiltertoconvergebecameratherlargeforlowSNRobservations.
Since we have ignored shape parameters of targets in this work, all of the instantaneous information from images
areinthepositionspace. Velocityspaceinformationarisessolelyfromthedynamics; particlesareheavilyweighed
over multiple images if their velocity is close to truth. Therefore, the particle cloud after the update step must have
enoughinstantiationsinthevelocityspaceforanobjecttobeconsistentlytracked. Therearemultiplepathswecan
taketoalleviatethisissue. Forexample,wearemotivatedtoconductfutureworkinacompiledlanguagewithmulti-
coreorGPGPUsupporttodrasticallyincreasethenumberofparticlesintheSMCwhileretainingfastcomputational
turnaround [49]. Alternatively, we can improve the implementation of the Bernoulli particle filter; e.g., add labeled
RFSs[23,59]orchangetheregularizationtechnique[13].
7. CONCLUSIONS
In this paper, a Bernoulli particle filter that directly processes sensor output from electro-optical systems is imple-
mented so as to simultaneously detect and track resident space objects. Compared to current operational practices
in ground-based optical space surveillance, the proposed method has the potential to reduce computational burden
for high-sensitivity fast-cadence imagery, provide consistent uncertainty quantification from the time of observation
through the delivery of actionable information, and lower the signal-to-noise ratio required to maintain object cus-
tody. Approachesrootedinfinitesetstatistics,liketheBernoullifilter,areadvantageousintheirabilitytoexplicitly
accountformultipletargetsorthelackthereof. Intheimplementation,filterparameterselectionandclutterrejection
through the measurement likelihood function were identified as having the largest influence on filter performance.
Figure3: EvolutionofparticlecloudforTrack1inred. Timeflowsfromlefttoright,thentoptobottom. Contourof
pixelintensityingrayscale. Truepositionstateindicatedbycyanlines. Inthebottomright,theestimatedprobability
oftargetexistenceasafunctionoftime.
Figure4: EvolutionofparticlecloudforTrack2inred. Timeflowsfromlefttoright,thentoptobottom. Contourof
pixelintensityingrayscale. Truepositionstateindicatedbycyanlines. Inthebottomright,theestimatedprobability
oftargetexistenceasafunctionoftime.
Figure5: EvolutionofparticlecloudforTrack3inred. Timeflowsfromlefttoright,thentoptobottom. Contourof
pixelintensityingrayscale. Truepositionstateindicatedbycyanlines. Inthebottomright,theestimatedprobability
oftargetexistenceasafunctionoftime.
2 x 2 x 4 x
1 y 1 y 2 y
[px]0 [px]0 [px]0
-1 -1 -2
-2 -2 -4
5 10 15 20 25 5 10 15 20 25 5 10 15 20 25
1 2
0.4
me]0.2 me]0.5 me]1
[px/fra--00..042 [px/fra-0.05 [px/fra-01
-1 -2
5 10 15 20 25 5 10 15 20 25 5 10 15 20 25
Image # Image # Image #
Figure 6: Dots indicate the error in the estimated mean position (top) and velocity (bottom) as a function of time.
Linesarethe 3-�boundsasderivedfromtheestimatedcovariance.
±