Conventional surface seismic data lack of usable low-frequency and long-offset signals, which leads to classical full-waveform inversion hardly retrieving long-wavelength contents of elastic models in the mid-to-deep parts. As an alternative approach, wave-equation reflection waveform inversion has drawn more attentions. However, almost all existing studies on reflection waveform inversion only use first-order optimization and thus have great room for improvement in convergence rate and accuracy. Here, in the framework of second-order optimization, we have derived reflection sensitivity kernels, functional gradient and Hessian operators with respect to background and perturbation components of the elastic models. We have demonstrated the debluring effect of Hessian operator for the functional gradient and revealed the working mechanism to provide optimal model updating and enhance inversion resolution. The synthetic example on an overthrust model shows that, compared with widely used conjugate-gradient method, approximate Hessian-based Gauss-Newton method significantly enhances the convergence rate and the capability of broadband model building. In the real data example of East China Sea, the proposed reflection waveform inversion outperforms widely used reflection travel-time tomography. It has improved velocity model building and thus supported reverse-time migration to delineate the fault systems and the basement in the Changjiang sag with higher resolution.