While restoration methods have been extensively studied in ultrasound (US) imaging, only few recent works have focused on modeling and understanding the blur from a physical point of view even in simple configurations such as lossless homogeneous media. Despite a highly non-stationary blur due to diffraction effects, many techniques rely on simplistic approximations based on shift-invariant models or sectional methods and their efficiency has not been demonstrated for plane-wave (PW) and diverging-wave (DW) imaging. In this work, we first propose a physical model of non-stationary blur in the context of PW and DW imaging. The blur operation is expressed as a composition of a US propagation operator and a delay-and-sum operator, each of which derived in previous works. We show that such a composition leads to a standard model of non-stationary blur as a Fredholm integral of the first kind. Secondly, we describe an approximation of the blur in the discrete domain based on the above decomposition coupled with an appropriate discretization of the latent space of the element-raw data. We show theoretically and empirically that its evaluation, using such an approximation, scales linearly instead of quadratically with respect to the grid size, better than shift-invariant approaches. Through simulations and in vivo experimental data, we demonstrate that using the proposed model in the context of maximum-a-posteriori image restoration results in a higher image quality than using state-of-the- art shift-invariant models.