Rather than the frequency domain, let's look at this in the time domain and particularly, the characteristic equation associated with a linear homogeneous 2nd order differential equation for some system:
\$r^2 + 2 \zeta \omega_n r + \omega^2_n = 0\$.
If the roots of the characteristic equation are real (which is the case if \$\zeta \ge 1\$), the general solution is the sum of real exponentials:
\$Ae^{\sigma_1 t} + Be^{\sigma_2t} \$
where
\$\sigma_1 = -\zeta \omega_n + \sqrt{(\zeta ^2 - 1)\omega^2_n} \$
\$\sigma_2 = -\zeta \omega_n - \sqrt{(\zeta ^2 - 1)\omega^2_n} \$
Since these are real exponentials, there is no oscillation in these solutions.
If the roots are complex conjugates (which is the case if \$\zeta < 1\$), the general solution is the sum of complex exponentials:
\$e^{\sigma t}(Ae^{j\omega t} + Be^{-j\omega t})\$
where
\$\sigma = -\zeta \omega_n\$
\$\omega = \sqrt{(1 - \zeta ^2)\omega^2_n}\$
This solution is a sinusoid with angular frequency \$\omega\$ multiplied by a real exponential. We say the system has a "natural frequency" of \$\omega\$ for a reason that I think is obvious.
Finally, setting \$\zeta = 0\$ (an undamped system) , this solution becomes:
\$Ae^{j\omega_n t} + Be^{-j\omega_n t}\$
which is just a sinusoid of angular frequency \$\omega_n\$.
In summary, a system may or may not have an associated natural frequency. Only systems with \$\zeta < 1\$ have a natural frequency \$\omega\$ and only in the case that \$\zeta = 0\$ will the natural frequency \$\omega = \omega_n\$, the undamped natural frequency.
Two possible methods of several:
A 2nd order CLTF model of a higher order system (viz. the \$\small\zeta\$ and \$ \omega _n\$ values), can be derived from the open-loop gain cross-over frequency, and phase margin, thus:
$$\small\left(\frac{\omega _c}{\omega _n}\right)^2 \approx \sqrt{(4\zeta ^4 +1)}\:-2\zeta^2 $$
and
$$\small\zeta \approx 0.01\times PM $$
where \$\small PM\$ is the phase margin in degrees, and \$\small \omega _c\$ is the frequency at which the open loop gain crosses \$\small 0\:dB\$ (i.e. the frequency at which the phase margin is defined).
This assumes that \$ \omega _c\$ (and consequently the phase margin) exists.
- Dominant poles: the dominant poles are located at \$\small s=0\$ and \$\small s=-5\$, since the other two poles are at least three times further from the origin than the \$\small (s+5)\$ pole. Hence ignoring the less dominant poles gives the OLTF:
$$\small G_c (s)=\frac{K(s+7)}{s(s+5)} $$
with CLTF:
$$\small G(s)=\frac{K(s+7)}{s(s+5)+K(s+7)}=\frac{K(s+7)}{s^2 +(K+5)s+7K} $$
Best Answer
You can estimate the damping ratio value from the Bode plot, of course. How good that estimate will be depends on the accuracy of your readings on the graph. (I redid the computation after carefully centering the reference lines and I got a slightly different result). But first and foremost you need to know what kind of transfer function you have plotted.
the pixel distance between 10k and 100k is 759
the pixel distance between 10k and f3dB is 211
the pixel distance between 10k and fn is 271
we can easily find
fn (kHz) = 10^(1+ 271/759) = 10^1.35765 = 22.75 (kHz)
so wn = 2 pi 22.75 kHz.
once we have traced the -3dB reference (just count the pixels in a 10 dB interval and move the horizonal reference 3/10 of that many pixel down from the low frequency plateau - or use the simulator utilities).
we can find the 3dB cutoff frequency. Again by counting pixels we find
f3dB (kHz) = 10^(1+ 211/759) = 10^(1.277998) = 18.967 (kHz)
we thus get w3dB = 2 pi 18.97 kHz. the fact that w3dB < wn tells us that this is a peakless transform (might be overdamped, critically damped or underdamped)
w3dB/wn = 18.967/22.75 = 0.8337...
(my previous lazy eyeball estimate was 20/22.53 = 0.8876. Kinda prove my point that the accuracy of the results relies on how carefully you get the values)
Note that the fact that w3dB > 0.644 wn tells us that we are in the underdamped case where, since w3dB < wn, there is no peak in the transfer function.
Relative position of 3dB cutoff frequency with respect to natural undamped frequency wn allows us to tell whether the second order transfer function without zeroes belongs to a: overdamped (w3dB < 0.644 wn), critically damped (w3dB = 0.644 wn) or underdamped (w3dB > 0.644 wn) system and even if the magnitude of the transfer function has a peak (w3dB > wn) or not ( 0.644 wn < w3dB < wn). The above graph shows the TF for two distinct real poles (whose partial contributes are shown) for zeta = 1.75, that is for an overdamped system, and dashed, the curve for a critically damped system with two identical real poles)
We have to expect a value of zeta greater than 1/sqrt(2) = 0.707...
Solving for w3dB/wn = 0.8337 we find
zeta = 0.828
(and wd = 0.56 wn)
This new, hopefully more accurate estimate for zeta is even farther away from the frontier between the no-peak vs peak regions of the underdamped zone. If the assumption we made about the transfer function is correct, then we are in the underdamped region (zeta < 1) where there is no peak (zeta > 0.707).
The boundary bewteen the "No peak" and "Resonance peak" regions in the underdamped zone is at zeta = 1/sqrt(2) = 0.707...
If we round the 3dB frequency at 19 kHz, with w3dB/wn = 10/22.75 = we get zeta = 0.827. In any case, that's a far cry from the critically damped case of zeta=1 and not even particularly close to the border between the "there is a peak" and "there is not a peak" set by zeta = 1/sqrt2 = 0.707... (corresponding to w3dB=0.644 wn). The poles are complex conjugates with an angle of 56° wrt to the imaginary axis and 34° with respect to the negative real axis (the peak no-peak border requires an angle of 45°)
I suggest you redo the computation after having extracted your values from the graph. Better yet, do it on a known function to assess the accuracy of the method.
See also my last answer to the question Finding resonant frequency or damping ratio from Bode Plot