What I'd probably do is to take a histogram of your data to determine where the "floor," "baseline," "background value" or whatever you want to call the flat part's value. Then I'd threshold slightly above that to get a binary map of where the two peaks are. Then I'd label the binary image with bwlabel. Then you can find out which pixels belong to peak 1 or peak 2 (using find() or regionprops()). Knowing those pixels (their values and x,y locations) I think you can take the log of your data and do a simple linear least squares to get all the parameters in the equation
g(x,y) = amplitude * exp(-((x-xctr)^2+(y-yctr)^2)/sigma^2)
in other words, you can get amplitude, xctr, yctr, and sigma. After you take the log, it's linear in all the coefficients so I don't see why any non-linear stuff is needed. Here's a snippet from a demo of mine that may help you:
% Do a least squares fit of the histogram to a Gaussian.
% Assume y = A*exp(-(x-mu)^2/sigma^2)
% Take log of both sides
% log(y) = (-1/sigma^2)*x^2 + (2*mu/sigma^2)*x + (log(A)-mu^2/sigma^2)
% Which is the same as
% lny = a1*x^2 + a2*x + a3
% Now do the least squares fit.
nonZeroBins = counts > 0;
lny = log(counts(nonZeroBins));
coeffs = polyfit(bins(nonZeroBins), lny, 2)
% Find the Gaussian parameters from the least squares parameters.
sigma = sqrt(-1/coeffs(1))
mu = coeffs(2) * sigma^2/2
amplitude = exp(coeffs(3) + mu^2/sigma^2)
% Now we have all 3 Gaussian parameters,
% so reconstruct the data from the Gaussian model.
xModeled = linspace(bins(1), bins(end), 100);
yModeled = amplitude * exp(-(xModeled - mu).^2 / sigma^2);
