Abstract: We investigate the anisotropic hybrid surface optical (SO) phonon-plasmon dispersion relations in monolayer and double-layer phosphorene systems located on the polar substrates, such as SiO2, h-BN, and Al2O3. We calculate these hybrid modes by using the dynamical dielectric function in the random phase approximation in which the electron-electron interaction and long-range electric field generated by the substrate SO phonons via Frohlich interaction are taken into account. In the long-wavelength limit, we obtain some analytical expressions for the hybrid SO phonon-plasmon dispersion relations which agree with those obtained from the loss function. Our results indicate a strong anisotropy in SO phonon-plasmon modes, which are stronger along the light-mass direction in our heterostructures. Furthermore, we find that the type of substrate has a significant effect on the dispersion relations of the coupled modes. Importantly, the hybrid excitations are apparently sensitive to the misalignment and separation between layers in double-layer phosphorene.