Skyrmion number

It is not a curl. This can be seen by expressing the curl in vector components. $$nabla times mathbf M=begin{pmatrix} partial_yM_z-partial_z M_y partial_zM_x-partial_x M_z partial_xM_y-partial_y M_x end{pmatrix}$$ Here $partial_x$ denotes the partial derivative with respect to $x$. The quantity $partial_xmathbf M$ is a vector just like $mathbf M$. It has components $$partial_x mathbf M=begin{pmatrix} partial_xM_x partial_xM_y partial_xM_z end{pmatrix}$$ Calculating the quantity $partial_xmathbf Mtimespartial_ymathbf M$ is then just a matter of applying the cross product. $$partial_xmathbf Mtimespartial_ymathbf M=begin{pmatrix} partial_xM_ypartial_yM_z-partial_xM_zpartial_yM_y partial_xM_zpartial_yM_x-partial_xM_xpartial_yM_z partial_xM_xpartial_yM_y-partial_xM_ypartial_yM_x end{pmatrix}$$ This is a daunting expression and you probably won't get a lot of intuition from looking at the components. What you can say about it is that $mathbf Acdot(mathbf Btimes mathbf C)$ forms the vector triple product. This gives the volume spanned by (the parallelepiped of) $mathbf A,mathbf B$ and $mathbf C$. So the quantity you're integrating is the volume spanned by $mathbf M,partial_x mathbf M$ and $partial_y mathbf M$.

To calculate the integral in your last equation is just a matter of plugging everything in in my last expression for $partial_xmathbf Mtimespartial_ymathbf M$. This is tedious but should be doable.

---

EDIT I'll add some more info to make the calculation less tedious. The partial derivatives can be expanded using the chain rule $partial_x=frac{partial r}{partial x}partial_r+frac{partial alpha}{partial x}partial_alpha$. These can be calculated to be $$partial_x=cosalphapartial_r-frac{sinalpha}rpartial_alpha partial_y=sinalphapartial_r+frac{cosalpha}rpartial_alpha$$ Next note that $partial_rmathbf M=frac{dtheta}{dr}partial_thetamathbf M$ and $partial_alphamathbf M=frac{dphi}{dalpha}partial_phimathbf M$. If we name these partial derivative vectors $mathbf e_theta=partial_thetamathbf M$ and $mathbf e_phi=partial_phimathbf M$ then the cross product becomes $$partial_xmathbf Mtimespartial_ymathbf M=left(cosalphafrac{dtheta}{dr}mathbf e_theta-frac{sinalpha}rfrac{dphi}{dalpha}mathbf e_phiright)timesleft(sinalphafrac{dtheta}{dr}mathbf e_theta + frac{cosalpha}rfrac{dphi}{dalpha}mathbf e_phiright)$$ Finally you can calculate that $mathbf e_thetatimes mathbf e_phi=sintheta ,mathbf M$ and you should be able to do this calculation without explicitly calculating all the components.

And yes you should add the factor $r$ when you switch to polar coordinates like you mentioned in your comment.