Cone Parameter Estimation from Simulated Range/Doppler Data William DeMeo and Jon Guerin June 9, 2005

1

The Forward Model

The “forward model” describes the process by which a given spinning cone produces a range/doppler image. The model accounts for various cone parameters, most notably aspect (cone orientation), irradiance, spin rate, and cone shape (diameter and length). To date, our model has not accounted for the effects of speckle or wobbling. Given a fixed aspect (angle and range), and no speckle, the range/doppler images are thought to be unique, but we did not prove this. The forward model works by sampling points on the cone and calculating the contribution each point makes to the range/doppler image. The sampling is done by taking successive points on the cone surface separated by a fixed arc length distance, resulting in a roughly “evenly sampled” cone. The goal is to use the forward model to estimate three of the parameters (namely, spin rate, length and diameter). We parameter values describing a cone that yields a range/doppler image that is closest in norm (i.e. Euclidean distance) to the data (the observed range/doppler image).

1.1

Problem Statement and Parameters of Interest

Given is a range/cross-range image that we assume was generated by a cone having a fixed, known orientation relative to the laser radar system. (We shall call the latter “the telescope.”) Without loss of generality, we fix a coordinate system with origin at the center of the base of the cone, letting the base of the cone occupy the “yz-plane” and taking the major axis of the cone to be the x-axis.1 (See Fig. 1) To specify the orientation of the cone relative to the telescope we use a vector A, which we shall call the “aspect vector,” and which extends from the origin at the base of the cone to the telescope position. Cone Orientation: some implementation details. To represent that the cone is oriented relative to the telescope in such that we are looking at the side of the cone – i.e. the cone’s major axis is orthogonal to the vector that extends from the telescope to the base of the cone – we use an “aspect vector” A = [0; -dist; 0];, where the entry dist is some distance from the telescope to the base of the cone. When the cone is “head-on” – i.e. the cone’s major axis is parallel to the vector that extends from the telescope to the cone – we represent its orientation with an aspect vector A = [dist, 0, 0];. A 45 degree view of the cone is specified with aspect vector A = [angle, -angle, 0];. We ran simulations for each of these three cases, each time taking the cone aspect vector as given, and varying the three other cone parameters (i.e. length, spin-rate, and base-diameter parameters). 1 In particular, this simplies the geometry of the problem by assuming that the cone is always fixed relative to the coordinate system. If we want to represent that the orientation of the cone relative to the telescope has changed, we simply allow the telescope position to vary and leave the cone coordinates.

1

1.2

Forward Model in Detail

First we construct a “forward model” which, given a cone with fixed orientation and parameter values (L0 , S0 , D0 ), specifies the range/cross-range image that such a cone would generate. We can then use this forward model to solve the inverse problem. That is, given a range/cross-range image, which cone (or cones) could have produce such an image, according to the forward model.

Z

Y

X

Figure 1: .

Velocity. First we compute the velocity of any point on the surface of the cone. This is the product of the scalar speed times the unit motion vector (i.e., a unit vector pointing in the direction of motion). First we compute the speed. Suppose the spin rate of the cone (i.e., the rate at which the cone spins around its major axis) is S0 revolutions per second. For a point p = (x, y, z) on the cone surface, the speed at p depends upon S0 and the radius of a cross-section of the cone at p. Since we have fixed the major axis of the cone to be the x-axis, the radius at p = (x, y, z) depends only on x. Indeed, it is easy to see in Fig. 1 that the radius of the cone at p is r0 r(p) = r0 − x, L0 where r0 = D0 /2. (Recall, D0 is the cone’s base diameter, so r0 = D0 /2 is the radius at the base.) Therefore, the point p travels the circumference 2πr(p) of the cone S0 times per second, with speed   r0 S(p) = 2πr(p)S0 = 2π r0 − x S0 L0 To compute the unit motion vector at p, assume the cone spins in a counter-clockwise direction relative to a viewer looking down the positive x-axis at the yz-plane. Then the unit motion vector has a 0 x-component and is tangent to the cone at p: T1 (p) = (0, − sin(tan−1 (z/y)), cos(tan−1 (z/y))). 2

Therefore, the velocity vector at p is:   V(p) = S(p)T1 (p) = 2πr0 (1 − xL−1 0 )S0 T1 (p). Doppler.

The doppler at any point p on the cone surface is given by the inner product D(p) =

V(p) · A , λ0 kAk

where λ0 is the wavelength of the outgoing laser pulse. Range.

