For most geophysical inverse problems, one solves for an Earth model that minimizes the misfit between the observed data and observational predications determined from forward modeling. However, the solutions to such problems are almost always nonunique, as some of the model parameters are poorly resolved by the data. We have developed a novel form of evolutionary programming (EP) that solves directly for a broad cluster of models that adequately fit the data. We obtain improved performance of the EP algorithm by rotating the model parameter axes at each generation so that the axes are aligned perpendicular and parallel to long valleys in the misfit surface. The rotated axes are defined by the eigenvectors of the covariance matrix of the set of models; the mutation rate in each new direction varies with the square root of the corresponding eigenvalue. This modified EP method finds a wide range of models satisfying the data, allowing one to determine the uncertainty in the model parameters. We illustrate this EP method by applying it to the problem of determining velocity structure from regional travel-time data. Application of this method to a travel-time data set generated from a model that includes a low-velocity zone (LVZ) illustrates that it successfully finds a wide range of models within the broad minimum that includes the initial solution. We apply the method to regional travel-time data recorded at the Kyrgyz broadband network (KNET) and find that velocities are well resolved, except at depths of 25 to 50 km.