OK, let's try one more idea. Minimize the sum of $f(x,y,z)^2+A(\partial_z f(x,y,z)-1)^2$ with some positive $A$ (you'll have to play with its choice to see what works best). The advantage is two-fold: first, now you have the true root separated (the $z$ derivative is not small) and second, you have two approximate equations for it ($f=0$ and $\partial_z f=1$), which should give you extra advantage, the hope being that the parasitic roots of $f$ won't be able to match the derivative too. Also, keep the power in $z$ well below that in $x$ and $y$. Note that you are still better off than with $z-f(x,y)$ because this explicit formula satisfies all the extra conditions we are trying to impose automatically.
And yes, it'll help to see the data, though, if possible, I'd prefer the one-dimensional case (the $x$ slice of your real data should have all the same problems already). Just post something reasonable (like $S(1),\dots, S(50)$) that we want to approximate by $f(x,z)=0$ with $f$ of some reasonable degree with which you currently have a problem. I'll try to play with it a bit when I have free time.