• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Modeling Impacts on Space Situational Awareness PHD Filter Tracking

    2016-12-12 08:52:41Frueh

    C.Frueh

    Modeling Impacts on Space Situational Awareness PHD Filter Tracking

    C.Frueh1

    In recent years,probabilistic tracking methods have been becoming increasingly popular for solving the multi-target tracking problem in Space Situational Awareness(SSA).Bayesian frameworks have been used to describe the objects’of interest states and cardinality as point processes.The inputs of the Bayesian framework filters are a probabilistic description of the scene at hand,the probability of clutter during the observation,the probability of detection of the objects,the probability of object survival and birth rates,and in the state update,the measurement uncertainty and process noise for the propagation.However,in the filter derivation,the assumptions of Poisson distributions of the object prior and the clutter model are made.Extracting the first-order moments of the full Bayesian framework leads to a so-called Probability Hypothesis Density(PHD)filter.The first moment extraction of the PHD filter process is extremely sensitive to both the input parameters and the measurements.The specifics of the SSA problem and its probabilistic description are illustrated in this paper and compared to the assumptions that the PHD filter is based on.As an example,this paper shows the response of a Cardinality only PHD filter(only the number of objects is estimated,not their corresponding states)to different input parameterizations.The very simple Cardinality only PHD filter is chosen in order to clearly show the sole effects of the model mismatch that might be blurred with state estimation effects,such as non-linearity in the dynamical model,in a full PHD filter implementation.The simulated multi-target tracking scenario entails the observation of attitude stable and unstable geostationary objects.

    Space Situational Awareness,Space Debris,Finite set statistics,Tracking,Modeling.

    1 Introduction

    The aim of multi-target tracking is to estimate the number and states of objects in a surveillance scene,when only non-perfect sensor measurements are available.In the multi-target tracking regime,in general,two different ways to explore this paradigm have been followed so far:track-based approaches and population-based approaches.The track-based approaches are intuitive,as they keep information on each of the different targets[Bar-Shalom(1978);Reid(1979);Blackman and Popoli(1999)].The measurements are associated explicitly to the different targets.However,the track-based approaches rely on expert knowledge for track acceptance and track rejection and cannot provide a description of the scene in the absence of new measurements.Multi-hypothesis filters(MHT)and joint probabilistic data association(JPDA)techniques are popular representations of the track-based approaches[Blair and Bar-Shalom(2000)].

    Finite set statistics,on the other hand,seeks to represent the population as a whole in a single random object[Mahler(2007)].As a result,new observations influence the whole population instead of just a single object track.This has the advantage that a probabilistic description of the complete scenario is available and regional statistics can be readily extracted without making heuristic assumptions.However,as there is no explicit track association,no track history is available,which makes the results of the filter less intuitive.In order to make the method computationally tractable,probability hypothesis density(PHD)and cardinality probability hypothesis density filter(CPHD),or Multi-Bernoulli filters became increasingly popular.In the context of space situational awareness(SSA),the objects of interest are Earth orbiting satellites and uncontrolled objects.Observations are mainly ground based;only non-resolved images are available for the vast majority of objects.The US Strategic Command currently tracks around 18,000 objects[US Air Force(2016a,b)];however,statistical measurements suggest a by orders of magnitude larger number of objects of interest in the near Earth space[Klinkrad and Sdunnus(1997)].The tracking of those objects is crucial to ensure a sustainable use of space and to protect the active space assets.Both track-based and population-based approaches have been explored in the solution of SSA tracking[Frueh and Schildknecht(2012);Frueh(2011);Kelecy,Jah,and DeMars(2012);Kreucher,Kastella,and Hero III(2004);DeMars,Hussein,Frueh,Jah,and Erwin(2015);Jia,Pham,Blasch,Shen,Wang,and Chen(2014);Gehly,Jones,and Axelrad(2014);Jones and Vo(2015);Cheng,DeMars,Frueh,and Jah(2013);Cheng,Frueh,and DeMars(2013);Singh,Poore,Sheaff,Aristoff,and Jah(2013)].

    In the current publication,we want to focus on a different aspect of the PHD filter.First,the specifics of the SSA tracking scenario are illuminated.Secondly,the well-known filter equations,as they are normally used,are recapitulated and the specific assumptions that are intrinsically made are made explicit.Their applicability in the context of SSA tracking and the effect of the mismatch between thefilter modeling assumptions and the SSA dynamical realities are investigated.Simulations are created for a cardinality only version of the PHD filter that estimates only the number of targets.A ground based optical observation scenario of geostationary controlled and uncontrolled objects has been chosen.The cardinality only PHD filter is utilized,because the effects of the modeling parameters can be shown more easily and allow for the determination of the effects independently from the particular state propagation and initial orbit determination method.The model assumptions that are made in the derivation are identical to the ones of the full PHDfilter.

    2 Tracking in the SSA regime

    In the SSA regime,as in other multi-target tracking problems,tracking has to solve three different sub-problems.First,how many objects are in the scenario of interest,the cardinality problem;second,the state of each object at each time,the state estimation problem;and third,the data association problem,which is the determination of which measurement belongs to which object.As mentioned before,the data association is not explicitly solved in population-based tracking.In the SSA regime,the situation arises that the movement of the objects is fully deterministic in principle.However,the precise knowledge of all the relevant parameters influencing especially the non-conservative forces are often not at all or only incompletely available.Such parameters could be the precise object shape,material properties,space weather parameters,etc.The active objects may maneuver.Both the maneuvers(even single thrust ones)and the effects and influences of not precisely known forces and torques,albeit uncertain,lead to a slow alternation of the orbit.In that sense,the problem may seem easier than other multi-target tracking scenarios.

    However,there are several specifics of the SSA multi-target tracking problem,which make it significantly harder than many other multi-target tracking problems.The measurement process and schedules and the lack of precise knowledge that is relevant for the dynamics pose severe challenges.In general,the exact knowledge is lacking that would allow to reliably predict the objects with small uncertainties over longer periods of time.The thrust information and object details for highfidelity non-conservative object-dependent force models are not available.The observation scene of interest,in the terminology of usual multi-target tracking,is the whole area from the Earth surface(for the decaying objects)to beyond geostationary orbit(42,000 km,resp.36,000 kilometers from the Earth’s surface),whereas the field of view(FOV)of the sensor is normally a few square degrees.Given the large populationand the limitedsensing capabilities,there are necessarily long time intervals of 30 days or more with no observation update on single objects[Flohrer,Schildknecht, R.and St?veken(2005);Musci(2006);Musci,Schildknecht,and Ploner(2004);Musci, Schildknecht, Ploner,and Beutler(2005);Oswald,Krag,Wegener,and Bischof(2004); Potter(1995)].The quantity of interest is classically the object state,i.e.position and velocity.For precise modeling,theoretically also attitude and object characterization parameters should be part of the state.Even for the classical state definition,only part of the state is observable.In optical observations with a single image,only two angles are measured.The situation is only slightly better in optical observations with long exposure images or when the information of two subsequent images are combined;then also angular rates are available in addition to the astrometric position.Radar measurements offer the slant range,and with larger uncertainties,the two angles.A Doppler radar also allows to measure the range rate.Laser range measurements offer very precise range measurements and angles.Prior to being able to propagate the state,either a classical first orbit determination has to be performed,which presupposes that more than one measurement can be associated to the same object,or an admissible regions approach has to be utilized.A further challenge is posed by the fact that only non-resolved object images are available.This means that no additional information,apart from the observed part of the state,is available for the data association scheme.The non-linear measurement process itself is imprecise;hence,the measured parts of the non-resolved object state bears uncertainties.For more details on the uncertainties and probability of detection of optical sensors see[Frueh and Jah(2013);Sanson and Frueh(2015,2016)].Fig.1 shows an optical object image and its spread over the pixel grid.

    Figure 1:Optical image of geostationary object.

    Figure 2:(a)and(b)examples of true object images,(c)example of cosmic ray trace(ZIMLAT,University of Bern AIUB).

    Furthermore,not all detections are true detections.Clutter contaminates the measurements.In optical measurements,e.g.cosmic ray events contaminate the measurements;they cannot readily be discriminated from true measurements,as can be seen in Fig.2.Because the objects move relative to the star background and their relative motion is a priori not known,simple stacking techniques to filter out cosmic ray events completely,cannot be used.Alternatively,cosmic filters can be applied to at least reduce the number of clutter measurements[Frueh and Schildknecht(2012)].

    Figure 3:(a)five subsequent images of the same object,spaced by 30 seconds,(b)fi ve geostationary objects in front of star background(ZIMLAT,University of Bern AIUB).

    Not every time that the object is actually in the field of view it is also detected.The signal can be below the noise level or e.g.occultation occurs.The attitude motion of objects can be very rapid,leading to large brightness changes in short time intervals;one example of this is shown in Fig.3a.Occultation in optical observations occur when the object is in front of a star,as illustrated in Fig.3b,which shows five geostationary objects in front of the star background[Frueh(2011)].This illustrates how easily occultation can occur in regular observations.

    3 PHD filter and Cardinality only PHD filter

    Figure 4:Derivation and data flow steps of the PHD filter.

    Instructive derivations of the PHD filter can be found in [Mahler (2007); Houssineau,Delande,and Clark(2013)].Gaussian mixture closed-form PHD filter derivations can be found in Clark,Panta,and Vo(2006);Vo and Ma(2006).Closely following the derivation in[Houssineau,Delande,and Clark(2013)],the formulation of the probability generating functionals(PGFls)Fare used to describe the multi-object population and then differentiated in order to find the first-order moments.

    A PGFl is a mapping of an arbitrary functionh:X→[0,1],generated by a series of probabilitiespΦ,associated to a point process Φ:

    wherexiare the different realizations of the point process in the spaceX,in our case the state space and the states of the different possible object realizations.A state is defined as position and velocity of the object.E[]is denoting the expectation value.It can easily be seen thatF(0)=pΦ({0}),which is the pdf in the absence of any objects(empty set),andF(1)=1.Fig.4 shows the differentiation and data flow steps of the PHD filter.

    For the prediction step,the PGFl has the following form:

    whereFt?1(h)is the PGFl describing the population at the previous time step,Fcrea(h)is the PGFl describing the object creation process,psuris the probability of survival,T(xt‖xt?1)is the state transition from the state att?1,denoted as xt?1,to the timet,denoted as xt.x is the set of the individual states of all objects.Eq.2 is derived as the superposition of the surviving objectsFsur(h)and the newly created onesFcrea(h).xt={x1,...,xk}tdescribes the set of the states of the objects at a given time t.It has to be noted that the PGFl is defined on the single object state.

    In order to arrive at the final simple form shown in Eq.2,it is assumed that the targets are independent(A1)and that the newborn targets and the surviving ones are independent (A2); the probability of survivalpsuris assumed to be the parameter in a Bernoulli process(A3).A Bernoulli process ? can be defined as follows.For a spatial probability density functionf(x)(Rf(x)dx=1):

    for a probability parameterp,which can take values between zero and one.φis the empty set.In the concrete case,this means the target either survives with a probabilitypsur,with a state distributed according tof,or dies with a probability 1?psur.That completely describes the survival-dying process statistically.

    The PGFl of the measurement process and the objects can be expressed as:

    with the arbitrary functionsh,g:X→[0,1].Gclutis the PGFl of the clutter process,pdis the probability of detection,Lis the likelihood function,and z is the measurement.In the derivation of Eq.5 it is assumed that the measurements are produced independently(A4),the probability of detectionpdcan be modeled a Bernoulli parameter(A5),and that the clutter process producesmeasurementsindependently from the actual objects in the scene(A6).

    Utilizing Bayes’rule and the fact that the actual number of measurementsZ=z1,z2,.....,zn(clutter and real object measurements) the sensor has received is known,the PGFl of the data update,derived as the directional derivative of the PGFl in E-q.5,has the following form:

    In order to derive the first moments,it is utilized that the first order moment density is the first-order directional derivative of the corresponding PGFl:

    This leads to the following filter equations for the first-order moment densities:

    whereμcrearepresents the first moment of the object creation.In order to derive these closed-form expressions,one has to assume that the predicted targets is described by a Poisson point process(A7)and that the point process describing the number of clutter measurements can be modeled as a Poisson process,too(A8).

    A simplified form of the PHD filter is the cardinality only version of the PHD filter,which only solves the cardinality problem(not to be confused with the Cardinality PHD filter,CPHD,which also,in addition to the first-order density,propagates the cardinality distribution of the random finite set.The filter equations for the cardinality only PHD filter are the same as in Eq.9,except for the state estimation part.

    wheremrepresents the number of measurements.

    4 Discussion of the Premises and Assumptions

    4.1 Independence of Objects

    In the regime of SSA,the assumption(A1)that the objects are independent is most of the time fulfilled.Independence means that the dynamics of the objects are not coupled;that is,one object does not influence the dynamics of the other.The single target differential equation can be written as:

    where x is the geocentric position of the object,G the gravitational constant,M the Earth mass and V(x)the Earth gravitational potential,third body gravitational perturbations of the Sun and Moon(k=1,2)with the states xk.Finally,∑a is the sum over all non-gravitational accelerations acting on the satellite.Popular nonconservative perturbations include solar radiation pressure and drag.An object interdependence could only be introduced via those accelerations;realistically,this could be the case for decaying objects and drag wake effects.A further case is on-orbit collisions at time scales where the dynamic effects from the collisions are still dominating.

    4.2 Independent Newborn and Surviving Objects

    The assumption(A2)that the surviving and newborn objects are independent is an assumption that might be in conflict with the real dynamical situation in the near Earth space.For sure,there is a weak coupling between the surviving targets and new launches,as new launches in part tend to be motivated to refurbish old assets.However,the dead objects still remain in orbit for extended periods of time.More crucial is the object creation that is not caused by new launches but via satellite aging and on-orbit collisions.Following Kessler and Cour-Palais(1978),the collision rate between satellites can be expressed as:

    where σ is the spatial density created by the other objects in the volume element dV,ˉv is the average relative velocity,andˉA is the average cross sectional area.Although Kessler’s theory might be controversial in some of his conclusions,the source of debris creation from collisions and impact motivated break-up is valid.Hence,there is a strong and direct coupling between the surviving targets.As with the assumption(A1),a violation can be found in times of tracking through collision events.

    4.3 Probability of Survival

    Assumption(A3)says that the object either survives with a well-defined probability value,or it dies;probability of survival is represented as the parameter of a(constant)Bernoulli process.Often,in the use of a PHD filter,this quantity is defined as a constant parameter;however,there is no restriction to make it time varying,which is more appropriate for the SSA problem.Object death is only immediately relevant in the strictest sense for objects in decaying orbits;however,sensing limits should not be neglected.In the absence of collisions,the model for the probability of survival can be modeled as the following in the absence of collisions:

    If the perigeerper=(1?e)a(eccentricitye,semi-major axisa)is higher than the upper bound sensing limit,of one or a combined sensor suite,ρhi(SNR(t)),or the apogeerap=(1+e)ais lower than the lower bound sensing limitρlo(SNR(t)),the object has left the sensing scene.No new measurements on the object can safely be expected;it hence is dropped.A precaution against orbital uncertainties can be made with an offsetε.A good choice ofεthen could be defined using the three sigma bound in eccentricity?eand the three sigma bound in semi-major axis?aasε= ?a+2?a·e.The lower sensing limit might coincide with the Earth’s surface.One could imagine the situation in which a target is sought to be keptalivealthough no new measurements can be expected.In that case,a maximum time span since it has left the scene(apogee,perigee condition)can be put as an additional constraint.It has to be noted that the detection limit depends on the signal-to-noise ratio(SNR).In a conservative approach,the maximum signal to noise ratio(SNRmax),in case it is known,should be used,which occurs for minimum phase angle and zenith observations under best conditions.More details on the SNR computation can be found in section 4.5.The probability of survival model as proposed here is hence linked to the probability of detection and can be written as:

    Theoretically,one could also take the specific sensor scheduling scenario into account,which might eliminate further objects;however,it would mean that the sensing schedule can be predicted for all times and excludes detections of objects under potentially growing?aand?e.

    The problem of object death in its current form in collisions that destroy the parent object is an instantaneous process that is hard to predict given the state uncertainties in the state propagation.Even when a collision seems certain,it is almost impossible to predict if the parent object remains sufficiently intact to allow for a survival in the classical sense,especially as object sizes and potential overlap in the collision are hardly known.However,as stated above,assumption A1 and A2 would be already violated in the case of a collision.As fatal on-orbit collisions are(fortunately)still very seldom events,with a combined reported and suspected number of five in the past sixty years of space fairing,see Klinkrad,Wegener,Wiedemann,Bendisch,and Krag(2006)and references therein,and best measures are taken to prevent further collisions,probability of survival even when taking collisions into account is still practically one for non-decaying objects within the sensing realm of the sensor suite.This means,even when processes like spawning are integrated in the filter process,it is extremely challenging to incorporate it without a specific trigger or external knowledge,because of the evanescent probability.

    4.4 Independence of Measurements

    Assumption(A4)requires independence of the measurements.Here,one needs to discriminate between optical measurements and radar measurements.In optical measurements,the measurement is the astrometric position(two angles,such as e.g.right ascension α and declination δ ).A potential coupling of the measurements of different objects that are in the same observation frame is given via the common noise level and the readout process.However,as the background noises are dominated by other celestial sources,such as zodiac light and the background stars and atmospheric scattering,this effect is negligible.Extracting the velocity information from two subsequent frames introduces a huge dependence between the two subsequent measurements:

    which is often done for two measurements of the same object that are very close in time together(difference only εt,a small time step relative to the orbital period of the object in question).The dependence can be circumvented in counting the two adjacent measurements only as a single angle and angular rate measurement(α(t1),δ(t1),˙α(t1),˙δ(t1))and not reusing the astrometric position measurement at time t2.Alternatively,angular rates may also be extracted from single exposures.The situation is different for radar measurements.Radar detection methods require coherent pulse to pulse integration,computing the match-function of a set of range and velocity values,see[Markkanen J.(2006);Krag,Klinkrad,Jehn,L.Leushacke,and Markkanen(2007)].Unlike in the optical case,those are normally treated as a single measurement,hence leaving the dependence of the measurements in the data,when they are interpreted the traditional way.

    4.5 Probability of Detection

    Assumption(A5)states that the probability of detection is modeled as the parameter of a Bernoulli process,hence that the object is either detected with a given probability or not.The probability of detection is directly dependent on the signal to noise ratio(SNR)[Frueh and Jah(2013);Sanson and Frueh(2015,2016)]and time-varying.Albeit it does not contradict the derivation of the PHD filter directly,the probability of detection is often modeled constant,which does not match the SSA situation.Moreover,very low probabilities of detection can occur,which,by construction,the PHD filter has difficulties to deal with.

    In order to model the probability of detection correctly,all light that is received at the detector has to be known,and how the detector converts it,and how the signal evaluation process is performed.Albeit,those steps are simple,they are often overlooked.The following discussion can in parts be found in[Sanson and Frueh(2016);Frueh and Jah(2013);Sanson and Frueh(2015)].The signal from the object is dominated by the observation geometry at the time of the observation,the object geometry,and surface properties.In the derivation here,the focus is placed on indirect illumination by the sun;however,the results are easily adapted to active illumination,by shifting the place and strength of the source accordingly.

    4.5.1Object Irradiation

    The reflected irradiationIat the location of the observer in[W/m2],is:

    ISunis the irradiation of the sun at the location of the object for a given wavelengthλ.In general,ISun(ˉλ)·?λmay be approximated with the solar constant,I0=1367.7W/m2;to be more precise,the deviation from the one astronomical unit scaling should be taken into account.ˉλwould be the mean wavelength.The mean should be chosen as the mean of the sensing spectrum of the sensor at hand.Ais the complete area off which the light is reflected,rtopois the topocentric distance(distance object to observer),and Ψ is the phase function.In the approximation,the precise solar spectrum and the spectral response of the object has been neglected and a combined white light response is modeled.For filter measurements,or color CCD sensors,the correct spectral response has to be taken into account.The phase function depends on the object properties,the shape and reflection properties;it defines how the light is reflected off the object for a given incident angle of the light and angle of observation.Ψ can be defined as:

    μ0is the incident light direction,μis the outgoing observer direction,andαis the phase angle between the incoming flux and the direction of the observer.In general,the bidirectional reflection model of choice can be applied.In the context of this paper,the surfaces are modeled as a mixture of Lambertian,diffuse reflection and absorption is adapted.For a Lambertian reflectionFor some simple shapes,closed-form solutions for the phase function can be obtained,such as the sphere;without loss of generality,the observer can be placed in the xy plane of the sphere;definingμ0=sinθcosφandμ=sinθcos(φ?α)leads to the well-known result,see Wertz(1978):

    for a diffuse reflection parameterCd(λ).αwould be the classical phase angle(angle between observer,object,sun).In general,real objects are modeled for a mixture between e.g.specular,Lambertian reflection,and absorption.In general,Cd+Cs+Ca=1 for opaque materials.Specular reflection on a sphere is not taken into account here as classical glints cannot be caught from ideally spherical space objects.Glints are only possible from flat surfaces.Spherical surfaces always return light to the observer as long as they are not in the Earth shadow or the illumination source is exactly opposite to the observer.For a flat plate,the integrals over the hemisphere have to be replaced by the integral over the illuminated surface,and are hence very simple.μ0=cosβandμ=cosγfor the flat plate case,whereβis the angle between the direction to the illumination source,S,and the normal vector of the surface,N.γis the angle between the observer O and the normal vector:

    4.5.2Celestial Background Irradiation

    Besides the irradiation of the object,other background sources also enter the detector.The irradiation of the object is not the only light that is reflected towards the observing sensor.In optical observations,several background sources,such as zodiac light,direct and stray moonlight,atmospheric scattering,and atmospheric glow,may be taken into account.Intensive studies on the background sources and tabulated values can be found in[Daniels(1977)Allen(1973)];application to SSA observations can be found in Frueh and Jah(2013).

    For celestial background sources Eq.17 has to be adapted as the following:

    Ii(λ)denotes the wavelength dependent exact irradiation values of the different sources.The signal expression in the following equations are derived using the approximation in Eq.22 at the mean optical wavelength of 600nm.

    Airglow spectral radiation,IAG(λ),is the faint glow of the atmosphere itself,which is caused by chemiluminescent reactions occurring between 80 and 100 km.The irradiation values can be directly computed from the irradiation valuesJi,as shown in the subsequent equations.Irradiation values are well documented in literature,e.g.in[Daniels(1977)Allen(1973)].Atmospherically scattered light is the sum of all light that is scattered by the atmosphere.Both can be modeled as the following:

    wheresis the angle under consideration,in case of the telescope.It can be chosen to be the field of view(FOV)of the sensor,a field averaged background would then be derived.This is legitimate for small to medium sized FOVs;for wide field telescopes,this approach is not recommended.The most precise modeling is achieved whensis chosen to be the angle that is integrated into a single pixel.

    Diffuse galactic light is scattered sunlight by small particles concentrated in the ecliptic.It can be modeled as the following:

    whereJGALis the spectral irradiation at zero galactic latitude(β).Zodiac light is the diffuse stray light from the sun,which is visible in the night sky even long after

    γ,δare the longitude and latitude in ecliptic coordinates andJZODI(γ,δ)is the total irradiation per unit angle.As zodiac light can be by far the strongest background source and has a large variability,the use of the exact tabular values for the total irradiation is recommended.

    Another major light source is the faint stars in the background.One way to include stars is to include them at the exact position as they appear in extensive star catalogs.However,this is a very time consuming procedure.Tables exist with the number of stars of given photographic magnitudes,see[Frueh and Jah(2013);Allen(1973)].Using these,they can be converted to radiance values,using the known spectral distribution of faint stars,as in the case of the galaxies.The conversion is done in the blue wavelength(440nm)to have the best equivalence with the photographic magnitudes denoted bym.This leads to the spectral star irradiation:

    wherenis the number of stars in the assigned bin. The irradiation values correspond to the irradiation without an atmosphere.Sometimes,one is not interested in a specific wavelength,but the total radiation;in this case,one can integrate or use approximations for the white light.

    4.5.3 Sensor response

    The response to any of the incoming celestial irradiation,from the object of interest or from any of the celestial sources,can be modeled as the following:

    with the speed of light,c,Planck’s constanth,the area of the aperture,D,and the obstruction of the aperture,d.R(ζ)is the van Rhijn factor,which can be approximatedto first order and describes the deviation from the zenith by angleζand the additional air mass and thickness that has to be accounted for in low elevations[Daniels(1977)].A very simple atmospheric extinction model has been adapted,withτbeing the atmospheric extinction coefficient.Note that atmospheric extinction is zero for the airglow and atmospheric scattering signals,because they are created in the atmosphere rather than being extinct via passing through it.

    The count rate is derived from the signal via the time integration,dt,during which the sensor is able to catch photons:

    whereQis the quantum efficiency.The photon count is a statistical process.A model for the signal and subsequent count at the detector is a Poisson distribution.The approximation neglects the shutter function itself and assumes that the integration time is the same over all the field of view of the sensor.The signal,however,is spread over several pixels,which report their count rates independently.The count rateCis hence interpreted as the mean of a Poisson process.Strictly speaking,quantization already happens at the signal level.However,this is without consequences,as only in the count,the time averaging takes place.The count rates are then transformed via the so-called gain,g,from the count level to the analog to digital units(ADU)that are read out electronically.The gain is not completely linear over all counts of the sensor,although a linearity of the gain is desired.Non-linearity normally occurs when reaching the saturation level of the sensor.In the process of the transformation,a truncation takes place as normally only integer ADU values can be reported.

    For simplicity,it is assumed that the complete irradiation is counted in a single pixel.Depending on the pixel size in combination with the setup of the optics,the Airy disk diffraction pattern that is formed is counted in several adjacent pixels.In this case,normally only the pixel counts of the first maximum are taken into account,which is around 84%of the total count of the complete diffraction pattern.This is especially true if the object is moving relative to the sensor.

    The sensor itself introduces additional counts,so-called dark noise,that stem from electron fluctuations that occur even when the shutter remains closed.The readout process itself creates readout noise.

    4.5.4Probability of Detection

    The probability of detectionPDcan be defined as the following(for further details,please also see Sanson and Frueh(2015,2016)):

    where SOBJis the object signal,modeled as the mean of a Poisson distribution,and N is the cumulative noise.Again,it is assumed that the signal is integrated in a single pixel.Following the derivation of[Merline and Howell(1995)],the noise can be represented as:

    where SSKYis the variance of the celestial background sources,modeled as a Poisson distribution with the count rate as the parameterCSKY.SDARKis the mean of the dark noise,modeled as a Poisson distribution with the parameter CDARK,is the readout noise,which by the central limit theorem is modeled as Gaussian distributed with the variance of a Gaussian distribution,and SUis the truncation noise variance,which approximately can be modeled as uniformly distributed,rendering the variance CU=g2/12,where g is the gain.nBis the number of pixels that have been used in the background determination.For the derivation of the CCD equation,including improved modeling for cases,when the source is not integrated in a single pixel,see Sanson and Frueh(2015,2016).Inserting it into the probability of detection formula:

    which allows to evaluate the different distributions in a cumulative density function for the Poisson distributed parts and conditional probability density function on the remaining terms:

    4.6 Clutter

    Assumptions(A6)and(A8)focus on the clutter process.The assumption that clutter and measurements are produced independently is easily fulfilled,as the physical processes that are generating the object signal and the clutter signal,such as the cosmic ray events in the optical observation,are independent.The signal of the object is the irradiation from the object,either stemming from the sun or from the observation station and is radiated back to the observer.Clutter,such as cosmic ray events,are charged particles that accidentally impinge on the sensor.No part of the measurement process excites or influences the charged particles that arrive at the detector.Cosmic ray events are a textbook example of a Poisson process that describes the arrival of the charged particles at no distinct order at this time integrated process.

    4.7 Predicted Object Distribution

    Assumption(A7),that the predicted object distribution is a Poisson process,is one of the most crucial ones.As we have seen in Eq.12,the underlying dynamics is of course purely deterministic,even though sometimes,exact knowledge on all parameters is missing.Together with the way objects are created or die,single outbursts in the population in crucial collisions,and a few launches and decays,the object population cannot readily be suspected to be in a close resemblance of a Poisson process.In the previous section,no weighting is applied in the filter equations,which means the filter is highly responsive at each new time step.A Poisson process implies that the actual number of targets can jump,same as it accurately applies to the clutter population that is detected.Actual jumps in the factual populations only happen in collisions,which are rare and problematic as discussed above already.The actual population is concentrated in certain orbits and not highly variable,hence the predicted object population on the observation frames is highly dependent on the specific observation scenario.In actively tracking specific objects,often called follow-up tracking,in contrast to so-called blind tracking or surveying,it is expected that the objects that are visible in the observation frames is more or less constant.

    It can be seen that the filter is put to its edges already with the birth and death processes,which takes on extreme values and does not settle on the middle ground.This is in connection to the Poisson process assumption for the object population,which might be highly inadequate.

    5 Simulation Results

    Simple tests to show the impact of the modeling parameters on the filter performance have been run.Five different variations of the same scenario have been created.For the tests,the cardinality only version of the PHD filter has been used,hence only the number of targets is estimated.In the filter there is hence no explicit state dependence.As explained above,this makes the analysis independent of specific solutions for problems of first orbit determination and propagation.However,it already illustrates the general filter performance under model mismatch.In order to adequately supply the Cardinality only filter,all state dependent values are assumed to be identical for each target.

    The scenario is chosen as the optical observation of the ASTRA satellite cluster.The ASTRA cluster is located in the geostationary ring.It is assumed that the telescope is in staring mode and the objects are in the field of view at all times,displayed in Fig.3b.The scenario comprises of 100 seconds and measurements every 10 seconds.The phase angle can be assumed to be nearly constant during such a small time period,the solar panels are roughly aligned with the sun direction,and a three degree offset is assumed.The direction to the observer is 80 degrees.The angles are assumed fixed in the filter.Motion of the sun and the observer are negligible in the 100 second scenario.ASTRA satellites are at a distance of 36,000 km.The solar panel area(which dominates the reflection)can be assumed to be around 45 square meters(based of the dimensions of retired ASTRA 1B).Solar panels have an absorption coefficient of around 0.3,a specular component of around 0.5,and a diffuse component of 0.2.Under these observation conditions,the ASTRA satellites have a magnitude of 12.ASTRA satellites are large space objects and they appear very bright on the night sky under small observation angles.The current large observation angle is chosen to illustrate the effect of different sensors even on relatively large,stabilized(constant angle to the sun)objects.In the simulations,a gain of 0.1 is assumed,a readout noise of 100 ADU,a sky background of 5,000,000 ADU,and an exposure time of 1 second is modeled;for the background estimation,it is assumed the whole detector of 2000 by 2000 pixels have been used.Each scenario is run 1000 times.In order to measure and illustrate the performance,the average of all runs is plotted,and the mean sample error can be evaluated:

    where m(t)is the true number of objects,five in the chosen scenarios,?m(t)is the estimated number of targets,and the number of runs in the simulation,n,is 1000.It has to be noted that within the filter,non-integer values for the number of objects are allowed.Accordingly,the sample variance is computed:

    Scenario 1 is the baseline scenario.No clutter is assumed,the number of true objects stays constant combined with a high probability of detection of 1.0,which corresponds to a large aperture telescope observation of one meter.In the second scenario,realistic large aperture observations are assumed with a high quality camera,which means that clutter is also present.In the third scenario,a small aperture is assumed and camera performance parameters are slightly decreased, where either of those factors decreases the probability of detection.Clutter is also present.In scenarios 4 and 5,variable probability of detection is modeled,which corresponds to the ASTRA satellites after attitude control has been disabled.Natural torques lead to attitude motion,which is reflected in a varying signal to noise ratio and,ultimately,in a non-static probability of detection.A mid-sized aperture is assumed.In scenario 4,the filter models the change in the probability of detection exactly according to the truth model,in scenario 5,the filter assumes a constant probability of detection,whereas the true probability of detections is time-varying,same as in scenario 4.

    5.1 Scenario 1

    Scenario one is a baseline scenario;a one meter aperture telescope is assumed with 10 cm obstruction and a very high quantum efficiency of 0.97.The probability of detection is 1.0,despite the large observation angle,and the signal to noise ratio at the detector is 23.2.A detection threshold of t=2.5 is assumed.The five objects in the cluster are in the field of view during the whole scenario,and no objects are created nor lost.Object birth is assumed to be practically zero(10?15),and probability of survival is 1.0.It is assumed that no clutter would be present and that the filter is initialized with the correct number of objects(5),both unrealistic assumptions.Apart from the fixed number of objects,all assumptions of the filter are met.Not surprisingly,the filter performance is also perfect;at every time of the scenario,the correct cardinality is estimated,also illustrated in 5.

    5.2 Scenario 2

    Figure 6:Scenario 2,one arbitrary solution in the Monte Carlo trials:True number of objects(red),the cardinality of the measurements(true measurements plus clutter every ten seconds(blue))and the filter estimate(dashed green).

    In the second scenario,clutter was added.The clutter is modeled to be Poisson distributed with a parameter of 4.Probability of detection is still 1.0.The remaining parameters are modeled like in scenario 1.This would correspond to observations with a large one meter aperture telescope and 10 centimeter obstruction,with a realistic clutter process caused by cosmic ray events.The actual number of cosmic ray events varies by location and sensor specifics,as some sensor materials themselves are very radioactive on a low level.A Monte Carlo simulation of 1000 trials has been conducted.Fig.6 shows one trial scenario outcome from the Monte Carlo simulations,where the estimated number of objects are rounded to the closest integer in the display(inside the filter,floating point numbers of the estimated number of objects for the computation are kept).It shows that in the mismatch of the object population distribution,despite the high and correctly modeled probability of detection of 1.0,the filter does not match the true object population.Fig.7 shows the average over all trials without rounding.The error,standard deviation,and the deviation from the true number of observations.As the probability of detection is one,the number of observations is overshooting.Not surprisingly,the average number of observations above the correct number of objects is very close to the Poisson parameter of 4.The error in the number of objects that are estimated is not stable.The one sigma standard deviation does not envelope the true averaged error that is found in the simulations.The error in the estimate reaches around six percent maximally in the scenario.

    Figure 7:Scenario 2:(a)sample mean error(red),one standard deviation(green);(b)difference in the average number of measurements and the true object number.

    5.3 Scenario 3

    In scenario three,it is assumed that a smaller aperture telescope of 10 centimeter is used.The obstruction is 0.5 centimeters.The camera is assumed to have a quantum efficiency of 0.8.The signal to noise ratio in this case would be 2.2.A detection threshold of t=2.5 is assumed.This leads to a probability of detection of around 0.3.Fig.8 shows one sample trial of the Monte Carlo simulation.Because of the low probability of detection,the number of objects is consistently underestimated.

    Figure 8:Scenario 3,one arbitrary solution in the Monte Carlo trials:True number of objects(red),the cardinality of the measurements(true measurements plus clutter every ten seconds(blue))and the filter estimate(dashed green).

    Figure 9:Scenario 3:(a)sample mean error(red),one standard deviation(green);(b)difference in the average number of measurements and the true object number.

    Fig.9a shows the error and standard deviation,and Fig.9b shows the difference in the number of measurements and the true number of objects that are in the field of view.Because of the low probability of detection,on average,six to seven observations are made(true objects and clutter)with a clutter Poisson distribution with parameter four.The error statistic shows that the estimate is consistently off and seems to worsen the longer the scenario lasts.The error is significantly larger as in the previous case and builds up to over 300 percent.The model mismatch consists of the static object numbers and the assumption of a Poisson process in the predicted number of objects.The filter is stressed by the very low probability of detection in the presence of significant clutter.However,the scenario is not unrealistic.Low probabilities of detection can,as outlined above,even occur for large objects,under large observation or phase angles,even when the objects are stabilized and the angle to the sun is optimal.Small apertures lead to lower probability of detections in optical sensors than larger apertures for the same observation conditions.The performance of the sensor also has an influence,albeit to a lesser degree.The ASTRA satellites are among the largest objects in the geostationary ring,and the area of the solar panels is tremendous.Any smaller object with comparable reflection values,or less reflective objects of the same size,will have a lower signal to noise ratio,which in turn then can lead to a low probability of detection.Exceptions are glint conditions.

    5.4 Scenario 4

    In the fourth scenario,a mid-sized aperture is assumed of 50 centimeters and a 3 centimeter obstruction.Instead of assuming stabilized satellites that have a constant angle towards the sun and a slowly changing angle towards the observer,it is assumed that the satellites slowly rotate around the axis perpendicular to the solar panels.In the absence of stable active attitude control,natural torques induce rotation.The situation is also comparable to spin stabilized satellites with a very heterogeneous outer surfaces.The satellite solar panels are assumed to move from an angle of 35 degrees to the observer and 25 degrees to the sun,continuously to 85 degrees to the observer and 75 degrees to the sun,respectively.With a quantum efficiency of the camera of 0.9,this leads to a change in the signal to noise ratios of 97 to 3.The probability of detection accordingly changes from initially 1.0 for seven of the ten measurements and then reaches values of 0.99,0.87,and 0.23.A detection threshold of 2.5 is assumed again.In the filter,the truth of the probability of detection is exactly matched.Fig.10 shows a sample trial of the scenario.It can be seen that the filter changes from underestimating to overestimating the number of objects when the probability of detection changes from values of one in the first measurements to lower values,although the filter is updated with the correct values for the probability of detection.Fig.11 shows the that consistently the error is increasing as the probability of detections are lowered values towards the end of the scenario.The number of measurements on average is stable around four,the Poisson parameter of the clutter process,and then drops as not all objects are detected any more.Overall,the performance of the filter is better than with consistently low probability of detection values.The error reaches only around 7 percent.

    Figure 10:Scenario 4,one arbitrary solution in the Monte Carlo trials:True number of objects(red),the cardinality of the measurements(true measurements plus clutter every ten seconds(blue))and the filter estimate(dashed green).

    Figure 11:Scenario 4:(a)sample mean error(red),one standard deviation(green);(b)difference in the average number of measurements and the true object number.

    5.5 Scenario 5

    Figure 12:Scenario 5,one arbitrary solution in the Monte Carlo trials:True number of objects(red),the cardinality of the measurements(true measurements plus clutter every ten seconds(blue))and the filter estimate(dashed green).

    Figure 13:Scenario 5:(a)sample mean error(red),one standard deviation(green);(b)difference in the average number of measurements and the true object number.

    Scenario five is a replica of scenario four,except that the filter is provided with a constant value of the probability of detection of the correct but constant average value for the scenario of 0.91.For the generation of the measurements,the correct probability of detection sequence of 1.0 for seven of the ten measurements and the values of 0.99,0.87,and 0.23 for the remaining three measurement times,respectively,is used. Fig.12 shows one trial of the Monte Carlo simulations.It can be seen that in this case the number of objects is constantly either over or under-estimated and that the estimation error increases in the last observation,where the true probability of detection and the averaged one used in the filter differ the most.This trend that is visible in the single trial is supported by the complete Monte Carlo simulation.The mean sample error is larger throughout the scenario compared(be aware of the different scales)to the same scenario with the correctly modeled probability of detection(scenario 4).Although the performance is still in the same order of magnitude at the beginning of the scenario when the large probability of detections are dominated,towards the end of the scenario,when the small probability of detections are reached,the error dramatically increases and reaches values of 40 percent.It appears that the filter is performing comparably with the averaged to the true probability of detection as long as truth and average are both high and relatively close to one another.The filter performance worsens when truth and average probability of detection differ significantly.

    6 Conclusions

    Finite set statistics based filters have advantages because they offer a fully probabilistic description of a tracking scenario.First order approximations,such as in the cardinality only PHD,have become increasingly popular.In this paper,the filter assumptions that are made in the PHD filter derivation have been compared to the reality of space situational awareness(SSA)tracking.It has been found that,especially in dedicated tracking,a potential mismatch between the predicted object distribution might exist.Probability of survival and object birth are put to their extreme in SSA,especially in the absence of collisions.Collisions pose specific challenges,such as,among others,a model mismatch in target independence.Measurement independence assumptions require careful measurement treatment and potential redefinition especially in radar measurements.A special focus has been laid on the correct definition of the probability of detection.Scenarios have been run,mimicking a simple SSA tracking of five objects of the ASTRA cluster in the geostationary ring.For the simulations,the cardinality only part of the PHD filter has been utilized.Although state dependence has been lost,it provides the advantage to investigate crucial parts of the full PHD filter performance without committing to and independent of a specific method for object propagation and initial orbit determination.

    The simulations suggest that,in the presence of realistic clutter,the mean sample error is of the order of ten percent,even when the probability of detection is one,and object birth and death are perfectly known and the predicted object population is the only potential model mismatch.Performance significantly worsens with low probability of detection,leading to mean sample errors of several hundred percent,even when perfectly matched in the filter.Using a time varying probability of de-tection leads to consistent results with the constant probability of detection.Results are significantly worsened when an average probability of detection is used,especially when the average and the true probability of detection differ significantly.In all cases,the standard variation was small relative to the mean sample error.

    Acknowledgement:I would like to acknowledge the fruitful collaboration with Prof.Dr.Daniel Clark and his research group,in particular Dr.Emmanuel Delande,from Heriot Watt University,Edinburgh,UK.I am very thankful for all the useful discussions.The work was supported by the VIBOT summer scholar program.

    Allen,C.(1973):Astrophysical Quantities.The Athlone Press,University of London,London.ISBN:0 485 11150 0.

    Bar-Shalom,Y.(1978): Tracking methods in a multitarget environment.Automatic Control,IEEE Transactions on,vol.23,no.4,pp.618–626.

    Blackman,S.S.;Popoli,R.(1999):Design and analysis of modern tracking systems.Artech House.

    Blair,W.D.;Bar-Shalom,Y.(Eds):Multitarget-Multisensor Tracking:Applications and Advances(Volume III).Artech House.

    Cheng,Y.;DeMars,K.J.;Frueh,C.;Jah,M.K.(2013):Gaussian mixture phdfilter for space object tracking.No.AAS 13-242 in Advances in the Astronautical Sciences.

    Cheng,Y.;Frueh,C.;DeMars,K.J.(2013):Comparisons of phd filter and cphdfilter for space object tracking.No.AAS 13-770 in Advances in the Astronautical Sciences.

    Clark,D.E.;Panta,K.;Vo,B.-N.(2006):The GM-PHD Filter Multiple Target Tracker.InInformation Fusion,Proceedings of the 9th International Conference on,pp.1–8.

    Daniels,G.(1977): A Night Sky Model for Satellite Search Systems.Optical Engineering No 1,vol.16,pp.66–71.

    DeMars,K.J.;Hussein,I.I.;Frueh,C.;Jah,M.K.;Erwin,R.S.(2015):Multiple-object space surveillance tracking using finite-set statistics.Journal of Guidance,Control,and Dynamics.doi:10.2514/1.G000987.

    Flohrer,T.;Schildknecht,T.;R.St?veken,E.(2005):Performance Estimation for GEO Space Surveillance.Advances in Space Research,vol.35,no.7,pp.1226–1235.

    Frueh,C.(2011):Identification of Space Debris.Shaker Verlag,Aachen.ISBN:978-3-8440-0516-5.

    Frueh,C.;Jah,M.(2013):Detection Probability of Earth Orbiting Using Optical Sensors.InProceedings of the AAS/AIAA Astrodynamics Specialist Conference,Hilton Head,SC.13-701.

    Frueh,C.;Schildknecht,T.(2012): Object Image Linking of Objects in Near Earth Orbits in the Presence of Cosmics.Advances in Space Research,vol.49,pp.594–602.

    Gehly,S.;Jones,B.;Axelrad,P.(2014):An AEGIS-CPHD Filter to Maintain Custody of GEO Space Objects with Limited Tracking Data.InAdvanced Maui Optical and Space Surveillance Technologies Conference,pg.25.

    Houssineau,J.;Delande,E.;Clark,D.E.(2013):Notes of the Summer School on Finite Set Statistics.Technical report,2013.arXiv:1308.2586.

    Jia,B.;Pham,K.;Blasch,E.;Shen,D.;Wang,Z.;Chen,G.(2014): Cooperative space object tracking using consensus-based filters.InInformation Fusion(FUSION),2014 17th International Conference on,pp.1–8.

    Jones,B.;Vo,B.-N.(2015):A Labelled Multi-Bernoulli Filter for Space Object Tracking.InProceedings of the 23rd AAS/AIAA Spaceflight Mechanics Meeting,Williamsburg,VA,January 11-15.

    Kelecy,T.;Jah,M.;DeMars,K.(2012): Application of a multiple hypothesisfilter to near GEO high area-to-mass ratio space objects state estimation.Acta Astronautica,vol.81,pp.435–444.

    Kessler,D.;Cour-Palais,B.(1978): The creation of a debris belt.Journal of Geophysical Research,vol.83,no.A6,pp.2637–2646.

    Klinkrad,H.;Sdunnus,H.(1997):Concepts and applications of the{MASTER}space debris environment model.Advances in Space Research,vol.19,no.2,pp.277–280.

    Klinkrad,H.;Wegener,P.;Wiedemann,C.;Bendisch,J.;Krag,H.(2006):Space Debris:Models and Risk Analysis,chapter Modeling of the Current Space Debris Environment.Springer Berlin Heidelberg,2006.

    Krag,H.;Klinkrad,H.;Jehn,R.;L.Leushacke;Markkanen,J.(2007):Detection of small-size space debris with the fgan and eiscat radars.In7th Workshop,Naval Postgraduate School Monterey,California.

    Kreucher,C.M.;Kastella,K.D.;Hero III,A.O.(2004): Tracking Multiple Targets Using a Particle Filter Representation of the Joint Multitarget Probability Density.InSignal and Data Processing of Small Targets,Proceedings of SPIE,volume 5204,pp.258–269.

    Merline,W.;Howell,S.(1995): A realistic model for point-sources imaged on array detectors:The model and initial results.Experimental Astronomy,vol.6,pp.163–210.

    Mahler,R.P.S.(2007):Statistical Multisource-Multitarget Information Fusion.Artech House.

    Markkanen J.,Postila M.,v.E.A.(2006):Small-size Space Debris Data Collection with EISCAT Radar Facilities.Technical report,2006.Final Report of ESOC Contract No.18575-/04/D/HK(SC).

    Musci,R.(2006):Identification and Recovery of Objects in GEO and GTO to Maintain a Catalogue of Orbits.Astronomical Institute,University of Bern.PhD thesis.

    Musci,R.;Schildknecht,T.;Ploner,M.(2004): Orbit Improvement for GEO Objects Using Follow-up Obervations.Advances in Space Research,vol.34,no.5,pp.912–916.

    Musci,R.;Schildknecht,T.;Ploner,M.;Beutler,G.(2005):Orbit Improvement for GTO Objects Using Follow-up Obervations.Advances in Space Research,vol.35,no.7,pp.1236–1242.

    Oswald,M.;Krag,H.;Wegener,P.;Bischof,B.(2004): Concept for an Orbital Telescope Observing the Debris Environment in GEO.Advances in Space Research,vol.34,no.5,pp.1155–1159.

    Potter,A.(1995): Ground-Based Optical Observations of Orbital Debris:A Review.Advances in Space Research,vol.16,no.11,pp.(11)35–(11)45.

    Reid,D.(1979):An Algorithm for Tracking Multiple Targets.Automatic Control,IEEE Transactions on,vol.24,no.6,pp.843–854.

    Sanson,F.;Frueh,C.(2015): Noise Quantification in Optical Observations of Resident Space Objects for Probability of Detection and Likelihood.InAAS/AIAA Astrodynamic Specialist Conference,Vail,CO.15-634.

    Sanson,F.;Frueh,C.(2016): Probability of Detection and Likelihood:Application to Optical Unresolved Space Object Observation.Celestial Mechanics and Dynamical Astronomy.

    Singh,N.;Poore,A.;Sheaff,C.;Aristoff,J.;Jah,M.(2013):Multiple Hypothesis Tracking(MHT)for Space Surveillance:Results and Simulation Studies.InAdvanced Maui Optical and Space Surveillance Technologies Conference,pg.16.

    US Air Force(2016):Celestrack Database based on two line elements,supported by CSSI.http://celestrak.com,2016a.

    US Air Force(2016): United States Strategic Command Website.http://www.stratcom.mil/factsheets/jspoc/,2016b.

    Vo,B.-N.;Ma,W.-K.(2006):The gaussian mixture probability hypothesis density filter.Signal Processing,IEEE Transactions on,vol.54,no.11,pp.4091–4104.

    Wertz,J.(1978):Spacecraft Attitude Determination and Control.Volume 73.D.Reidel Publishing Company,Dordrecht:Holland.ISBN:90-277-0959-9.

    1Purdue University,West-Lafayette,IN,USA.

    欧美激情久久久久久爽电影| 内射极品少妇av片p| 亚洲成a人片在线一区二区| 久久99精品国语久久久| 欧美成人一区二区免费高清观看| 欧美日韩在线观看h| 欧美日本亚洲视频在线播放| 日产精品乱码卡一卡2卡三| 久久久久久久久久久免费av| 欧美潮喷喷水| 99久国产av精品| 国产国拍精品亚洲av在线观看| av天堂中文字幕网| 久久鲁丝午夜福利片| 两性午夜刺激爽爽歪歪视频在线观看| 午夜精品一区二区三区免费看| 亚洲久久久久久中文字幕| av女优亚洲男人天堂| 一区二区三区免费毛片| 精华霜和精华液先用哪个| 中文字幕av在线有码专区| 99久久精品国产国产毛片| av在线播放精品| 天美传媒精品一区二区| 亚洲国产欧洲综合997久久,| 日本-黄色视频高清免费观看| 少妇的逼水好多| 欧美高清性xxxxhd video| 精品日产1卡2卡| 国产片特级美女逼逼视频| 欧美成人一区二区免费高清观看| 一个人免费在线观看电影| 亚洲最大成人av| 精品久久久久久久久av| 99久久中文字幕三级久久日本| 亚洲激情五月婷婷啪啪| 亚洲第一区二区三区不卡| 波野结衣二区三区在线| 黄色一级大片看看| 看免费成人av毛片| 亚洲欧美日韩东京热| 国产精品国产高清国产av| 久久鲁丝午夜福利片| 99视频精品全部免费 在线| 亚洲成人精品中文字幕电影| 亚洲av.av天堂| 少妇裸体淫交视频免费看高清| 99国产精品一区二区蜜桃av| 国产精品一二三区在线看| 色尼玛亚洲综合影院| 赤兔流量卡办理| 尤物成人国产欧美一区二区三区| av女优亚洲男人天堂| 老女人水多毛片| 人妻制服诱惑在线中文字幕| 精品人妻视频免费看| 99久国产av精品国产电影| 国产精品一区二区性色av| 插逼视频在线观看| 天堂网av新在线| 亚洲成a人片在线一区二区| 精品久久久久久久末码| 亚洲精品国产成人久久av| 欧美bdsm另类| 国产一区二区在线观看日韩| 亚洲欧洲国产日韩| 丰满的人妻完整版| 男女做爰动态图高潮gif福利片| 亚洲成人精品中文字幕电影| 亚洲国产高清在线一区二区三| 女同久久另类99精品国产91| 国产极品精品免费视频能看的| 国产极品天堂在线| 边亲边吃奶的免费视频| 午夜a级毛片| 婷婷六月久久综合丁香| 日日干狠狠操夜夜爽| 美女 人体艺术 gogo| 又黄又爽又刺激的免费视频.| 白带黄色成豆腐渣| 国内精品美女久久久久久| 九九久久精品国产亚洲av麻豆| av在线天堂中文字幕| 1024手机看黄色片| 国产探花在线观看一区二区| .国产精品久久| 国产乱人视频| 久久人人爽人人爽人人片va| 99久国产av精品国产电影| 别揉我奶头 嗯啊视频| 变态另类成人亚洲欧美熟女| 久久婷婷人人爽人人干人人爱| 亚洲av一区综合| 欧美xxxx性猛交bbbb| 国产色爽女视频免费观看| 麻豆成人av视频| 黄片无遮挡物在线观看| 超碰av人人做人人爽久久| 天天躁日日操中文字幕| 午夜爱爱视频在线播放| 日韩 亚洲 欧美在线| 极品教师在线视频| kizo精华| av天堂在线播放| 又粗又硬又长又爽又黄的视频 | h日本视频在线播放| 久久99蜜桃精品久久| 色哟哟·www| 国产亚洲精品久久久com| 一本精品99久久精品77| 免费av观看视频| 最好的美女福利视频网| 久久国产乱子免费精品| 长腿黑丝高跟| 最近的中文字幕免费完整| 日韩视频在线欧美| 国产爱豆传媒在线观看| 日韩三级伦理在线观看| 禁无遮挡网站| 久久久久久久久久久丰满| 国产精品一区二区在线观看99 | 免费电影在线观看免费观看| 亚洲国产欧美在线一区| 久久中文看片网| 久久久精品94久久精品| 麻豆成人午夜福利视频| 日韩欧美在线乱码| 免费观看a级毛片全部| 舔av片在线| 国产精品嫩草影院av在线观看| 国产精品麻豆人妻色哟哟久久 | 成人高潮视频无遮挡免费网站| 日本与韩国留学比较| 色综合站精品国产| 国产黄色视频一区二区在线观看 | 观看免费一级毛片| 日本黄色视频三级网站网址| 国产精品免费一区二区三区在线| 在线免费十八禁| 日韩视频在线欧美| 2022亚洲国产成人精品| 精品午夜福利在线看| 高清毛片免费观看视频网站| 日韩大尺度精品在线看网址| ponron亚洲| 黄片wwwwww| 岛国在线免费视频观看| 久久草成人影院| 不卡视频在线观看欧美| 精品一区二区免费观看| 成人漫画全彩无遮挡| 女人十人毛片免费观看3o分钟| 亚洲精品影视一区二区三区av| 久久草成人影院| 两性午夜刺激爽爽歪歪视频在线观看| 精品人妻视频免费看| 久久99热这里只有精品18| 欧美一区二区亚洲| 日韩av在线大香蕉| 精品欧美国产一区二区三| 国产综合懂色| 成年版毛片免费区| 国产精品av视频在线免费观看| 如何舔出高潮| 91狼人影院| 日韩精品青青久久久久久| 又粗又硬又长又爽又黄的视频 | 日本在线视频免费播放| 在线观看免费视频日本深夜| 中文精品一卡2卡3卡4更新| 99riav亚洲国产免费| 日本三级黄在线观看| 内射极品少妇av片p| 少妇被粗大猛烈的视频| 校园人妻丝袜中文字幕| 直男gayav资源| 国产精品麻豆人妻色哟哟久久 | 日韩 亚洲 欧美在线| 中文字幕制服av| 你懂的网址亚洲精品在线观看 | 美女xxoo啪啪120秒动态图| 精品欧美国产一区二区三| kizo精华| 久久精品久久久久久噜噜老黄 | 一级毛片我不卡| 禁无遮挡网站| 亚洲欧美中文字幕日韩二区| 日日啪夜夜撸| 变态另类成人亚洲欧美熟女| 我的老师免费观看完整版| 欧美不卡视频在线免费观看| 国产精品伦人一区二区| 色视频www国产| 国产精品爽爽va在线观看网站| 久久久精品欧美日韩精品| 亚洲内射少妇av| 免费看美女性在线毛片视频| 嫩草影院精品99| 亚洲欧洲日产国产| 三级毛片av免费| 97超碰精品成人国产| 国产精品久久久久久亚洲av鲁大| 亚洲av成人av| 亚洲欧美日韩高清专用| 久久久久免费精品人妻一区二区| 黄片wwwwww| 国产高清有码在线观看视频| 老女人水多毛片| 你懂的网址亚洲精品在线观看 | 免费观看精品视频网站| 国产免费一级a男人的天堂| 久久精品影院6| 亚洲欧美中文字幕日韩二区| av专区在线播放| 日日啪夜夜撸| 亚洲国产精品成人综合色| 国内久久婷婷六月综合欲色啪| h日本视频在线播放| 国产一区亚洲一区在线观看| 给我免费播放毛片高清在线观看| 少妇高潮的动态图| 国产精品日韩av在线免费观看| 国产成人aa在线观看| 不卡视频在线观看欧美| 日韩欧美一区二区三区在线观看| 欧美色视频一区免费| 日韩一区二区视频免费看| 国产一级毛片七仙女欲春2| 精品久久久久久成人av| 亚洲av电影不卡..在线观看| 欧美精品一区二区大全| 一级毛片久久久久久久久女| 免费人成视频x8x8入口观看| 国产久久久一区二区三区| 成年版毛片免费区| 亚洲一级一片aⅴ在线观看| 亚洲成人中文字幕在线播放| 国产亚洲5aaaaa淫片| 亚洲自拍偷在线| 中文精品一卡2卡3卡4更新| 久久韩国三级中文字幕| 国产精品伦人一区二区| 精品国内亚洲2022精品成人| 国产精品1区2区在线观看.| 日本五十路高清| 亚洲自拍偷在线| 国产一区二区激情短视频| 久久午夜福利片| 国产69精品久久久久777片| 亚洲国产精品国产精品| 伦理电影大哥的女人| 国产精品久久久久久亚洲av鲁大| 国产精品免费一区二区三区在线| 亚洲va在线va天堂va国产| 色视频www国产| 国产成人精品一,二区 | 最近视频中文字幕2019在线8| 狠狠狠狠99中文字幕| 亚洲欧洲日产国产| 久久久久性生活片| 国产精品国产高清国产av| 婷婷六月久久综合丁香| 99久久精品国产国产毛片| 中文字幕av在线有码专区| 少妇丰满av| 老司机影院成人| 91狼人影院| 亚洲综合色惰| 国产成人aa在线观看| 国产av一区在线观看免费| 干丝袜人妻中文字幕| 亚洲av.av天堂| av视频在线观看入口| 欧美日韩在线观看h| 亚洲精品粉嫩美女一区| 久久久久久久午夜电影| 亚洲国产欧美人成| 久久久久九九精品影院| 如何舔出高潮| 免费看日本二区| 蜜桃久久精品国产亚洲av| 国产av麻豆久久久久久久| 亚洲av二区三区四区| 亚洲国产高清在线一区二区三| 久久精品国产亚洲av香蕉五月| 99在线人妻在线中文字幕| 最近视频中文字幕2019在线8| 精品欧美国产一区二区三| 最近的中文字幕免费完整| 午夜福利成人在线免费观看| 亚洲成av人片在线播放无| 免费观看a级毛片全部| 99九九线精品视频在线观看视频| 国产日韩欧美在线精品| 国产精品久久久久久亚洲av鲁大| 狂野欧美激情性xxxx在线观看| 免费看av在线观看网站| 国产极品精品免费视频能看的| 哪个播放器可以免费观看大片| 久久九九热精品免费| 日日摸夜夜添夜夜添av毛片| 亚洲欧美成人综合另类久久久 | 国产精品三级大全| 国产精品久久久久久精品电影小说 | www.色视频.com| 一级毛片久久久久久久久女| 麻豆av噜噜一区二区三区| 欧美精品国产亚洲| 国产伦在线观看视频一区| 婷婷精品国产亚洲av| 亚洲欧美精品自产自拍| 亚洲中文字幕日韩| 久久99热这里只有精品18| 色播亚洲综合网| 男人狂女人下面高潮的视频| 蜜桃久久精品国产亚洲av| 亚洲国产日韩欧美精品在线观看| 两个人的视频大全免费| 久久综合国产亚洲精品| 秋霞在线观看毛片| 成人欧美大片| 人妻夜夜爽99麻豆av| 色哟哟·www| 免费在线观看成人毛片| 国产成人freesex在线| 日韩在线高清观看一区二区三区| 亚洲精品自拍成人| 亚洲av免费在线观看| 精品久久久久久久末码| 免费在线观看成人毛片| 日日摸夜夜添夜夜添av毛片| av在线观看视频网站免费| 精华霜和精华液先用哪个| 日本爱情动作片www.在线观看| 最近2019中文字幕mv第一页| 国产日本99.免费观看| 亚洲欧美精品综合久久99| 国产亚洲精品av在线| 伊人久久精品亚洲午夜| 黄片无遮挡物在线观看| 一级毛片久久久久久久久女| 中文字幕免费在线视频6| 最近最新中文字幕大全电影3| 岛国毛片在线播放| 欧美+亚洲+日韩+国产| 亚洲四区av| 午夜久久久久精精品| 久久鲁丝午夜福利片| 国产极品天堂在线| 69人妻影院| 日韩欧美国产在线观看| 22中文网久久字幕| 欧美一区二区亚洲| 干丝袜人妻中文字幕| 国语自产精品视频在线第100页| 97超视频在线观看视频| 高清毛片免费观看视频网站| 国产大屁股一区二区在线视频| 成年女人看的毛片在线观看| 97超视频在线观看视频| 国产人妻一区二区三区在| 小蜜桃在线观看免费完整版高清| 国产蜜桃级精品一区二区三区| 乱系列少妇在线播放| 91精品一卡2卡3卡4卡| 99热网站在线观看| 免费看日本二区| 在线观看av片永久免费下载| 亚洲国产色片| 色综合色国产| 偷拍熟女少妇极品色| 午夜久久久久精精品| 91aial.com中文字幕在线观看| 欧美色视频一区免费| 少妇熟女aⅴ在线视频| 成人特级av手机在线观看| 久久99热6这里只有精品| 日产精品乱码卡一卡2卡三| 午夜老司机福利剧场| 久久久久免费精品人妻一区二区| av又黄又爽大尺度在线免费看 | 97人妻精品一区二区三区麻豆| 狠狠狠狠99中文字幕| 如何舔出高潮| 久久久久九九精品影院| 亚洲精品成人久久久久久| 男女做爰动态图高潮gif福利片| 国产精品,欧美在线| 97热精品久久久久久| 亚洲成a人片在线一区二区| 午夜福利成人在线免费观看| 国产又黄又爽又无遮挡在线| 久久久午夜欧美精品| 欧美日韩乱码在线| 九九爱精品视频在线观看| 美女脱内裤让男人舔精品视频 | 99久久精品热视频| 国产成人影院久久av| 日本在线视频免费播放| 全区人妻精品视频| 国产高清视频在线观看网站| 老司机影院成人| 少妇高潮的动态图| av又黄又爽大尺度在线免费看 | 国内揄拍国产精品人妻在线| 91在线精品国自产拍蜜月| 天堂网av新在线| www.色视频.com| 性欧美人与动物交配| 国产精品1区2区在线观看.| 欧美xxxx黑人xx丫x性爽| 伦理电影大哥的女人| 偷拍熟女少妇极品色| 久久久久久久久久成人| 乱人视频在线观看| 97超碰精品成人国产| 91麻豆精品激情在线观看国产| 婷婷色av中文字幕| 午夜福利在线观看免费完整高清在 | 亚洲婷婷狠狠爱综合网| 淫秽高清视频在线观看| 性欧美人与动物交配| 久久久久久久久久成人| 97超视频在线观看视频| 成人亚洲精品av一区二区| 色哟哟哟哟哟哟| av卡一久久| 久久久久免费精品人妻一区二区| 丰满人妻一区二区三区视频av| 久久久久性生活片| 欧美精品一区二区大全| av卡一久久| 可以在线观看毛片的网站| 人妻夜夜爽99麻豆av| 国产精品三级大全| 在线播放国产精品三级| 亚洲成人久久性| 综合色av麻豆| 欧洲精品卡2卡3卡4卡5卡区| 99久久久亚洲精品蜜臀av| 综合色丁香网| 夜夜夜夜夜久久久久| 久久久久免费精品人妻一区二区| 热99在线观看视频| 在线免费十八禁| 中文字幕av在线有码专区| 国内久久婷婷六月综合欲色啪| 日韩,欧美,国产一区二区三区 | 91午夜精品亚洲一区二区三区| avwww免费| 午夜精品在线福利| 老熟妇乱子伦视频在线观看| 日韩欧美精品免费久久| 伊人久久精品亚洲午夜| 亚洲av一区综合| 久久这里有精品视频免费| 九九热线精品视视频播放| 亚洲在线观看片| 丰满人妻一区二区三区视频av| 色尼玛亚洲综合影院| 黄片wwwwww| 午夜精品国产一区二区电影 | 亚洲中文字幕日韩| 乱人视频在线观看| 91精品国产九色| 国产极品精品免费视频能看的| 亚洲内射少妇av| 51国产日韩欧美| 黄色配什么色好看| 日韩欧美三级三区| 男人和女人高潮做爰伦理| a级一级毛片免费在线观看| 日韩精品青青久久久久久| 高清日韩中文字幕在线| 午夜精品在线福利| 免费人成在线观看视频色| 亚洲真实伦在线观看| 69av精品久久久久久| 成人综合一区亚洲| 自拍偷自拍亚洲精品老妇| 一区二区三区免费毛片| 亚洲欧美成人精品一区二区| 久久这里只有精品中国| 成人av在线播放网站| 伊人久久精品亚洲午夜| 亚洲自偷自拍三级| 国产v大片淫在线免费观看| 成人特级黄色片久久久久久久| 一级毛片电影观看 | 一边摸一边抽搐一进一小说| 看十八女毛片水多多多| 亚洲五月天丁香| 免费看光身美女| 久久综合国产亚洲精品| 免费av不卡在线播放| 精品久久久久久久末码| 亚洲第一区二区三区不卡| 中国美女看黄片| 三级男女做爰猛烈吃奶摸视频| 国产一区二区激情短视频| 99热这里只有是精品在线观看| 小说图片视频综合网站| 国产黄a三级三级三级人| 在线免费观看不下载黄p国产| 午夜免费男女啪啪视频观看| 美女脱内裤让男人舔精品视频 | 亚洲无线在线观看| 欧美变态另类bdsm刘玥| 欧美性猛交╳xxx乱大交人| 国产精品不卡视频一区二区| 国产黄色小视频在线观看| www.av在线官网国产| 卡戴珊不雅视频在线播放| 在线天堂最新版资源| 亚洲欧美清纯卡通| 免费观看在线日韩| 黄色日韩在线| 老司机影院成人| 午夜福利在线在线| 晚上一个人看的免费电影| 97超视频在线观看视频| 久久久a久久爽久久v久久| 99久久中文字幕三级久久日本| 波多野结衣巨乳人妻| 噜噜噜噜噜久久久久久91| 日日摸夜夜添夜夜添av毛片| 欧美zozozo另类| 久久久久久久久久久免费av| 亚洲最大成人手机在线| 久久久久免费精品人妻一区二区| 日韩强制内射视频| 日本熟妇午夜| 午夜免费激情av| 精品一区二区三区人妻视频| 人妻少妇偷人精品九色| 国内揄拍国产精品人妻在线| 高清在线视频一区二区三区 | av在线蜜桃| 久久国内精品自在自线图片| 人妻夜夜爽99麻豆av| 成年免费大片在线观看| 日韩欧美精品v在线| 精品日产1卡2卡| 18禁黄网站禁片免费观看直播| 日本撒尿小便嘘嘘汇集6| 国产真实伦视频高清在线观看| 婷婷色av中文字幕| 美女黄网站色视频| 啦啦啦观看免费观看视频高清| 久久久精品94久久精品| 黄色配什么色好看| 性色avwww在线观看| 寂寞人妻少妇视频99o| 99久久精品国产国产毛片| 日韩精品青青久久久久久| 99热这里只有是精品在线观看| 亚洲精品日韩av片在线观看| 精品国内亚洲2022精品成人| 一级二级三级毛片免费看| 亚洲精华国产精华液的使用体验 | 黄色日韩在线| 午夜老司机福利剧场| 色哟哟哟哟哟哟| 如何舔出高潮| 99久久中文字幕三级久久日本| 国内久久婷婷六月综合欲色啪| 99热网站在线观看| av专区在线播放| 亚州av有码| 精品国内亚洲2022精品成人| 精华霜和精华液先用哪个| 少妇人妻一区二区三区视频| 亚洲国产色片| 成人毛片60女人毛片免费| 免费人成视频x8x8入口观看| a级毛片免费高清观看在线播放| 波多野结衣高清作品| 51国产日韩欧美| 不卡一级毛片| 中国国产av一级| 国产精品人妻久久久久久| 天美传媒精品一区二区| 色5月婷婷丁香| 丝袜美腿在线中文| 免费av不卡在线播放| 国产精品嫩草影院av在线观看| 蜜臀久久99精品久久宅男| 亚洲人与动物交配视频| 1000部很黄的大片| 久久国产乱子免费精品| 最近视频中文字幕2019在线8| 色尼玛亚洲综合影院| 亚洲精品日韩av片在线观看| 一本久久中文字幕| 看免费成人av毛片| a级毛片a级免费在线| 国产在视频线在精品| 此物有八面人人有两片| 男人的好看免费观看在线视频| 久久久欧美国产精品| 中文精品一卡2卡3卡4更新| 国产一级毛片在线| 最好的美女福利视频网| 午夜激情福利司机影院| 欧美+亚洲+日韩+国产| 亚洲人与动物交配视频| 99久久久亚洲精品蜜臀av| 春色校园在线视频观看| 在线天堂最新版资源| 亚洲五月天丁香| .国产精品久久| 欧美成人a在线观看| 久久鲁丝午夜福利片| 99九九线精品视频在线观看视频| 99精品在免费线老司机午夜|