We present a variant of a traditional genetic algorithm, known as a niching genetic algorithm (NGA), which is effective at multimodal function optimization. Such an algorithm is useful for geophysical inverse problems that contain more than one distinct solution. We illustrate the utility of an NGA via a multimodal seismological inverse problem: the inversion of teleseismic body waves for the source parameters of the MW 7.2 Kuril Islands event of 2 February 1996. We assume the source to be a pure double-couple event and so parametrize our models in terms of strike, dip, and slip, guaranteeing that two global minima exist, one of which represents the fault plane and the other the auxiliary plane. We use ray theory to compute the fundamental P and SH synthetic seismograms for a given source-receiver geometry; the synthetics for an arbitrary fault orientation are produced by taking linear combinations of these fundamentals, yielding a computationally fast forward problem. The NGA is successful at determining that two major solutions exist and at maintaining the solutions in a steady state. Several inferior solutions representing local minima of the objective function are found as well. The two best focal solutions we find for the Kuril Islands event are very nearly conjugate planes and are consistent with the focal planes reported by the Harvard CMT project. The solutions indicate thrust movement on a moderately dipping fault--a source typical of the convergent margin near the Kuril Islands.