The range at any point p on the cone surface is given by ˜ R(p) = kA − pk.

However, for the purpose of creating the range/cross-range image, we center the range values about the origin and use the relative range, R(p) = kA − pk − kAk. Irradiance. The irradiance at any point p on the cone surface is proportional to the inner product of the outward-pointing unit normal to the surface at p and the unit vector direction of the light source (which, in this case, is the laser). Let n(p) denote the unit normal to the surface at p. We derive a formula for n(p) and then compute the irradiance as (proportional to) the inner product of n(p) with the normalized aspect vector A/kAk. We previously derived an equation for the tangent vector T1 . Once we compute a second vector T2 , that is tangent to the surface and orthogonal to T1 , the cross product T1 × T2 will yield the normal to the surface. It should be clear from Fig. 1 that the following tangent vector is the one we need: T2 (p) = T2 (x, y, z) =

(L0 − x, −y, −z) (L0 − x, −y, −z) =p kT2 (p)k (L0 − x)2 + y 2 + z 2

Therefore, using “the right-hand rule” to ensure that we get the outward -pointing vector, the normal to the surface is given by e1 e2 e3 n(p) = (T1 × T2 )(p) = T1 (1) T1 (2) T1 (3) T2 (1) T2 (2) T2 (3) where, for lack of better notation, Ti (j) denotes the jth entry of the vector Ti (p). Substituting for T1 and T2 and performing the cross product results in the unit normal vector   z sin(tan−1 (z/y)) + y cos(tan−1 (z/y))  (L0 − x) cos(tan−1 (z/y)) n(p) = α  −1 (L0 − x) sin(tan (z/y))) p where α = 1/ (L0 − x)2 + y 2 + z 2 . Finally, the irradiance is proportional to the inner product of the unit normal and the unit aspect. That is, I(p) ∝ n(p) · A/kAk. (1) The constant of proportionality is not important since it depends on the reflectivity characteristics of the cone surface at the point p. As long as we assume all points are the same in this regard, the right hand side of (1) will give irradiances values for all points that are correct relative to one another.

3

Range/Cross-range Image. For a given point (r, c) in the range/cross-range domain, we compute the range/cross-range value as follows: X G(r, c) = I + (p), p∈S(r,c)

where I

+

is the positive irradiance, I + (p) = max{0, I(p)},

and, if C ⊂ R3 denotes the set of points on the cone surface, S(r, c) = {p ∈ C : |R(p) − r| < r , |D(p) − c| < c }. We use the positive irradiance I + because negative values of I correspond precisely to those points on the cone surface that are hidden from view (i.e., on the back-side of the cone relative to the detector). The set S(r, c) is the collection of points on the cone surface that contribute irradiance values to this particular range/cross-range, (r, c).

1.3

Parameter Estimation

Since we assume the aspect vector A is known (and since the object is assumed to be a cone) there remain only three unknown parameters which must be estimated. These are the cone’s length L (meters), spin rate S (Hertz), and base diameter D (meters). We carried out the estimation in two different ways. First we used a Matlab optimization function that does not require us to supply the gradients of the range/cross-range function with respect to each parameter. Second, to compare with and verify the results of the Matlab optimizer, we wrote an optimization function (in spin cone estimate.m – Sec. 3.2 describes the algorithm) that relies on gradient information. For this purpose, rather than derive closed form expressions, we approximated the required gradients using a basic finite difference scheme. It may be possible to derive closed form gradient equations. However, as the dependence of the range/doppler image values on the parameters is defined implicitly in terms of an integral,2 the derivation is non-trivial.

2

Detected Data

The detected data is from RME experiments on 4/8/98. The pulse-burst laser radar data can be used to generate range/doppler images. This tests were images of a cone on the ground and not in the sky, so the front and back stop feedback must be removed first. Then you can proceed to ”rack and stack” the return signal to get the range/doppler image. A pulse burst signal consists of numerous micropulses that form a macropulse in the shape of a pulse tone. The duration of the macropulse is 16 microseconds. The duration of the micropulse is 2 nanoseconds. There are 40 nanoseconds between micropulses. The rack and stack refers to chopping up the return signal to align the micropulses. So this is done in 40 nanosecond intervals which corresponds to a 6 meter range. Then for each range you perform a Fourier transform to get the doppler information. The doppler frequency range turns out to be +/- 12.5 MHz. That is why the local oscillator on the laser is set to 6.25 MHz. It is to try to center the micropulse at that frequency but the location is not guaranteed.

3

Simulations and Experimental Results

