Conventional electromagnetic and elastic methods are inapplicable to quantitative detection of voids under the floor of sluice operating underwater. Motivated by the difficulty in detecting internal defects in complex underwater structures, a novel numerical method for modeling elastic wave propagation in heterogeneous media has been proposed by deriving the k-space form of the first-order stress-velocity equation. The spatial-temporal rotated staggered grid and the perfectly matched layer are introduced to improve accuracy and efficiency. The k-space method is used to simulate impact elastic waves induced by voids under the floor of sluice operating underwater. The time-domain waveform and spectral response of impact elastic waves are investigated by comparing with the realistic data. Our results demonstrate that compared with the finite difference, finite element, and pseudo-spectral methods, the proposed k-space method can comprehensively simulate the reflected and scattered elastic waves as well as fluid-solid coupling effects of two-phase multi-layer structures with significantly improved accuracy and efficiency. The void causes a significant larger and more slower decay amplitude and a longer lasting time-domain waveform. The void is characterized by the third peak frequency in the spectrum of the impact elastic waveform. The proposed k-space method has potential impact on the efficient and accurate elastic-wave detection of internal defects in complex underwater structures such as sluices, dikes, dams and bridges.