Migration attempts to produce an image of the subsurface by reversing the propagation effects in seismic data. Although in principle this requires the inverse of a modelling operator, in practice the adjoint of the modelling operator is used instead. In cases where the data are subject to significant aliasing, truncation, noise, or are incomplete, the adjoint of a modelling operator is not a good approximation to the inverse, and this degrades the resolution of the final migrated image. An improved approximation to the inverse operator can be obtained using a least-squares approach. In this work we present a least-squares formulation of reverse-time migration scheme that is based on an explicit matrix representation of generalised diffraction-stack migration. Our implementation uses a modified source wavelet to perform the forward and inverse steps at each iteration, and is regularised using a roughness penalty constraint.