Using the median isn’t a good way to exclude outliers, especially if a significant number of pixels are outliers (stars in a busy field). If, let’s say, 10% of pixels contain stars, then those outliers at the top of your distribution effectively mean that your median gives you the 55% percentile of the sky background distribution, which might be rather a noticeable skew towards the brighter end of the distribution.
Of course, that would be far better than taking the mean of a distribution containing stars, because the latter really would give you a nonsensical result. But to do this well, you really want to entirely get rid of the stars from your sample of pixels.
As Nick says, the way to do that is via “sigma clipping”. Take the mean and standard deviation of your pixel brightnesses. Exclude all pixels which are more than X standard deviations away from the mean. Iterate a few times.
I’ve come across this problem in the context of doing background subtraction from the images from my meteor cameras. In that context, the median has an additional issue that it’s quantised. If you take the median of a large sample of integer values between 0 and 255, you can’t get the noise below 1 ADU unit. With the mean of large numbers of samples, my sky background model sometimes reaches a formal precision of about ~ 0.2 ADU unit.
I would assume Nick has enough bit depth in his data that quantisation is unlikely to be a worry, though.