A practical suggestion is to divide your data into two independent datasets, then use one set to produce empirical COV function and estimate the same input functional on the points of the second dataset.
Then you can check the accuracy of your covariance function by comparison of computed results with original quantities.
You can use geostatistical techniques to create covariancie matrix and after to produce covariance function, for example it can be exponential, gaussian , spherical ,or other. Another option could be using artificial neural network
The comment by Sabah is more practical however data used for covariance function dev must be much greater in number than points to be used for varification. the data must be well distributed for both sets.
Do you need the covariance function to be isotropic? You usually do. But you will discover that your experimentally determined covariance function is anything but isotropic. What do you do then?
Yes, I need the covariance function to be isotropic where in order to perform least square collocation(LSC), we must have a model of covariance function to the success for application of LSC.
An often used approach is to compute empirical covariances. Subsequently, these values might be fitted to a pre-selected model covariance functions e.g. to the Tscherning-Rapp 1974 model