you derived the equation from the measurement values having uncertainties. For this reason, the coefficients of your functions are also uncertain. Several approaches are known to transfer the measurement uncertainties to the coefficients, for instance, propagation of uncertainty or the Monte-Carlo simulation. Please take a look to the "Guide to the Expression of Uncertainty in Measurement" to get further information, cf. https://www.bipm.org/en/committees/jc/jcgm/publications
Michael Lösler correctly provided the best known methods for propagating measurement uncertainties to predictions/estimates/interpolated. In addition, there is also the possibility of using methods of resampling repetitions.
Please consider our work below that applied resampling methods to compute the both uncertainties and bias in the context of Artificial Neural Networks (ANN). Although it has been applied to a problem involving agronomic variables and ANN, the principle is the same for interpolation using Kriging.
Article Resampling in neural networks with application to spatial analysis
there was one publication coming into my mind when reading your question:
Article A weighted total least-squares algorithm for fitting a straight line
The authors concluded that " ... The complete covariance matrix of the fitting parameters is calculated which is absolutely necessary for an evaluation of uncertainties according to the Guide to the expression of uncertainty in measurement. For cases where the straight line is not vertical, conversion equations to the usual parameters (slope and y-axis intersection) and their uncertainty matrix are provided."
I completely agree with that, and that's exactly what I'm trying to emphasize. However, there is no need to use the so-called total least squares technique because any least-squares technique is suitable to provide such as dispersion matrix. An evaluation/comparison study of the total least squares was carried out, for instance, by Frank Neitzel:
Article Generalization of total least-squares on example of unweight...
or by Malissiovas et al. in the provocative contribution:
I very much appreciate this fruitful suggestions across our own fields of expertise. The above more or less provocative contributions are a good start to be prepared for a possible twilight of the gods (of TLS). ;)
I remembered one paper:Article Two Alternative Methods for Height Transformation
in which the problem is solved with the method of Martin Vermeer. The related study material from prof. Martin Vermeer is: Vermeer, M. 2003. Heights and Height Systems in Geodesy. Postgraduate Course, Helsinki University of Technology, Helsinki.
Karin Kollo used one method for variances in section: 2. Theoretical aspects 2.1. Bilinear transformation approach.