This note studies the iterative solution to the coupled quaternion matrix equations $[\sum_{i=1}^{p}\mathbb{T}_{1i}(X_i)$, $\sum_{i=1}^{p}\mathbb{T}_{2i}(X_i)$, $\cdots$, $\sum_{i=1}^{p}\mathbb{T}_{pi}(X_i)]$= $[M_1, M_2,\cdots ,M_p]$, where $\mathbb{T}_{si}, s=1,2,\cdots,p;$ is a linear operator from $Q^{m_i\times n_i}$ onto $Q^{p_s\times q_s}$, $M_s \in Q^{p_s\times q_s}, s=1,2,\cdots,p. i=1,2,\cdots,p$. by making use of a generalization of the classical complex conjugate graduate iterative algorithm. Based on the proposed iterative algorithm, the existence conditions of solution to the above coupled quaternion matrix equations can be determined. When the considered coupled quaternion matrix equations is consistent, it is proven by using a real inner product in quaternion space as a tool that a solution can be obtained within finite iterative steps for any initial quaternion matrices $[X_1(0),\cdots, X_p(0)]$ in the absence of round-off errors and the least Frobenius norm solution can be derived by choosing a special kind of initial quaternion matrices. Furthermore, the optimal approximation solution to a given quaternion matrix can be derived. Finally, a numerical example is given to show the efficiency of the presented iterative method. Huang Wenxue 8031