For the experiments, we did not use the detected data because we had yet to model the speckle. Instead we used the forward model to generate the ideal result and tried to estimate the parameters from there. It is useful to know if the ideal model can be optimized. In general, we found that it is possible to optimize 2 To compute the value of the range/doppler image at a specific range/cross-range coordinate (r,c), we integrate the irradiance over a region on the cone surface that corresponds to (r,c). It is this region that depends on the parameters.

4

the ideal model using the basic Matlab fminsearch. This required using enough sample points on the cone. Too few points lead to terrible results. For example, we saw huge differences in 1 cm vs 5 mm arc length sampling. We also made graphs to visualize the function to optimize by varying two parameters at a time and the 5 mm graphs were considerably smoother. We also expected and observed that when the cone is coming directly at you, it is not possible to estimate the spin rate because you are not getting any spin doppler readings. We also started implementing the use of better optimizers such as LBFGSB and an optimizer that allows for parameter constraints and confidence regions. We finished a version allowing for parameter constraints3 , but have not computed confidence regions.4

3.1 3.1.1

Guerin’s Experiments Experimental Settings

code Cone distance Aspect Arc length (S)Spin rate (D)Base Diameter (L)Cone length Outgoing wavelength Range/doppler graph range doppler 3.1.2

mfiles/flttd/spin test.m 19,000 meters 90 (sideways), 0 (head-on), or 45 degrees 0.005 meters 2 rps 0.75 meters 2 meters 10.15 × 10−6 (80x512 points) [-3, 3] [−12.5 × 10−6 , 12.5 × 106 ]

Experimental Results

E0 - initial parameter estimates EF - final parameter estimates Aspect 90 deg E0 S=3 D=1.75 L=3 EF S=2 D=0.75 L=2 Aspect 45 deg E0 S=3 D=1.75 L=3 EF S=2 D=0.75 L=2 Aspect 0 deg E0 S=3 D=1.75 L=3 EF S=4.96 D=0.75 L=2 Aspect 45 deg E0 S=3 D=1.75 L=1 EF S=2 D=0.75 L=2 Aspect 45 deg (only unexpected result) 3 in

spin cone estimate.m, but see also spin cone estimate test.m briefly considering this idea, it became less clear that confidence regions generated by the non-linear least squares estimation procedure in spin cone estimate.m would provide any meaningful information. Most likely, rather than tell us about our final parameter estimates, such statistics would indicate levels of confidence in the estimated parameter steps at each successive iteration. 4 Upon

5

E0 S=1 D=0.25 L=1 EF S=1.79 D=0.709 L=2.64

3.2

DeMeo’s Experiments

Using our forward model, simulated range/cross-range data are generated for an ideal cone having parameters values θt = (Lt , St , Dt )0 . (Examples of the parameter values we used are given below.) The resulting image, D = D(θt ), is stored in an 80 × 512. (In our Matlab code the array storing this “real data” is named rd). A new range/cross-range image is produced with the forward model using some arbitrary initial guesses, denoted θ0 , for the length, spin-rate, and base diameter parameters. This new image, G(θ0 ), is stored in a second 80 × 512 array (called g in our Matlab programs), and is evaluated based on how close it is, in Euclidean distance, to the given data D. More precisely, we consider the norm difference, kD − G(θ0 )k, and proceed with an iterative scheme for estimating parameter values θˆ which satisfy ˆ = min kD − G(θ)k kD − G(θ)k θ

Starting with θ0 , D, and G(θ0 ), to arrive at the next estimate, θ1 , we proceed as follows: 1. Approximate G(θ1 ) by a first order Taylor expansion about θ0 , G(θ1 ) ≈ G(θ0 ) + ∇G(θ0 ) (θ1 − θ0 ), where ∇G(θ0 ) is the gradient of G evaluated at θ0 . Then kD − G(θ1 )k is approximately kD − G(θ0 ) − ∇G(θ0 ) (θ1 − θ0 )k = kZ − ∇G(θ0 )ηk, where η = θ1 − θ0 . (We approximate the gradient of G by computing its finite differences using the forward model.) 2. Solve the OLS problem kZ − ∇G(θ0 )ˆ η k = min kZ − ∇G(θ0 )ηk. η

3. The new parameter estimate is θ1 = θ0 + ηˆ. 3.2.1

Experimental Settings

True Parameter length spin-rate base-diameter

Values5 L = 2m, S = 2Hz, D = 0.75m.

Initial Parameter length spin-rate base-diameter

