The matrix formulae of beam higher-order moments through an optical system in atmospheric turbulence are derived by using the Wigner distribution function (WDF). It is shown that the fourth-order moment matrix depends on two turbulence parameters Tμν and μν, the third-order and the second-order moment matrixes depend on only Tμν, and the first-order moment matrix is independent of turbulence. The results obtained in this paper are general, for example, the matrix formulae of beam higher-order moments through Kolmogorov turbulence, non-Kolmogorov turbulence, optical system or free space can be given as special cases. The advantages of matrix formulation are that the expressions for higher-order moments are brief, and the physical meanings of higher-order moments are clear. The results obtained in this paper are of considerable theoretical and practical interest.