Full waveform inversion reconstructs subsurface structures by matching the synthetic waveform to the observed waveform. Inaccuracy of the source wavelets can, thus, easily lead to an inaccurate model. Simultaneously updating source wavelets and model parameters is a conventionally used strategy. However, when the initial model is very far from the true model, cycle skipping exists, and estimating a reliable source wavelet is very difficult. We propose a combinatory inversion workflow based on seismic events. We apply a Gaussian time window around the first break and gradually increase its width to include more seismic events. The influence of inaccurate source wavelets is alleviated by applying a Gaussian time window around the first break to evaluate the normalized cross-correlation-based objective function. There are inevitable small model artifacts caused by inter-event interactions when calculating cross-correlations. As a result, we switch to the optimal transport function to clean the model and update the source wavelets simultaneously. The combinatory strategy has been applied to models with different types of geological structures. Starting from a crude initial model, we recovered a high-resolution and high-fidelity model and the source wavelets in two synthetic experiments. Finally, we apply our inversion strategy to a real-land seismic dataset in Southeast China and obtain a higher-resolution velocity model. By comparing an inversion velocity profile with well log information and the recorded data with the simulated data, we conclude that our inversion results for the field data are accurate and this new strategy is effective.