1 Introduction

Disk-resolved images can be obtained only of a limited number of the small atmosphereless bodies of our solar system. This is why disk-integrated photometry, i.e., measurements of total brightnesses, is and will remain a major source of information on these objects. The total brightness of an object as a function of time (as the viewing/ illumination geometry changes) is called its lightcurve; the parameters to be determined from lightcurve observations are the object’s shape, its rotational state, and the scattering properties of its surface.

The problem of lightcurve inversion has been studied for nearly a century; however, restricted choices of the observing geometry and scattering law in the early analytical studies have caused erroneously negative preconceptions about the uniqueness and stability of the solution. Partly because of this, and partly because of the relative scarcity of accurate lightcurve observations, the shape models used in lightcurve inversion have until recently been extremely simple – typically modifications of the triaxial ellipsoid (see, e.g., [6]). The first general and (semi)analytical solution was offered in [2]. This study proved the uniqueness of the solution for convex bodies (when various observing geometries are available), but the apparent ill-posedness of the problem requiring heavy regularization remained an obstacle in practice.

Some probe and radar images of asteroids have been obtained during the last decade (see, e.g., [9], [1], [10] and references therein). The most conspicuous feature is the existence of substantial nonconvexities: there are globally significant valleys and indentations as well as large local craters on the surfaces. The number and accuracy of ground-based lightcurve measurements have also increased significantly. The new observations indicate that one should study a) how useful and descriptive convex inversion is in general (even strongly nonconvex bodies can often be reasonably well represented by their convex hulls or similar shapes), and b) whether nonconvex inversion is feasible.

In the following, we describe new practical methods for posing and solving the inverse problem, and we present some applications to astronomical data. The scattering law and the size scale of the object can be chosen quite freely, so our techniques can be applied to various other problems in, e.g., general remote sensing as well. A detailed description of the approach is to appear elsewhere ([3], [4]), so we only outline the basic concepts here.

2 Convex inversion

Once it is known that a surface patch dσ is both visible and illuminated, its contribution dL to the total brightness L is given by (omitting trivial scale factors such as the squares of distances)

where S and ϖ are the scattering law (in a simple form here; more arguments can naturally be included) and albedo; µ=E·n and µ0=E ₀·n, where E and E ₀ are, respectively, unit vectors towards the observer (Earth) and the Sun, and n is the surface unit normal. Lambert law, for example, is S _L=µµ ₀, while Lommel-Seeliger law is S _LS=S _L/(µ+µ ₀).

The convex inverse problem can be cast in the form

where L is the vector of the observed brightnesses, related through the matrix A to the vector g that contains the shape parameters to be solved.

2.1 Polyhedron model

Let g describe the areas of the facets of the polyhedron approximating the surface; their surface normals are fixed using, e.g., a suitable grid. The matrix elements A_ij are then obtained by applying Eq. (1) to each facet j (representing dσ) and observing geometry i in a frame of reference fixed to the object.

The standard solution of Eq. (2) by minimizing the square norm

using least-squares normal equations or singular value decomposition would usually produce negative g_j s as the system is ill-conditioned. On the other hand, any standard regularization method would produce a smoothing element. Apart from the convexity constraint, the only requirement is that the facet areas be positive; no correlation between neighbours must be assumed. The easiest way to guarantee positivity is to represent each g_j exponentially, the optimization parameter being now the exponent a_j : g_j =exp(a_j ). It is straightforward to verify that the solution vector a minimizing χ² is unique, and that the global minimum is the only local one. Thus any nonlinear optimization method will ultimately find the correct solution; we use conjugate gradients [8] as the number of facets is large (typically≥1000).

Once the areas of the facets are known, the vertices of the facets can be obtained by representing the so-called Minkowski problem as a nonlinear optimization problem ([2], [3]) that can be solved with standard methods. An important fact is that the Minkowski problem is very stable; this stability ensures that the shape solution remains close to the ‘correct’ one (i.e., one that approximates the original nonconvex shape well) even if the scattering model is not very accurate.

