Hint:
Since
\begin{align*}
I_z&=\int\int_S\left(x^2+y^2\right)\mathrm ds
\end{align*}
and the shape $S$ of the lamina is given by $z=x^2+y^2$ with $0\leq z\leq h$, so
\begin{align*}\mathrm ds &=\left[\left(\frac{\partial z}{\partial x}\right)^2+\left(\frac{\partial z}{\partial y}\right)^2+1\right]^{1/2}\mathrm d A\\
&=\sqrt{4x^2+4y^2+1}\;\mathrm d A
\end{align*}
It follows
\begin{align*}
I_z&=\int\int_S\left(x^2+y^2\right)\sqrt{4x^2+4y^2+1}\;\mathrm d A
\end{align*}
We can work with polar coordinates setting $x=r\cos\theta$, $y=r\sin\theta$, for $0\leq r\leq \sqrt{h}$, $0\leq \theta \leq 2\pi$ and $\mathrm d A=r\mathrm dr\mathrm d \theta$.
Then
\begin{align*}
I_z&=\int_{\theta=0}^{\theta=2\pi}\int_{r=0}^{r=\sqrt{h}}{r^2\sqrt{4r^2+1}\,r\mathrm dr\mathrm d \theta}\\
&=\int_{\theta=0}^{\theta=2\pi}\int_{r=0}^{r=\sqrt{h}}{r^3\sqrt{4r^2+1}\,\mathrm dr\mathrm d \theta}
\end{align*}