Estimates (Guesses) θ1 = 1m, θ2 = 1Hz, θ3 = 1m.

Other Notable Parameter Settings Arc-length for the sampling grid on the cone: Differential used for computing gradients: Distance from laser to center of base of cone:

0.01m (arc = 0.01;). 5% (delta = .05;). 19km (dist = 19000).

These and other parameter values were held fixed for the three experiments below. Their Matlab definitions were set (in spin cone estimate test.m) as follows: 5 These are the parameter values used to generate the simulated raw data. Success is achieved if our parameter estimates converge to these values.

6

arc = 0.01; % cone sampling rate (arc length, circle spacing) L = 2; % length of the cone S = 2; % spin rate D = 0.75; % base diameter Lambda0 = 10.15*10^-6;%outgoing wavelength min_range = -3; max_range = 3; min_doppler = -12.5*10^6; max_doppler = 12.5*10^6; theta = [1; 1; 1]; % Initial guess: length, spin rate, base diameter delta = .05; % finite difference step-size (in percent of param value) TOL = 10^-8; % optimizer error tolerance for convergence LN = 1; % contrain parameters? (0,1) (no,yes) LB = [.25; .25; .25]; % lower bounds on parameter values UB = [ 5; 5; 5]; % upper bounds on parameter values MAXITERS=50; % maximum number of iterations

3.2.2

Experimental Results

A range/cross-range image produced with the forward model using the “true” parameter values (L=2, S=2, D=0.75) is taken as “real data,” and stored in the array rd. A second range/cross-range image is produced using the initial guesses for length, spin-rate, and base diameter, and the result is stored in the array g. The iterative estimation procedure was performed (in the Matlab script spin cone estimate test.m), yielding the results given below. The results correspond to three different scenarios, distinguished only by cone orientation. All other parameters we left unchanged. In particular, the length, spin-rate, and base diameter parameter values of the “true” cone were fixed at 2m, 2Hz, and 0.75m, respectively. 1. Aspect vector: A = [0; -dist; 0] (90 degrees) After 50 iterations, estimates were: L = 2.00483, S = 2.0008, D = 0.750095 with residual: |rd-g| = 0.000363976. Here rd is the range/cross-range image for the simulated data and g is that produced by the forward model using the cone parameter estimates, (L, S, D) = (2.00519, 1.99915, 0.750233). Convergence of the iterative procedure is illustrated in Fig. 2. 2. Aspect vector: A = [0; -dist; 0] (“head-on”) After 50 iterations, estimates were: L = 2.00024, S = 1, D = 0.750043 with residual: |rd-g| = 0.000017615 Convergence of the iterative procedure is illustrated in Fig. 3. 3. Aspect vector: A = [angle, -angle, 0] (45 degrees) where angle = sqrt(dist^2/2), and dist set at 19km. After 50 iterations, estimates were: L = 2.00069, S = 1.99767, D = 0.750129 with residual: |rd-g| = 0.000439862 Convergence of the iterative procedure is illustrated in Fig. 4.

7

Head On Cone Orientation

3

L = 2m

90 degree angle orientation

3

2 L = 2m

2 1 3

1 1

0

10

20

30

40

50

60

1

2 0

10

20

30

40

0.9 0.8 0.7 0.05

0

50

60

0.2 10

20

30

40

50

10

20

30

40

50

60

0

3

45 degree angle orientation

0

10

20

30

40

50

60

0

10

20

30

40

50

60

0

10

20

30

40

50

60

0

10

20

30

40

50

60

0.9 0.8 0.7 0.1 0.05 0

60 S = 2Hz

0

10

20

30

40

50

60

D = 0.75m

0

10

20

30

40

50

60

0

10

20

30

40

50

60

Figure 3: Convergence of the non-linear least squares estimation of cone parameters for head-on (0 degree) view. Algorithm converges to the correct parameter values for length and base diameter. However, spin rate estimate does not change from the original guess of 1Hz, which is what we expect for this orientation.

2 1

50

60

2 1

40

0.1

Figure 2: Convergence of the non-linear least squares estimation of cone parameters for perpendicular (90 degree) view. Algorithm converges to the correct parameter values for length, spin rate, and base diameter.

3

30

norm of residuals

error = 0.00036

0

20

0.9 0.8 0.7

D = 0.75m

0

10

1

S = 2Hz

1

0

Figure 4: Convergence of the non-linear least squares estimation of cone parameters for diagonal (45 degree) view. Algorithm converges to the correct parameter values for length, spin rate, and base diameter.

8