OK. Let's first figure out a necessary condition on $a_x>0; a_{xy}, a_{xz}, a_{yt}\in(\frac 1p,1)$ for the quantity
$$
Q=\int_\Omega \frac{f_1(x)f_2(y)f_3(z)f_4(t)}{x^{a_x}(1-x^2-y^2)^{a_{xy}}|x-z|^{a_{xz}}|y-t|^{a_{yt}}}\,dx\,dy\,dz\,dt
$$
to be bounded for all $f_j\in L^{p'}$. To this end, take a small $\alpha>0$ and put
$$
f_1=f_3=\alpha^{-1/p'}\chi_{[0,\alpha]},\qquad f_2=f_4=\alpha^{-2/p'}\chi_{[1-2\alpha^2,1-\alpha^2]}\,.
$$
Then $Q$ is at least $\frac{\alpha^6\alpha^{-6/p'}}{4\alpha^{a_x+2a_{xy}+a_{xz}+2a_{yt}}}$, so we must have
$$
a_x+2(a_{xy}+a_{yt})+a_{xz}\le \frac 6p
$$
to avoid a blow-up as $\alpha\to 0$.
Now let us show that the strict inequality in the above condition is sufficient to bound $Q$. WLOG all $f_j\ge 0$. First, let's get rid of variables $z,t$. Notice that for every function $g\ge 0$ on the real line with $\int g^{p'}=1$,
$$
\int\frac{g(z)}{|x-z|^a}\approx \sum_{\gamma} \gamma^{-a}\int_{z:|x-z|\le \gamma}\le
\sum_\gamma \gamma^{-a}\min[2 Mg(x)\gamma, (2\gamma)^{1/p}]\,,
$$
where the summation is taken over all $\gamma=2^k, k\in\mathbb Z$ and $M$ is the Hardy-Littlewood maximal function. If $\frac 1p<a<1$ (which is our case), the last sum is a sum of a double-sided geometric progression whose maximum is attained at $\gamma\approx Mg(x)^{-p'}$ yielding the bound $Mg(x)^{1-p'+ap'}$. Note that this power is positive for $a>\frac 1p$ but short of $1<p'-1$ for $a<1$, $p<2$.
Using it, we reduce the estimate of $Q$ to that of the two-dimensional integral
$$
P=\iint_{0<x<\frac12<y<1\\x^2+y^2<1}\frac{f_1(x)Mf_3(x)^{1-p'+p'a_{xz}}
f_2(y)Mf_4(y)^{1-p'+p'a_{yt}}}{|x|^{a_x}(1-x^2-y^2)^{a_{xy}}}\,dx\,dy
$$
Denote the products of two functions of the same variable appearing in $P$ by $F_1(x)$ and $F_2(y)$ respectively. Then
$$
F_1\in L^{\frac{p'}{2-p'+p'a_{xz}}}=L^q,\quad F_2\in L^{\frac{p'}{2-p'+p'a_{yt}}}=L^r\,.
$$
Now do the usual trickery of partitioning the domain of integration into pieces on which the factors in the denominator are essentially constant (warning: this technique usually gives you the right interior of the admissible parameter domain but is often inconclusive on the boundary, so it is great for the back of envelope computations, but not for drawing the ultimate conclusions). More precisely, look at the region $x\approx\alpha, 1-x^2-y^2\approx\beta$ where $\alpha,\beta$ are $2^{-k}$, $k\ge 0$.
If $\beta>\alpha^2$, this region can be covered by a single rectangle $I\times J$ with $|I|\approx\alpha$, $|J|\approx\beta$. Making the trivial Holder estimate for $\int_I F_1$ and $\int_J F_2$, we conclude that the integral over such region is bounded by a constant multiple of
$$
\alpha^{-a_x}\beta^{-a_{xy}}\alpha^{1-\frac 1q}\beta^{1-\frac 1r}
=\alpha^{-a_x+\frac 2p-a_{xz}}\beta{-a_{xy}+\frac 2p-a_{yt}}\,.
$$
Since the power of $\beta$ is negative, summing over $\beta\ge \alpha^2$ is equivalent to the evaluation at $\beta=\alpha^2$, which yields
$$
\alpha^{\frac 6p-a_x-2(a_{xy}+a_{yt})-a_{xz}}\,.
$$
Our condition it that this power of $\alpha$ is positive, so summing over $\alpha$ presents no problem now.
If $\beta<\alpha^2$, we can cover our region by $N\approx\frac{\alpha^2}{\beta}$ rectangles $I_k\times J_k$ with $|I_k|\approx\frac\beta\alpha$, $|J_k|\approx \beta$ and both systems $I_k$ and $J_k$ being morally disjoint (say, having bounded covering number if you want a formal statement).
Recall that if $u_k, v_k>0$ and $q,r>1$, then
$$
\sum_{k=1}^N u_kv_k\le A(N,q,r)\left[\sum_{k=1}^n u_k^q\right]^{1/q}
\left[\sum_{k=1}^n v_k^r\right]^{1/r}
$$
where
$$
A(N,r,q)=\begin{cases}1, &\frac 1q+\frac 1r\ge 1\,;
\\
N^{1-\frac 1q-\frac 1r}, &\frac 1q+\frac 1r\le 1\,.
\end{cases}
$$
Also, if $G\in L^s$ ($s>1$), $\int G^s=1$, and $I_1,\dots,I_N$ are morally disjoint intervals of length $\ell$ each, then
$$
\sum_{k=1}^N \left[\int_{I_k} G\right]^s\le\ell^{s-1}
$$
Thus,
$$
\sum_{k=1}^N \left[\int_{I_k} F_1\right]\left[\int_{J_k} F_2\right]\le
A(N,q,r)(\beta/\alpha)^{1-\frac 1q}\beta^{1-\frac 1r}\,.
$$
If $\frac 1q+\frac 1r\ge 1$, then $A(N,q,r)=1$ and the total power on $\beta$ (taking into account $a_{xy}$ in the denominator) is
$$
\frac 4p-a_{xz}-a_{yt}-a_{xy}\,.
$$
If this power were $\le 0$, then, since $\frac 2p-a_{xz}-a_{yt}<0$, our condition would be violated, so the power is positive and, therefore, summation over $\beta<\alpha^2$ reduces to evaluating one term for $\beta=\alpha^2$, which has been considered before.
If $A(n,q,r)=N^{1-\frac 1q-\frac 1r}$, then, recalling that $N=\frac{\alpha^2}\beta$, we get the total power on $\beta$ equal to $1-a_{xy}$ with the same conclusion.
It remains to see when this sufficient condition can be satisfied for some $a_x>0; a_{xy}, a_{xz}, a_{yt}> 1-\frac 1{2p'}$. Multiplying by $p$, we see that we need the inequality
$$
6-5(p-\tfrac{p-1}2)>0
$$
to hold, which results in $\frac{p+1}2<\frac 65$ or $p<\frac 75$ (initially claimed $\frac{11}9$ was the result of a computational error: arithmetic is a much harder subject than analysis).
