Only recently, OD, the deuterated isotopolog of hydroxyl, OH, has become accessible in the interstellar medium; spectral lines from both species have been observed in the supra-Terahertz and far infrared regime. Here we study rotational lines of OD and OH towards 13 Galactic high-mass star forming regions, with the aim to constrain the OD abundance and infer the deuterium fractionation of OH in their molecular envelopes. We used the Stratospheric Observatory for Infrared Astronomy (SOFIA) to observe the $^2\Pi_{3/2}$$J=5/2-3/2$ ground-state transition of OD at 1.3 THz ($215~\mu$m) and the rotationally excited OH line at 1.84 THz ($163~\mu$m). We also used published high-spectral-resolution SOFIA data of the OH ground-state transition at 2.51 THz ($119.3~\mu$m). Our results show that absorption from the $^2\Pi_{3/2}$ OD $J=5/2-3/2$ ground-state transition is prevalent in the dense clumps surrounding active sites of high-mass star formation. We performed detailed radiative transfer modelling to investigate the OD abundance profile in the inner envelope for a large fraction of our sample. Our modelling suggests that part of the absorption arises from the denser inner parts, while the bulk of it as seen with SOFIA originates in the outer, cold layers of the envelope for which our constraints on the molecular abundance suggest a strong enhancement in deuterium fractionation. We find a weak negative correlation between the OD abundance and the bolometric luminosity to mass ratio, an evolutionary indicator, suggesting a slow decrease of OD abundance with time. A comparison with HDO shows a similarly high deuterium fractionation for the two species in the cold envelopes, which is of the order of 0.48% for the best studied source, G34.26+0.15. Our results are consistent with chemical models that favour rapid exchange reactions to form OD in the dense cold gas.Read less <