In this paper, by using the real representations of quaternion matrices, the particular structure of the real representations of quaternion matrices, the Kronecker product of matrices and the Moore–Penrose generalized inverse, we obtain the expressions of the minimal norm least squares solution, the pure imaginary least squares solution, and the real least squares solution for the quaternion matrix equation AXB+CXD=E, respectively. Our resulting formulas only involve real matrices, and therefore are simpler than those reported in Yuan (2014). The corresponding algorithms only perform real arithmetic which also consider the particular structure of the real representations of quaternion matrices, and therefore are very efficient and portable. Numerical examples are provided to illustrate the efficiency of our algorithms.