A boundary element formulation problem

I am working on some boundary integral equation formulation and I am currently stuck with some mathematics. I wish anyone can help me out with this. 

I have an anisotropic (sometimes called generalized) biharmonic differential operator which  takes the form

L = k11 D1^4 +  k12 D1^2 D2^2 + k22 D2^4 

where D1 = d/dx, D2 = d/dy, my problem is two dimensional.

I need to find a fundamental solution (Green's function) for this operator, that is 

L(u) = -delta

where delta is Dirac delta function.

I thought of two solutions but they both seem to fail, first, i thought of making some coordinate transformation so that the coefficients of the terms of the operator become the same, and hence, use the traditional biharmonic Green's function (1/8/pi r^2 ln(r) ), however, my proposed transformation takes the form 

xi = a11 x + a12 y

eta = a21 x + a22 y

there a11, a12, a21 and a22 are to be determined according to the above requirement (all operator terms have equal coefficients), however, using this transformation got me more terms in the final operator form (terms involving D1^3 D2 and D1 D2^3)

My second solution was to decompose the anisotropic biharmonic operator into two anisotropic laplacian operators and solving two anisotropic laplacian equations instead of one, I succeeded in the decomposition and I've solved the first anisotropic laplacian equation, however, I cannot solve the seconf one. the solution I got was the following

L = (a1 D1^2 + a2 D2^2)*(b1 D1^2 + b2 D2^2)

let L1 =  (a1 D1^2 + a2 D2^2) 

and

L2 =  (b1 D1^2 + b2 D2^2)

now, L1(L2(u)) = -delta

let L2(u) = v

then L1(v) = -delta

this gives v = -1/sqrt(a1*a2)*ln(r')

where r' = sqrt(x^2/a1 + y^2/a2)

now, we have L2(u) = v

thus L2(u) =  -1/2/sqrt(a1*a2)*ln(x^2/a1 + y^2/a2)

which gives  (b1 D1^2 + b2 D2^2)(u) = -1/2/sqrt(a1*a2)*ln(x^2/a1 + y^2/a2)

I am stuck at this point, I have no clue how to solve this

any ideas ??

thanks

Ahmed