How to calculate the double integral of RBFInterpolater?

Question

I need to use multiquadric to calculate the volume of the point cloud. When using the Rbf class, the dblquad function can be used to calculate the double integral. But with RBFInterpolator class, I don't know how to achieve similar functionality and the bugs appeared (picture below) : Here is the Sample code:

from scipy import interpolate
import numpy as np
from scipy.integrate import dblquad

if __name__ == '__main__':
    points = np.array([[1, 4, 5], [2, 5, 7], [3, 6, 9]])
    # ** 1.Rbf ** #
    # 3d points, eg.[x,y,z]
    func1 = interpolate.Rbf(points[:, 0], points[:, 1], points[:, 2], smooth=0, function='multiquadric')
    volume1, _ = dblquad(func1, 1, 2, 3, 4)
    print(volume1)

    # ** 2.RbfInterpolator ** #
    func2 = interpolate.RbfInterpolator(points[:, 0:2], points[:, 2], smoothing=0, kernel='multiquadric', epsilon=1)
    volume2, _ = dblquad(func2, 1, 2, 3, 4)
    # print(volume2)

The picture below is a point cloud example. Grey points are original points sensor collected. Blue points are the points I entered in Rbf(multiquadric), and green points are the output results. But if you watch closely, you will find the green curved surface does not go through part of the blue points. The double integral in example code represents space volume above the green curved surface. My initial idea is using RbfInterpolator with smoothing=0 to see whether they can across. I am wondering whether there is a better way can solve this problem. Looking forward to your reply, really appreciate it!