sorry for necromancing this thread but I wanted to give an elementary (with physics style abuse of notation) answer to this:
Yes there are some simple generalizations that support this. We need to recall where the Euler-Maclaurin formula comes from:
Let $I$ denote the identity operator. Let $H$ denote the shift operator $H[f] = f(x+1)$. Then (with abuse of notation) we can compute the following:
$$ \frac{I}{I-H} = I + H + H^2 + H^3 ... = f(x)+f(x+1)+f(x+2)+f(x+3)...$$
So then it follows that
$$ - \frac{I}{I-H} |_{x=a}^{x=b} = f(a)+f(a+1)+ \ ... \ + f(b-1) $$
Now if you recall that $H = e^{\frac{d}{dx}} $
We then have
$$ \sum_{k=a}^{b-1} f(x) = - \frac{I}{I-e^{\frac{d}{dx}}} |_{x=a}^{x=b} $$
Now recall on the right hand side that
$$ - \frac{1}{1-e^x} = \frac{1}{x} - \frac{1}{2} + \frac{1}{12}x + ... $$
So we conclude then that:
$$ - \frac{I}{I-e^{\frac{d}{dx}}} |_{x=a}^{x=b} = \int_{a}^{b} f dx - \frac{1}{2}(f(b)-f(a)) + \frac{1}{12} (f'(b)-f(a)) + ... $$
Which is the Euler Maclaurin formula.
Now we go to the multivariable case. Suppose we have a function $f(x,y)$ we have shift operators $H_x = H(x+1,y)$ and $H_y = H(x,y+1)$
It's easy to see that $H_x, H_y$ commute so that with abuse of notation we can write
$$ \sum_{n=a}^{b-1} \left[ \sum_{k=c}^{d-1} f(n,k) \right] = \frac{I}{I-H_x } \frac{I}{I-H_y}|_{x=a}^{x=b} |_{y=c}^{y=d} = \frac{I}{I-e^{\frac{\partial}{\partial x}}} \frac{I}{I-e^{\frac{\partial}{\partial y}}}|_{x=a}^{x=b} |_{y=c}^{y=d} $$
Thereby giving you that:
$$ \sum_{n=a}^{b-1} \left[ \sum_{k=c}^{d-1} f(n,k) \right] = \int_{c}^{d} \left( \int_{a}^{b} f(x,y) \partial x - \frac{f(b,y)-f(a,y)}{2} + \frac{ \partial_x[f]_{x=b} - \partial_x[f]_{x=a} }{12} + ... \right) \partial y + \\ -\frac{1}{2} \left( \int_{a}^{b} f(x,y) \partial x - \frac{f(b,y)-f(a,y)}{2} + \frac{ \partial_x[f]_{x=b} - \partial_x[f]_{x=a} }{12} + ... \right)|_{y=c}^{y=d} + ... $$
Which is just the Euler Maclaurin Formula for variables $x$ and $y$ composed with each other.
In general moving to $\mathbb{Z}^d$ follows naturally from this by just consider each variable $x_1 ... x_d$ separately and composing them just like we did for $\mathbb{Z}^2$ above.
As a more general corollary. If we want to evaluate
$$ f(x+c_0) + f(x+c_1) + f(x+c_2) + ... $$
Then we can let $T(x) = \sum_{n=0}^{\infty} x^{c_n}$. Then the operator $T(e^{\frac{d}{x}})$ when expanded via the Laurent series for $T(e^x)$ will give you the corresponding euler maclaurin formula.
As an example if we want to evaluate
$$ f(x) + f(x+1) + f(x+\sqrt{2}) + f(x+\sqrt{3}) + ... $$
We look at $\sum_{n=0}^{\infty} x^{\sqrt{n}} = \frac{\Gamma(1+2)}{(-z)^2} + \sum_{k=0}^\infty \zeta\left(-\frac{k}{2}\right)\frac{z^k}{k!}$ as per here: and therefore we conclude
$$ f(x+\sqrt{a}) + f(x+\sqrt{a+1}) + f(x+\sqrt{a+2}) + ... f(x+\sqrt{b-1}) = \\ \underbrace{-2\int \int f(x) dx}_{\text{bounds selected carefully}} - \frac{f(b)-f(a)}{0!2} - \frac{\zeta(-\frac{1}{2})(f'(b)-f'(a))}{1!} + \frac{(f''(b)-f''(a))}{2!12} + ... $$