The generalized statistical dynamical theory of x-ray scattering by imperfect single crystals with randomly distributed Coulomb-type defects has been extended to characterize structure imperfections in the real multilayers of arbitrary thickness in Bragg diffraction geometry. The recurrence relations for the coherent amplitude reflection and transmission coefficients of such multilayers, which consist of any number of layers with constant strain and randomly distributed defects in each one, have been derived within the concept of the dynamical wave field, i.e., the so-called Ewald-Bethe-Laue approach, with rigorous accounting for boundary conditions at layer interfaces. The analytical expression for the differential dynamical diffuse component of the reflection coefficient of an imperfect multilayer system has been obtained as well. This expression establishes direct connection between the distribution of the diffuse scattering intensity in a momentum space and statistical characteristics of defects in each layer. In addition, the integrals from differential diffuse scattering intensity over the Ewald sphere and over vertical divergence have been found, which correspond to diffuse components in measurements of rocking curves and reciprocal space maps, respectively. The developed theory provides the dynamical description for the one- and two-dimensional angular distributions of the mutually consistent coherent and diffuse components of x-ray scattering intensities, which are measured by the high-resolution double- and triple-crystal diffractometers, respectively, from imperfect films, multilayer structures, superlattices, etc. Examples of simulated rocking curves for the imperfect superlattice with defects of several types and reciprocal space map for the ion-implanted sample of yttrium iron garnet film with defects are given and discussed.