Weird result of fitting a 2D Gauss to data

The function you are fitting can never account properly for your data: your function is a probability density function (PDF) of a bivariate Gaussian distribution – as such, the volume under it (the integral of $x$ and $y$ over the whole plane) is by definition of the PDF equal to unity. Without going into any interpolation, integration, or any fancy estimation – the volume covered by your data can be estimated as follows: the height is of order of $h=1$, and making a very crude approximation of the profile as a triangle with base $2r=100$, you get a cone with volume $\frac{1}{3}\pi r^2 h = \frac{1}{3}\pi\cdot 50^2\cdot 1 \approx 2618$ – that is not even close to unity.

Therefore, you need to fit the function $A\cdot f(x,y)$, not just $f(x,y)$. And providing reasonable starting values is generally a good habit, too (if not provided explicitly, Mathematica assumes all starting values are 1):

f[x_, y_] := Exp[-(((x - μ1)^2 + (y - μ2)^2)/(2*σ^2))]/(2*Pi*σ^2)

ListPointPlot3D[v]

fit = NonlinearModelFit[v, A f[x, y], {{A, 1000}, {μ1, 50}, {μ2, 100}, {σ, 10}}, {x, y}]

fit["ParameterTable"]

Plot3D[fit[x, y], {x, 0, 150}, {y, 0, 150}, PlotRange -> All]

Maybe you do not need to fit, just integrate. Analytically we have

Integrate[x f[x, y], {x, -∞, ∞}, {y, -∞, -∞}, Assumptions -> σ > 0]
(*μ1*)
Integrate[y f[x, y], {x, -∞, ∞}, {y, -∞, -∞}, Assumptions -> σ > 0]
(*μ2*)
Integrate[((x - μ1)^2 + (y - μ2)^2)f[x, y],{x, -∞, ∞},{y, -∞, -∞},Assumptions -> σ > 0]
(*2 σ^2*)

This is sufficient to get the parameters you need by numerical integration of your data.

Using your data we obtain

 n = Total[v[[All, 3]]];
μ1 = Total[v[[All, 1]] v[[All, 3]]/n];
μ2 = Total[v[[All, 2]] v[[All, 3]]/n]
mxy = Total[((v[[All, 1]] - mx)^2 + (v[[All, 2]] - my)^2) v[[All, 3]]/n];
σ = Sqrt[mxy/2];

{n, mx, my, σ}
(*{2284.83, 43.3037, 97.2928, 23.0803}*)

f[x_, y_] := n/(2 Pi σ^2) Exp[-((x - μ1)^2 + (y - μ2)^2)/(2 σ^2)];
GraphicsRow[{ListDensityPlot[v], DensityPlot[f[x, y], {x, 1, 128}, {y, 1, 160}]}]