That's really the best paper. :-( It computes just what you're looking for - or very close.
The digital object identifier (doi) for this paper is:
http://dx.doi.org/10.1119/1.11834, but it may be behind a paywall.
I don't think I can do the paper justice in a post, but I can provide a quick summary:
The total photon arrival from the source, integrated over the entire view, scales as r. So the photon arrival rate increases if you're moving towards the object.
You seem to have missed a square root in your presentation of r, I'll assume it's just a typo unless otherwise argued about.See
http://en.wikipedia.org/wiki/Relativistic_Doppler_effect for example, for the derivation of r Just In Case.The total energy delivered scales as r^2, because ##E = h \nu## and ##\nu## gets doppler shifted. The usual defintion of intensity is via delivered energy, so the intensity scales as r^2.
Relativistic aberation makes the angle subtended by the object smaller - this causes the object's apparent area to shrink by a factor of r^2. This would make the surface brightness scale as r^4 (the same energy is delivered in a smaller area). But if you are far enough away so that the object is smaller than the optical resolution of the telescope that you look at it through, this effect won't matter. YOu'll be limited by your optics, and you'll only see a r^2 increase.
Usually, stars don't show a disk , the telescope can't resolve the surface, and we talk about the "stellar magnitude" based on the total energy received, the entire visual field maps to what's effectively a point. So that's my starting assumption. With this assumptoin we get a r^2 brightness enhancent - and a doppler shift.
The authors actually work out the received brightness in the human visual range, by using a crude model of the eye's black and white frequency response.
I'll upload a few screenshot, which I think constitutes "fair use" for educational purposes, so you can at least get some information.
One post graphs their results.
Doppler factors > 1 represent motion towards the star. You can see that eventually the brightness starts to decline.
The other screen snapshot represents the equation for the unweighted spectrum S(r). You multiply this by your "response function", to filter out invisible frequencies, and that integral gives your brightness. The authors used an gaussian weighting function to represent the sensitivity of the human eye, rather than a crude square step function (which has a value of 1 if the frequency is visual and 0 if it isn't).
But I don't feel up to presenting it at the level of detail that includes their approximate visual weighting function specifically, you'll need to track down the original paper for that.
You might be wondering why the factor in front of the integral is r^-2, if the intensity scales as r^2. If you perform the integral though (the paper does it), you'll see that the end result DOES scale as r^2.
WHile you can see that the authors integrated in terms of wavelength rather than frequency, it appears that my earlier mistake was a failure to scale ##d\nu## properly under the transform. I scaled what multiplied it properly, but not ##d\nu## itself. The paper doesn't use the same approach, they use ##\lambda## instead.