If you only need to make predictions within Mathematic, then because of the low variability about the apparent smooth curve, you should use Interpolation. (It would be certainly more desirable if there was some theoretical curve function to consider.)

An alternative to interpolation, would be to use @AntonAntonov 's [Quantile regression]1.

But if you need a prediction equation to use outside of Mathematica and if two or more linear segments give what you would consider an adequate fit, then you should consider piecewise linear regression. Below is some code to estimate two linear segments along with estimating the break point.

(* Generate some data *)
data = Table[{x, 500/(x^2 + 1000)}, {x, 1, 300}];

(* Make sure that the two line segments join at the unknown cut point xc *)
sol = Solve[a1 + b1 xc == a2 + b2 xc, a2][[1]]
(* {a2 -> a1+b1 xc-b2 xc} *)

nlm = NonlinearModelFit[Log[data], 
   Boole[x <= xc] (a1 + b1 x) + Boole[x > xc] (a2 + b2 x) /. sol,
   {{a1, -0.5}, {b1, -0.13}, {b2, -2}, {xc, 3}}, x];
nlm["BestFitParameters"]
(* {a1 -> -0.36738734072224805, b1 -> -0.24323881420018847,
    b2 -> -1.8168252905594362, xc -> 3.5484803234032367} *)

Show[ListLogLogPlot[data],
 LogLogPlot[Exp[nlm[Log[x]]], {x, 1, 300}, PlotStyle -> Red]]

If you only need to make predictions within Mathematica, then because of the low variability about the apparent smooth curve, you should use Interpolation. (It would be certainly more desirable if there was some theoretical curve function to consider.)

An alternative to interpolation, would be to use @AntonAntonov 's [Quantile regression]1.

But if you need a prediction equation to use outside of Mathematica and if two or more linear segments give what you would consider an adequate fit, then you should consider piecewise linear regression. Below is some code to estimate two linear segments along with estimating the break point.

(* Generate some data *)
data = Table[{x, 500/(x^2 + 1000)}, {x, 1, 300}];

(* Make sure that the two line segments join at the unknown cut point xc *)
sol = Solve[a1 + b1 xc == a2 + b2 xc, a2][[1]]
(* {a2 -> a1+b1 xc-b2 xc} *)

nlm = NonlinearModelFit[Log[data], 
   Boole[x <= xc] (a1 + b1 x) + Boole[x > xc] (a2 + b2 x) /. sol,
   {{a1, -0.5}, {b1, -0.13}, {b2, -2}, {xc, 3}}, x];
nlm["BestFitParameters"]
(* {a1 -> -0.36738734072224805, b1 -> -0.24323881420018847,
    b2 -> -1.8168252905594362, xc -> 3.5484803234032367} *)

Show[ListLogLogPlot[data],
 LogLogPlot[Exp[nlm[Log[x]]], {x, 1, 300}, PlotStyle -> Red]]