In this paper, we consider the equation $$plus _{B} ^{k}u(x)=um_{r=o}^{m}c_{r}plus _{B}^r ẹlta,$$ where $\oplus _{B} ^{k}$ is the operator iterated $k$-time and is defined by $$plus _{B} ^{k}=eft[eft(B_{x_{1}}+B_{x_{2}}+\cdots+B_{x_{p}}\right)^{4}-eft(B_{x_{p+1}}+B_{x_{p+2}}+\cdots+B_{x_{p+q}}\right)^{4}\right]^{k},$$ where $p+q=n, x=(x_{1},\ldots , x_{n})\in \mathbb{R}^{+}_n$, $B_{x_{i}}=\frac{\partial ^{2}}{\partial x_{i}^{2}}+ \frac{2v_{i}}{x_{i}}\frac{\partial }{\partial x_{i}}$, $v_{i}=2\alpha _{i}+1$, $\alpha _{i}>-\frac{1}{2}$, $x_{i}>0$, $i=1,2,\ldots,n$, $c_{r}$ is a constant, $k$ is a nonnegative integer, $\delta$ is the Dirac-delta distribution, $\oplus _{B} ^{0}\delta =\delta$ and $n$ is the dimension of $\mathbb{R}^{+}_n$. It is shown that, depending on the relationship between $k$ and $m$, the solution to this equation can be ordinary functions, tempered distributions, or singular distributions.