Why Zeta regularization is not valid for multiple-loops? Why zeta regularization only valid at one-loop? 
I mean there are zeta regularizations for multiple zeta sums.
Also we could use the zeta regularization iteratively on each variable 
to obtain finite corrections to multiple loop diagrams.
 A: The standard use of zeta-regularization in quantum field theory is to replace the ill defined one-loop functional determinant with something finite:
$$\begin{align} 
   \det(H) &= \exp(\mathrm{tr}\log H) := \exp(-\zeta'_H(0))\,, \\
   \zeta_H(s) &= \mathrm{tr} H^{-s} = \frac1{\Gamma(s)}\int_0^\infty\mathrm{d}t\, t^{s-1}\mathrm{tr}(\exp(-H t))
\end{align}$$
This is a generalization of the 
Riemann zeta-function $\zeta(s)=\zeta(s,q)$ and 
Hurwitz zeta-function $\zeta(s,q)=\sum_{n=1}^\infty (q+n)^{-s}$ 
since $$\zeta_H(s) = \mathrm{tr} H^{-s} = \sum_{n} \lambda_n^{-s}\,.$$ For example, for the harmonic oscillator, $\lambda_n=n+\frac12$ so $\zeta_H(s)=\zeta(s,\frac12)$.
This obviously only works at one-loop because the quantum corrections only look like a functional determinant at one-loop.
But related ideas can be extended to multiple loops. (How does the following compare to your scheme? It's hard to tell exactly what you're saying in your paper...)
In your previous question I mentioned operator regularization (e.g. PhysRevD.35.3854) as an extension of zeta-regularization.
Here, the functional determinant is written using 
$$\begin{align} 
   \det(H) &= \exp(\mathrm{tr}\log H) \\
 \log H &:= \lim_{s\to0}\Big[-\frac{d^m}{d s^m}\Big[\frac{s^{m-1}}{m!}H^{-s}\Big]\Big]\,,
\end{align}$$
for some large enough (only actually need $m=1$) integer $m$. With $m=1$ and writing $H^{-s}=\frac1{\Gamma(s)}\int_0^\infty \mathrm{d}t\, t^{s-1}\exp(-H t)$ you recover the zeta-regularized form given above.
Similarly you can also regularize the propagators $H^{-1}$ as
$$    H^{-1} = \frac{d}{d H}\log H := \lim_{s\to0}\Big[\frac{d^m}{d s^m}\Big[\frac{s^{m}}{m!}H^{-s-1}\Big]\Big]\,.
$$
This is the underlying basics of operator regularization. 
In arXiv:1006.1806, this was generalized to 
$$    H^{-n} := \lim_{s\to0}\Big[\frac{d^m}{d s^m}\Big[
      \big(1+\alpha_1s+\dots+\alpha_ms^m\big)\frac{s^{m}}{m!}H^{-s-n}\Big]\Big]\,.
$$
where the $\alpha_i$ are initially arbitrary finite constants and $m$ is the loop order (or higher).
Comparisons were made with some dimensionally regularized Feynman integrals and it was found that for a choice of $\alpha_i$ they matched. 
These constants are fixed according the renormalization scheme that is used. I'm not sure if  a "minimal subtraction" scheme exists for operator regularization (unless you do a comparison with dimensional regularization for each diagram!) or if a physical/on-shell renormalization scheme has to be used to ensure consistent results.
A much more solid paper is Operator regularization and multiloop Green's functions.
There they define the subtraction operator (with out the $\alpha_i$ constants)
$$ \mathcal S^{(m)}\big(\prod_i A_i^{-1}\big) 
   = \lim_{s\to0}\frac{d^m}{d s^m}\Big[\frac{s^{m}}{m!}\Big(\prod_i A_i^{-s-1}\Big)\Big]
$$ 
and show that its naive application to obtain finite amplitudes breaks unitarity.
(The results of arXiv:1006.1806 suggest that if you include the $\alpha_i$'s and choose them correctly, you should be able to get correct results...)
Then they actually use Bogolibov's recursion formula to show how to construct a consistent $\mathcal R$-operator from $\mathcal S$.
Detail computations are given for both the naive and sophisticated use of operator regularization in massive $\phi^4_4$ theory.
Note that using a BPHZ-like scheme essentially reintroduces counterterms (hidden in the subtractions in each diagram) to the method and definitely breaks the simplicity of the original $\zeta$-regularization.
