Estimate Variance of Additive White Gaussian Noise (AWGN) Given Multiple Realizations with Different Mean
Basically each pixel is a realization so all you need is to work in the 3rd dimension (Though you can also get better by using the Spatial Data).

Method 1

So the trick here is to use the multiple images to estimate the Mean (True value) of each pixel and then calculate the STD on all samples (numRows * numCols * numRealizations).

Assuming we have single channel image. So we pack all give images into a tensor tI with dimensions: numRows * numCols * numRealizations.

We calculate the mean per pixel (Averaging on the 3rd dimension) and then subtract each image from the calculated average image.
Then we're left with many realizations of the noise (Well, noise and the left over from the estimated mean error).

Method 2

Estimate the STD in the 3rd dimension per pixel. Them average all numRows * numPixels estimations.

Method 3

Just as method 2. But since the property which obeys to linear operatiosn is the Variance calculate the variance along the 3rd dimension, average it over all pixels and then take the sqrt() of the average Variance.

Summary

Here is a graph showing the estimated noise STD as a function of the number of realizations:



As can be seen, for an image with many pixels (Say more than 20,000 or so) Method 3 becomes almost perfect.

The full code is available on my StackExchange Signal Processing Q61273 GitHub Repository.

Remarks


Method 1 - Basically the Maximum Likelihood Estimator for the problem. It is biased as this is the property of the ML for the Variance / STD.
Method 2 - Works nice, really intuitive.
Method 3 - Basically the MLE with the non biased estimation for the variance. As can be seen. It is also the best in real world.


Yep, Method 3 is much better only because it divides the data by N while Method 1 divides by N - 1. But it has significance impact on performance.