Sjoerd's answer applies the power of Mathematica's very general model fitting tools. Here's a more low-tech solution.
If you want to fit a Gaussian distribution to a dataset, you can just find its mean and covariance matrix, and the Gaussian you want is the one with the same parameters. For the Gaussian I believe this is the optimal fit in a certain precise sense that I can't quite remember now -- for other distributions this may not work equally well, so there you will want NonlinearModelFit.
One wrinkle is that your data doesn't fall off to zero but to something like Min[data] $\approx 0.0387$, but we can easily get rid of that by just subtracting it off. Next, I normalize the array to sum to zero, so that I can treat it like a discrete probability distribution. (All this really does is allow me to avoid dividing by Total[data, Infinity] every time in the following.)
min = Min[data];
sum = Total[data - min, Infinity];
p = (data - min)/sum;
Now we find the mean and covariance. Mathematica's functions don't seem to work with weighted data, so we'll just roll our own.
{m, n} = Dimensions[p];
mean = Sum[{i, j} p[[i, j]], {i, 1, m}, {j, 1, n}];
cov = Sum[Outer[Times, {i, j} - mean, {i, j} - mean] p[[i, j]], {i, 1, m}, {j, 1, m}];
We can easily get the probability distribution for the Gaussian with this mean and covariance. Of course, if we want to match the original data, we have to rescale and shift it back.
f[i_, j_] := PDF[MultinormalDistribution[mean, cov], {i, j}] // Evaluate;
g[i_, j_] := f[i, j] sum + min;
The match is not too bad, although the data has two peaks where the Gaussian has one. You can also compare the fit through a contour plot. (Use the tooltips to compare contour levels.)
Show[ListPlot3D[data, PlotStyle -> None],
Plot3D[g[y, x], {x, 1, 9}, {y, 1, 9}, MeshStyle -> None,
PlotStyle -> Opacity[0.8]], PlotRange -> All]
Show[ListContourPlot[data, Contours -> Range[0.05, 0.25, 0.025],
ContourShading -> None, ContourStyle -> ColorData[1, 1], InterpolationOrder -> 3],
ContourPlot[g[y, x], {x, 1, 9}, {y, 1, 9}, Contours -> Range[0.05, 0.25, 0.025],
ContourShading -> None, ContourStyle -> ColorData[1, 2]]]
The variances along the principal axes are the eigenvalues of the covariance matrix, and the standard deviations (which I guess you're calling the semi-axes) are their square roots.
Sqrt@Eigenvalues[cov]
(* {1.86325, 1.50567} *)