Write the TF as \$ \small G(s)=\large\frac{ab}{(s+a)(s+b)}\$, where \$a\$ and \$b\$ are real and distinct.
This factorises to: \$\small G(s)=\large\frac{ab}{(b-a)} \large(\frac{1}{s+a}-\frac{1}{s+b}\large)\$
From which the inverse LT, i.e. the unit impulse response, is:
\$g(t)= \large\frac{ab}{(b-a)} \normalsize (e^{-at}-e^{-bt})\$
Differentiating and setting to zero determines the time to peak value, \$t_P\$:
\$t_P =\large\frac{ln(\frac{b}{a})}{(b-a)}\$, with corresponding peak value: \$g(t_P)=be^{-b t_P}\$
Hence, measuring \$t_P\$ and \$g(t_P)\$ will yield the value of \$b\$; and substituting back into the expression for \$t_P\$ will yield \$a\$
Note that the unit impulse response can be obtained by differentiating the unit step response, but the differentiator must be carefully designed to avoid introducing excessive noise.
Typical overdamped step and impulse responses are shown below.
TL;DR: NO, you can't use the underdamped settling time formula to find out the settling time of an overdamped system. And you can't use it for a critically damped system either.
LONG FORM answer follows...
Critically damped case
For the critically damped case (\$\zeta=1\$), the step response is:
$$
v_{out}(t) = H_0 u(t) \lbrack 1 - (1+\omega_0 t) e^{-\omega_0 t} \rbrack
$$
If we define the settling time \$T_s\$ using the same "within 2% of final response" criteria, then:
$$
0.02 = (1+\omega_0 T_s) e^{-\omega_0 T_s}\\
$$
Solving numerically for \$\omega_0 T_s\$ (by simply using Excel's solver) we obtain:
$$
T_s \approx \frac{5.8335}{\omega_0}
$$
Overdamped case
For the overdamped case (\$\zeta>1\$), the step response is:
$$
v_{out}(t) = H_0 u(t) \left[ 1 - \frac{s_2}{s_2-s_1}e^{s_1 t} - \frac{s_1}{s_1-s_2}e^{s_2 t} \right]
$$
where \$s_1, s_2\$ are the real roots of the transfer function denominator:
$$
s_1 = -\zeta \omega_0 + \omega_0 \sqrt{\zeta^2-1} \\
s_2 = -\zeta \omega_0 - \omega_0 \sqrt{\zeta^2-1}
$$
For convenience we define:
$$
\begin{align}
\Delta &= \frac{s_2-s_1}{2} = - \omega_0 \sqrt{\zeta^2-1} \\
\Sigma &= \frac{s_1+s_2}{2} = - \zeta \omega_0 \\
K &= \frac{\Sigma}{\Delta} = \frac{\zeta}{\sqrt{\zeta^2-1}}
\end{align}
$$
So that:
$$
\begin{align}
s_1 &= \Sigma-\Delta \\
s_2 &= \Sigma+\Delta
\end{align}
$$
If we define the settling time \$T_s\$ using the same "within 2% of final response" criteria, then:
$$
\begin{align}
0.02 &= \frac{s_2}{s_2-s_1} e^{s_1 T_s} + \frac{s_1}{s_1-s_2} e^{s_2 T_s} = \\
&= \frac{\Sigma + \Delta}{2 \Delta} e^{(\Sigma - \Delta) T_s} - \frac{\Sigma - \Delta}{2 \Delta} e^{(\Sigma + \Delta) T_s} = \\
&= \frac{e^{\Sigma T_s}}{\Delta} \left[ \frac{\Sigma+\Delta}{2} e^{-\Delta T_s} - \frac{\Sigma-\Delta}{2} e^{\Delta T_s} \right] = \\
&= \frac{e^{\Sigma T_s}}{\Delta} \left[ \frac{\Delta}{2} \left( e^{\Delta T_s} + e^{-\Delta T_s} \right) - \frac{\Sigma}{2} \left( e^{\Delta T_s} - e^{-\Delta T_s} \right) \right] = \\
&= \frac{e^{\Sigma T_s}}{\Delta} \left[ \Delta \cosh{(\Delta T_s)} - \Sigma \sinh{(\Delta T_s)} \right] = \\
&= e^{K \Delta T_s} \left[ \cosh{(\Delta T_s)} - K \sinh{(\Delta T_s)} \right] = \\
&= e^{-K |\Delta| T_s} \left[ \cosh{(-|\Delta| T_s)} - K \sinh{(-|\Delta| T_s)} \right]
\end{align}
$$
And finally:
$$
0.02 = e^{-K |\Delta| T_s} \left[ \cosh{(|\Delta| T_s)} + K \sinh{(|\Delta| T_s)} \right] \\
$$
Now that we have rewritten the expression in term of \$ |\Delta| T_s\$ and \$K\$ (instead of in terms of \$s_1\$ and \$s_2\$), we can numerically solve for \$ |\Delta| T_s\$, (by simply using Excel's solver) for any arbitrary given \$\zeta>1\$.
Example 1: a moderately overdamped system with \$\zeta = 1.1\$. Thus \$K = \frac{1.1}{1.1^2-1} \approx 2.4\$, and then solving numerically:
$$
T_s \approx \frac{3.172}{|\Delta|} = \frac{3.172}{\omega_0 \sqrt{1.1^2-1}} \approx \frac{6.922}{\omega_0}
$$
Example 2: a heavily overdamped system with \$\zeta = 5\$. Thus \$K = \frac{5}{\sqrt{24}} \approx 1.0206\$, and then solving numerically:
$$
T_s \approx \frac{190.21}{|\Delta|} = \frac{190.21}{\omega_0 \sqrt{24}} \approx \frac{38.827}{\omega_0}
$$
There is also an approximation for heavily overdamped (\$\zeta \gg 1\$) systems based on the dominant pole:
$$
v_{out}(t) \approx H_0 u(t) \left[ 1 - e^{s_1 t} \right]
$$
If we define the settling time \$T_s\$ using the same "within 2% of final response" criteria, then:
$$
0.02 \approx e^{s_1 T_s}
$$
and:
$$
T_s \approx \frac{\ln(0.02)}{s_1} = \frac{-\ln(0.02)}{\omega_0 (\zeta-\sqrt{\zeta^2-1})}
$$
We can compare this approximation with the exact results that we have derived before.
For \$\zeta = 5\$:
$$
T_s \approx \frac{38.725}{\omega_0}
$$
An estimation error just about -0.25%. Quite good indeed.
For \$\zeta = 1.1\$:
$$
T_s \approx \frac{6.096}{\omega_0}
$$
An estimation error of approx -12%. Not bad taking into account that \$\zeta = 1.1\$ is just marginally above the critically damped case!.
Bonus
We can write a generic settling time expression for \$\zeta>1\$ as follows
$$
T_s = \frac{\psi}{\omega_0}
$$
where \$\psi\$ is a coefficient roughly proportional to the damping factor \$\zeta\$.
I've numerically calculated the value of \$\psi\$ for a range of \$1<\zeta<9\$ using the expression previously derived for settling within 2% of the final value,
$$
0.02 = e^{-K |\Delta| T_s} \left[ \cosh{(|\Delta| T_s)} + K \sinh{(|\Delta| T_s)} \right]
$$
Then I've calculated (for comparison purposes) 1) the dominant pole approximation, 2) a 3rd order polynomial regression on my numerically calculated dataset, and 3), 4) the relative error due to these two approximations.
Here is an Excel plot with the results:
Best Answer
For systems with real left-half plane poles, you can usually estimate it by only considering the dominant pole (the pole with the lowest frequency). In your case this would be \$p_d=-2\$. The result gets more accurate as the non-dominant pole (\$p_{nd}\$) moves further away from the dominant pole.
By only considering the dominant pole, you get a rather simple equation:
$$\begin{align} \frac{5}{4}\cdot e^{-2t}&=0.02\cdot \frac{5}{8} \\ t &= -\frac{1}{2}\cdot \ln\left( 0.02\cdot \frac{5}{8}\cdot \frac{4}{5} \right) \approx 2.30258509299s\\ \end{align}$$
The idea is that the non-dominant pole at \$p_{nd}=-4\$ leads to a term \$e^{-4t}\$ which will be damping so quickly that it doesn't affect the overall settling time. The advantage is the simplicity of the equation, and the fact that is actually a pretty common occurrence to have a very dominant pole and far-away non-dominant poles in electronic circuits.
In your specific case, it is possible to analytically calculate the settling time. The time it takes for the time-dependent terms to damp to 2% of the final value can be calculated using (similar to Andy's answer, but using the absolute value):
$$\begin{align} \left| e^{-4t}-2\cdot e^{-2t} \right| &=0.02 \\ &\Updownarrow (y=e^{-2t})\\ y^2-2\cdot y &= \pm 0.02 \\ &\Updownarrow (\text{There are 4 distinct solutions, but I only take the relevant one}) \\ y = e^{-2t} &= 1 - \frac{7}{5\sqrt{2}} \\ &\Updownarrow \\ t &= -\frac{1}{2}\cdot \ln\left(1 - \frac{7}{5\sqrt{2}}\right) \approx 2.30006613189s \end{align}$$
So a factor of 2 for \$p_{nd}/p_d\$ leads to about an error of 0.1% on the calculated settling time when using the dominant-pole approximation. Whether or not this is sufficient I leave to you.