For the surface to be convex, we must have ∑ _j n _j g_j =0. If this residual vector is nonzero, a convex solution can be found by augmenting the matrix A by three new rows corresponding to the three linear constraints (and their weight factors) while adding three zeros to L; with similar regularization, we also get a (feasible but not unique) map of the albedo distribution over the surface.

2.2 Function series model

A convex surface can also be uniquely described by its curvature function G (the inverse of Gaussian curvature) as a positive definite function of the direction of the surface normal. The positivity constraint is automatically fulfilled if we write G as

where (ϑ, ψ) are the spherical coordinates of the surface normal and $Y_l^m$ are spherical harmonics. If albedo variegation is absorbed into G, the observed brightness is

where the integration region A ₊ on the unit sphere of normal directions is the part on which µ, µ ₀≥0; S is the scattering law, and dσ is the surface element of the sphere.

The numerical integration of Eq. (5) is easy to perform on a densely triangulated unit sphere, using its facets to approximate dσ. The triangulation returns Eq. (5) to Eq. (2) and Eq. (1), with

where Δσ_j is the area of the sphere facet corresponding to (ϑ_j, ψ_j ). Since we are minimizing a nonlinear least-squares function Eq. (3) and the number of the coefficients a_lm to be solved for is not large (typically from, say, 40 to 100), it is advantageous to use the Levenberg-Marquardt optimization method [8] that usually converges robustly.

 

Fig. 1. The shape model and a lightcurve of asteroid 6489 Golevka at two observing geometries [two movies (690 kB, 690 kB)].

Download Full Size | PDF

In Fig. 1 we present representative frames of two video clips showing the convex model of asteroid 6489 Golevka rotating in the viewing/illumination geometries corresponding to two observed lightcurves. Also shown is the time evolution of the model lightcurve against the measured brightnesses. α is the phase angle, θ ₀ the polar angle of illumination (Sun) measured from the asteroid’s north pole, and θ the corresponding angle of the observer (Earth).

3 Scattering properties

Any parameters other than those of shape (or albedo) can be changed to free ones and optimized together with the rest without changing the method; this is especially easy with the function series model. The rotation parameters (period and axis direction) are readily included in rotation matrices operating on direction vectors. Another important set of parameters consists of the coefficients of the functions in the scattering law.

Accurate, universal few-parameter descriptions of the scattering properties of the surface material are hard to come by. The existing models, such as those by Hapke or Lumme and Bowell, are notorious for producing ambiguous and unrealistic parameter values in inverse problems. A typical example of a scattering phenomenon yet to be properly modelled is the opposition effect, i.e., the brightening near zero phase angle, caused by coherent backscattering and shadowing. This phenomenon alone often forces us to use relative brightness measurements instead of absolute ones: if the effect of solar phase angle on surface brightness cannot be modelled accurately, fitting absolute brightnesses is useless. This problem is further exacerbated by the unfortunate fact that the absolute brightnesses of observations are typically not known accurately enough. In order to obtain proper line-over-points fits, we need to rewrite the objective function in the form

where < A ⁽ⁱ⁾ g > is the mean brightness of the ith model lightcurve fitting the ith observed lightcurve (through which α, due to the slow orbital motion, usually remains virtually constant) whose mean brightness is L̄ ⁽ⁱ⁾. Thus both the observed and the model lightcurves are renormalized to mean brightnesses of unity.

We have found that the solution based on relative brightnesses is for all practical purposes equivalent to the one obtained from absolute photometry. This is quite understandable: it is primarily the shapes of the lightcurves that are strongly connected with the pole, the period, and the shape of the asteroid. The absolute brightnesses are principally connected with the scattering properties; thus it is actually useful to employ Eq. (7) to decouple one set of parameters from the rest as much as possible. To make this decoupling manifest, we employ the scattering model

