We discuss anharmonicity of the multi-octupole-phonon states in $^{208}$Pb based on a covariant density functional theory, by fully taking into account the interplay between the quadrupole and the octupole degrees of freedom. Our results indicate the existence of a large anharmonicity in the transition strengths, even though the excitation energies are similar to those in the harmonic limit. We also show that the quadrupole-shape fluctuation significantly enhances the fragmentation of the two-octupole-phonon states in $^{208}$Pb. Using those transition strengths as inputs to coupled channels calculations, we then discuss the fusion reaction of $^{16}$O+$^{208}$Pb at energies around the Coulomb barrier. We show that the anharmonicity of the octupole vibrational excitation considerably improves previous coupled-channels calculations in the harmonic oscillator limit, significantly reducing the height of the main peak in the fusion barrier distribution.