You have not provided your data, so I cannot show you how to solve the problem on your data. I guess I could make up some data. But then I would need to write the code too. :-(
The solution that I have usually found to work to create such an envelope is not that incredibly difficult, but difficulty is defined by the eyes of the beholder. Really, it is a simple enough trick:
Solve for a least squares surface, subject to the constraints that the residuals are all of the same sign. Thus here we would have
zhat - z >= 0
Apply this constraint to EVERY data point. Now use whatever model you wish to fit the surface, but you will need to employ a fitting tool that can allow inequality constraints. If you have 1000 data points, then you will have 1000 inequality constraints on the solution.
My recommendation would be to fit a surface like that I use for gridfit, so perhaps a tensor product linear surface. Were it me solving the problem, I'd add in a regularization term to smooth the surface, again like gridfit. But I would solve the problem using lsqlin, because lsqlin allows you to provide linear inequality constraints.
Simplest would be to hack gridfit, adding in the linear inequality constraints, then call lsqlin. Sadly, gridfit does not allow this as a solution. I do offer the option of such an envelope constraint in my SLM toolbox, where it works nicely enough. But SLM applies only to 1-d problems.
To show you that the idea does work very nicely, here is an example using SLM.
x = rand(1,100);
y = x.^2 - rand(size(x))/10;
slm = slmengine(x,y,'knots',0:.2:1,'plot','on','increasing','on','concaveup','on','envelope','supremum')
Since I know that the noise is always additive and negative, a supremum (least upper bpound) fit seems correct.
You need to do the same thing, but applied to a 2-d surface. Since the choice of how to do it would depend on the specific model involved, either lsqlin or fmincon would be necessary.