which combines single (Lommel-Seeliger term S _LS) and multiple scattering (Lambert term S _L) with a weight factor c for the latter. For the sake of convenient inversion, the phase function f(α) is taken to multiply the sum of the single and multiple scattering terms. Such a simplification can be justified by the well-known ambiguity in what should be called single or multiple scattering in a particulate medium consisting of inhomogeneous particles. In this formulation f(α) can be determined afterwards from a set of scale factors for each lightcurve (see Fig. 2) while it does not have to be known when solving for the other parameters.

 

Fig. 2. The phase function f(α) for asteroid 433 Eros.

Download Full Size | PDF

The philosophy behind Eq. (8) is that, for lightcurve inversion, the scattering law must be simple: too many parameters and possibilities cause instability and unrealistic results. Also, we want to avoid detailed physical parameters and simply try to describe the general photometric properties of the surface. In the same vein, we keep f(α) simple by fitting an exponential and linear model to the scale factors. This has proven to be a versatile choice when modelling the intensity of scattered light as a function of phase angle (cf. [5]), although any other function can be used without altering the pole/shape result in any way. Since f(α) can be multiplied by any constant (when the absolute scales of the size and darkness of the asteroid need not be fixed), we can normalize the linear part to unity at opposition after the fit to be able to compare the shapes of different phase functions. Thus we report f(α) in the three-parameter form

where A ₀ and D are the amplitude and scale length of the opposition effect, and k is the overall slope of the phase curve. In Fig. 2, c=0.1,A ₀=1.0,D=4.5°, and k=-0.009.

Note that the present empirical modelling does not require hypotheses about the physics behind the opposition effect. In fact, we do not even have to know whether the opposition effect exists in the first place: if Eq. (9) does not fit the scale factors well, we can use any other function or just plot the points as such. We have found that Hapke, Lumme-Bowell, and the above scattering model all give very closely the same results for the shape and the rotational state; examples of this are shown in [4]. The drawbacks of an empirical model are that it can seldom be used for extrapolation, and the points to be fitted have to be dense for proper interpolation. Physical models usually constrain the shape of the phase curve more efficiently with fewer points.

4 Nonconvex inversion

The main complications in nonconvex modelling are that all uniqueness theorems are lost, and the parameter space is usually plagued by local minima. There are not very many methods available; after some experimenting, we have found it best to use a short functional series describing the locations of the vertices of a triangulated surface. One possibility is to use an expression such as Eq. (4) for the radii of the vertices in given directions (cf. [7]).

The trial lightcurves are now computed using Eq. (1) and a fast ray-tracing algorithm; the resulting χ² is minimized with Levenberg-Marquardt or genetic algorithms. The series is typically truncated at an early order; the contributions of detailed nonconvexities are usually drowned in the noise.

The weak effect of macroscopic nonconvex features even at large solar phase angles is somewhat surprising. Indeed, we have not yet found lightcurve data that could not be explained well with a convex model, even when produced by shapes known to have large concave features. This is a result of the smoothing property of lightcurves and the fact that the inversion result is a ‘photometric’ rather than geometric convex hull of the body. Nonconvex models do not really reach lower χ²s (in fact, their fits are usually worse due to the numerous ‘blind alleys’ in parameter space). Minkowski stability no longer applies to them: thus, even in theory, reliable general nonconvex inversion (without a priori fixed types of nonconvexities) requires highly accurate observations, very favourable observation geometries, and an accurate scattering model. Luckily, large flat areas on the convex solution already indicate the presence and locations of major nonconvex features.

In short: pole, period, and the photometric convex hull of the object overwhelmingly dominate lightcurve morphology. Concavities and scattering properties are, in a way, comparable to noise: even though changing them does change the lightcurves somewhat in the direct problem, the solution of the inverse problem does not change much since there are no notably different convex shapes that could model the lightcurves better. The well-posedness and stability of the inversion procedure is thus a result of the convex formulation.