Here's an idea that works if you have 1D interpolation routines and the query points are on a regular grid.
You can compute the SVD of the data matrix, $Y = U\Sigma V^T$, interpolate the columns of $U$ and $V$, to obtain, $\hat{U}$ and $\hat{V}$, and then compute the interpolated values as $\hat{U} \Sigma \hat{V}^T$. If $Y$ is approximately low-rank, you can discard singular vectors corresponding to small singular values.
Note that when interpolating the columns of $U$ and $V$, determining which pair of nodes is closest to the query point only needs to be done once for all columns.
If the number of query points is large, forming the output via this matrix-matrix product can be faster than direct bilinear interpolation.
Here's some matlab code that demonstrates this idea:
x = linspace(-1, 1, 100);
[X,Y] = meshgrid(x);
A = exp(-(X.^2+Y.^2)).*sin(2*pi*X.*Y);
xq = linspace(-1, 1, 1000);
[Xq,Yq] = meshgrid(xq);
tic
Aq = interp2(X, Y, A, Xq, Yq);
toc
tic
[U,S,V] = svd(A);
sv = diag(S);
Uq = interp1(x, U, xq);
Vq = interp1(x, V, xq);
AqSVD = Uq*(sv.*Vq');
toc
norm(AqSVD - Aq, inf)/norm(Aq, inf)
Elapsed time is 0.033942 seconds.
Elapsed time is 0.012933 seconds.
ans =
1.29298230462322e-15