There is a simple way to test whether your scaling is correct: the amplitude of the $0$ frequency y_m(0,0)
should be proportional to the mean of the windowed data (the one you apply the FFT too). Something similar was discussed in What are the units of my data after an FFT?, and in DCT and mean difference of an image for the constant of proportionality, related to the number of samples, or its square root.
Hence, the amplitude of the PSD (at $0$-frequency) should be proportional to the square of the mean of the windowed data. This recipe should help you in general for normalizations with linear transformations: track the average and the $0$ or DC frequency and compensate for the scalar term (and the power).
However, I am troubled by your 1D normalization. Two procedures are detailed in Spectrum and spectral density estimation by the Discrete Fourier transform (DFT), including a comprehensive list of window functions and some new flat-top windows (Section 9: Scaling the results).
The first one, called 'power spectrum' (PS) is with your $S_1$. I am not too familiar with it [and welcoming comments by other folks on SE.DSP], so you could go for the square of the sum ($S$) of the 2D window coefficients:
$$S^2 = \left(\sum_{i=0}^m \sum_{j=0}^n w_i w_j\right)^2\,.$$
In my experience, in most of the cases in the frequency domain, one normalizes by the energy of the window (called $S_2$ and PSD in the above paper), or its square root. Here you normalize by the square of the sum of the window coefficients. Quite non-standard for a PSD to me! And this would not make sense with a $0$-mean window (like a Gaussian derivative). In that case, I would normalize by the energy of the product:
$$\mathcal{S} = \sum_{i=0}^m \sum_{j=0}^n \left(w_i w_j \right)^2\,.$$
So you proposal is unlikely: you can see an inhomogeneity in the units: square on top, fourth-power on bottom.
The most important question remains: why do you want to normalize (several possible answers)?
And a final comment: one should call it an Hann window (not Hanning), and I am not sure the tensor-product of two Hann windows does qualify for a "2D Hann" (circular versions probably more appropriate).