In this paper, an iterative method is proposed to solve the least-squares problem of matrix equation AXB+CYD=E over unknown matrix pair [X, Y]. By this iterative method, for any initial matrix pair [$X_1,;Y_1$], a solution pair or the least-norm least-squares solution pair of which can be obtained within finite iterative steps in the absence of roundoff errors. In addition, we also consider the optimal approximation problem for the given matrix pair [$X_0,;Y_0$] in Frobenius norm. Given numerical examples show that the algorithm is efficient.