Hi zetafunction,
Sector decomposition gives a recursive approach to disentangling overlapping divergences.
For a more direct approach you could modify the parametric BPH-procedure in the final reference below.
Sector decomposition:
By making the change of variables x->(1-x)/x and sim for y, you can make the integrals go from 0 to 1 with the entangled divergence at x=y=0.
Then, following the first example in the Heinrich paper below and splitting the integral into a x > y = t x part and a y > x = t y part, you can disentangle the divergences.
You can do the t-integral first, using a BPH-like Taylor series subtraction at t=0 to get a finite result for the first integral. Then in the x>y sector, the second integral is divergent at x=1. Another BPH subtraction gives a final result.
Following this procedure I find the integral = 2. Of course, this is completely subtraction dependent. To find what the subtractions from the original integrand are you can reverse the procedure. (Although I have to admit, I struggled with that last bit)
I can email you my Mathematica notebook if you wish.
I've put a couple of references below. The final one is a nice explanation/proof of BPH procedure in parametric space - something that I found very useful last year.
G. Heinrich, "Sector Decomposition," International Journal of Modern Physics A, vol. 23, 2008, p. 1457.
A.V. Smirnov and V.A. Smirnov, "Hepp and Speer Sectors within Modern Strategies of Sector Decomposition," Strategies, 2008.
M.C. Bergère and J.B. Zuber, "Renormalization of Feynman amplitudes and parametric integral representation," Communications in Mathematical Physics, vol. 35, 1974, pp. 113-140.
More references at http://www.mendeley.com/research-papers/collections/15779/FeynmanIntegrals/