Vertical seismic profile (VSP) data are usually acquired with three-component geophones of unknown azimuthal orientation. The geophone orientation must be estimated from the recorded data as a prerequisite to processing such as P- and S-wave separation, calculation of wave-incident directions, and 3D migration. We compare and combine two methods for estimating azimuthal orientation by least-squares fitting over a large number of shots. Combining the two methods can be done in an automated manner, which provides more accurate estimates of the geophone orientations than previous methods. In the polarization-plane method, we calculate the polarization plane of the first P-wave arrival. Then we subtract the source azimuth to determine the geophone orientation, independently for each geophone, with an angular uncertainty of [Formula: see text], and with no accumulated errors. In the relative-angle method, we obtain relative angles between adjacent geophone pairs using trace crosscorrelations, and operate on all coherent signals (even noise). Swapped geophone components can be detected automatically using the polarization-plane method. The main limitation of these (and all other known) methods is that uncertainties associated with path refraction are not estimated, unless some geophones have a priori known orientations, or we have a known earth model to correct for refraction.