Abstract: The depth-domain seismic wavelet is the dynamic response function of the subsurface velocity field, and its waveform exhibits significant depth-variant characteristic. This characteristic makes the depth-domain seismic wavelet fail to satisfy the linear depth-invariant assumption in the conventional convolution model, resulting in the inability to apply conventional time-domain seismic wavelet extraction methods to the depth domain. To address this problem, this paper proposes a depth-domain seismic wavelet extraction method that integrates the depth-domain weighted stacking method and the iterative separate parameter estimation (ISPE) method. First, based on the depth-domain weighted stacking method, vertical continuity of the wavelet morphology is achieved through piecewise cubic Hermite interpolation. Second, combined with the non-stationary convolution model, depth-variant wavelet clusters and their Circulant matrices are constructed to directly synthesize depth-domain seismic records. Subsequently, within the ISPE framework, the depth-domain generalized seismic wavelet model is introduced, and the Pearson correlation coefficient is used as a quantitative metric for waveform similarity. This reduces the parameter dimension of seismic wavelet extraction while precisely controlling the iteration termination conditions, achieving accurate extraction of depth-domain seismic wavelets. Model tests and practical data applications show that even under conditions of strong noise interference, complex structural conditions, and limited seismic and well-logging data, this method can still extract reliable depth-domain seismic wavelets.