Here is a look at the default behavior of ImageAdjust on a color image:
img = ExampleData[{"TestImage", "Lena"}];
adj = ImageAdjust@img;
Column[ Table[
inout =
Transpose@(Flatten[ImageData@ColorSeparate[#][[channel]], 1] & /@
{img, adj});
fit = LinearModelFit[ inout , x , x];
Row[{fit,
Show[{Plot[fit[x], {x, 0, 1}, PlotStyle -> Red],
ListPlot[inout]},
ImageSize -> 300]
}] , {channel, 3}]]

What you see is a simple linear transform on each channel. We could have got at those numbers just by looking at the data range for each channel:
Table[
range = Sort[
Flatten[ImageData@ColorSeparate[img][[channel]], 1]][[{1, -1}]];
({range[[1]] # , -#} &@ 1/(Subtract @@ range)) , {channel, 3}]
{{-0.268657, 1.26866}, {-0.0122449, 1.04082}, {-0.0368664, 1.17512}}
From this I suppose the answer to the question should be to first pass through all the images and find the absolute range over all the images for each channel, then use ImageApply to do the scaling.
Here is a look at the default behavior of ImageAdjust on a color image:
img = ExampleData[{"TestImage", "Lena"}];
adj = ImageAdjust@img;
Column[ Table[
inout =
Transpose@(Flatten[ImageData@ColorSeparate[#][[channel]], 1] & /@
{img, adj});
fit = LinearModelFit[ inout , x , x];
Row[{fit,
Show[{Plot[fit[x], {x, 0, 1}, PlotStyle -> Red],
ListPlot[inout]},
ImageSize -> 300]
}] , {channel, 3}]]

What you see is a simple linear transform on each channel. We could have got at those numbers just by looking at the data range for each channel:
Table[
range = Sort[
Flatten[ImageData@ColorSeparate[img][[channel]], 1]][[{1, -1}]];
({range[[1]] # , -#} &@ 1/(Subtract @@ range)) , {channel, 3}]
{{-0.268657, 1.26866}, {-0.0122449, 1.04082}, {-0.0368664, 1.17512}}
From this I suppose the answer to the question should be to first pass through all the images and find the absolute range over all the images for each channel, then use ImageApply to do the scaling.
Assuming all your images are the same size and channel depth you can create a per pixel, per channel scale factor base on what ImageAdjust has done and multiply each image by it.
adjustImages[keyImage_, images_] :=
With[{scaleFactor = Flatten@ImageData@ImageAdjust@keyImage/Flatten@ImageData@keyImage},
Image@ArrayReshape[scaleFactor Flatten@ImageData@#,
Dimensions@ImageData@keyImage]&/@images]
And apply so:
adjustImages[kim,{im1,im1,im3,...,imn}]
You may want to decide what the the result should be if rescaling using the key image pushes values beyond the normal range and saturates.
It may also be worth exploring what it actually is you want to adjust in the images to improve performance on whatever pipeline follows, rather than trusting ImageAdjust to optimise things to your requirements.
This is stated in the documentation (third example in Properties & Relations)
ImageAdjust[image,{c,b,γ}] is equivalent to ImageApply[Clip[(b+1)(c+1)#^γ - c/2, {0, 1}]&, image]
Consider these examples:
ref = N @ Subdivide @ 255;
img = ref~Partition~10 // Image;
ImageAdjust[img, {1, 0, 1}];
ListLinePlot[{Flatten @ ImageData @ %, ref}]
ImageAdjust[img, {0, 1, 1}];
ListLinePlot[{Flatten @ ImageData @ %, ref}]
ImageAdjust[img, {0, 0, 2}];
ListLinePlot[{Flatten @ ImageData @ %, ref}]

Major update at the bottom. First part may be obsolete.
A brute force approach:
Define a function that provides a measure of the difference between the automatically adjusted image and an image with given contrast, brightness and gamma adjustments (for now, this only works for images that are made of a raster of color triplets):
ClearAll[f];
f[c_?NumericQ, b_?NumericQ, γ_?NumericQ, im_Image] :=
Total[Map[#.# &, ImageData@ImageAdjust[im] - ImageData@ImageAdjust[im, {c, b, γ}], {2}], 2]
Get an image:
img=Import["https://i.sstatic.net/nAaCm.png"]

Now minimize the image differences:
adj = FindMinimum[f[c, b, γ, img], {{c, 0.1, 0, -1, 1}, {b, 0.1, 0, -1, 1}, {γ, 1.1, 1, 0.1, 3}}]
{14.19915417, {c -> -0.06143890985, b -> 0.09520280553, γ -> 1.094933574}}
Comparison of original image and automatically adjusted image:
Row[{img, ImageReflect[ImageAdjust[img], Left -> Right], img}]

Comparison of the semi-automatically adjusted image and fully automatically adjusted image:
Row[{ImageAdjust[img, {c, b, γ} /. adj[[2]]],
ImageReflect[ImageAdjust[img], Left -> Right],
ImageAdjust[img, {c, b, γ} /. adj[[2]]]}
]

The differences are small:
ImageSubtract[ImageAdjust[img, {c, b, γ} /. adj[[2]]], ImageAdjust[img]]

And 'small' is this:
ImageSubtract[ImageAdjust[img, {c, b, γ} /. adj[[2]]], ImageAdjust[img]] // ImageData // Max
0.0431372549
Using ImageAdjust itself to enhance those differences:
ImageSubtract[ImageAdjust[img, {c, b, γ} /. adj[[2]]], ImageAdjust[img]] // ImageAdjust

Not quite an exact match but close.
EDIT
Originally, I believed that ImageAdjust does a simple histogram stretch, scaling the maximum pixel value to 1 and the minimum to 0, basically this:
imgd = img // ImageData;
min = Min[imgd];
max = Max[imgd];
Rescale[imgd, {min, max}, {0, 1}]
But that simply didn't work out. If I use the measure introduced in the beginning to determine the difference between the {c,b,γ}-adjusted pictures and the automatically ImageAdjusted pictures on the one hand and the rescaled pictures and the automatically ImageAdjusted pictures on the other hand, the {c,b,γ} seems to be better 4 out of the 5 test images I used it on.
(
img = #;
imgd = img // ImageData;
min = Min[imgd];
max = Max[imgd];
adj = FindMinimum[f[c, b, γ, img],
{{c, 0.1, 0, -1, 1}, {b, 0.1, 0, -1, 1}, {γ, 1.1, 1, 0.1, 3}}];
{Total[Map[#.# &,
ImageData@ImageAdjust[img]-ImageData@ImageAdjust[img, {c, b, γ} /. adj[[2]]], {2}],
2],
Total[Map[#.# &, ImageData@ImageAdjust[img]-Rescale[imgd, {min, max}, {0, 1}], {2}], 2]
}
) & /@ { ...images...}
For the five pictures on the ImageAdjust doc page I get:
{{50.5394233, 81.76899059}, {14.19915417, 14.46839005}, {30.34660515, 46.0902785}, {6.459238754, 2.528418174}, {537.6580546, 1011.243348}}
It seems that ImageAdjust does more than a simple stretching of the values. It probably uses histogram equalization, which linearizes the cumulative distribution function of the pixel intensity values. This is a rather non-linear transformation. It probably requires a transformation to another color space first, to separate color values from luminance values (CIELab for instance). The transformation will then be done on the luminance component only after which an inverse color space transformation occurs.
Yet another edit
Trying to test the hypothesis Mathematica is using histogram equalization.
Load an image and turn it into grayscale (we won't do this colorspace hence and forth stuff here):
img = Import["https://i.sstatic.net/xGMge.png"]
imgG = ImageData[ColorConvert[img, "Grayscale"]];

Do histogram equalization:
ClearAll[cdf];
SetAttributes[cdf, Listable]
Scan[(cdf[#[[1]]] = #[[2]]) &,
Transpose[MapAt[Accumulate, Transpose[Sort@Tally[Flatten[imgG]]], 2]]]
imgH = (cdf[imgG] - cdf[Min[imgG]])/(Times @@ ImageDimensions[img] - cdf[Min[imgG]]);
Accumulated counts before equalization:
ListPlot@Transpose[MapAt[Accumulate, Transpose[Sort@Tally[Flatten[imgG]]], 2]]

Accumulated counts after equalization:
ListPlot@Transpose[MapAt[Accumulate, Transpose[Sort@Tally[Flatten[imgH]]], 2]]

Image comparisons:
Original grayscale image:
imgG // Image

Automatically adjusted image:
imgG // Image // ImageAdjust

Histogram equalized image:
imgH // Image

Clearly, ImageAdjust does not do a histogram equalization.
Final edit
For grayscale images at least, Mathematica seems to do a simple rescaling after all. The rescaling I used above used the minimum and maximum values measured over all color channels. In a grayscale image, rescaling with
Rescale[imgG, {Min[imgG], Max[imgG]}, {0, 1}]
yields an image that is virtually identical to the adjusted image:
(imgG // Image // ImageAdjust // ImageData) -
Rescale[imgG, {Min[imgG], Max[imgG]}, {0, 1}] // Abs // Max
2.220446049*10^-16
So, conclusion: there are no {c, b, γ} parameters to be estimated
because ImageAdjust, called without additional parameters, doesn't use
them.
With only one argument, ImageAdjust[img] performs histogram stretching, ie. pixel values, which run from min to max in img, are rescaled so they run from 0 to 1 in the ouput image.
The setting {contrast, brightness, gamma} = {0, 0, 1} does not change anything:
ImageAdjust[img, {0, 0, 1}] === img gives True.
Histogram Mean Color Balance
Your example photographs
imgList = Import /@ {"http://i.imgur.com/qzFttRD.jpg",
"http://i.imgur.com/fh9tAK1.jpg", "http://i.imgur.com/b2kWM3y.jpg",
"http://i.imgur.com/amNRIhh.jpg"};
ImageAssemble[imgList~Partition~2]

The following code balances the colors of an image by transforming the histogram of each color channel to the same common mean.
meanColorBalance[img_] := Module[{rgb, mean, dist, tDist},
rgb = ColorSeparate[img, "RGB"];
dist = HistogramDistribution[Flatten@ImageData[#]] & /@ rgb;
mean = Mean[Mean /@ dist];
tDist = TransformedDistribution[x - (Mean@# - mean), x \[Distributed] #] & /@ dist;
Inner[HistogramTransform, rgb, tDist, ColorCombine[{##}, "RGB"] &]]
mcbList = meanColorBalance /@ imgList;
ImageAssemble[mcbList ~Partition~ 2]

Applied to the "Unprocessed Color" Mars image
meanColorBalance@Import["https://i.sstatic.net/bUvXg.png"]

White Area Mean Balanced Colors
The last processed photograph can be used to determine the positions of the pixels that belong to white areas of the image.
whitePixels = PixelValuePositions[mcbList[[-1]], White, 0.2][[3000 ;;]];
The first 3000 pixel positions aren't used, as these are in the background and therefore their exact position might easily change from picture to picture.
HighlightImage[mcbList[[-1]], whitePixels, Method -> {"Compose", 0.8}]

With the following code the histogram of each color channel is shifted to get a color neutral gray for the whitePixels.
whiteAreaBalanced[img_, whitePixelsPosition_] := Module[{rgb, meanShift},
rgb = ColorSeparate[img, "RGB"];
meanShift = Mean[(# - Mean /@ #) &@PixelValue[img, whitePixelsPosition]];
Inner[Image[ImageData[#1] - #2] &, rgb, meanShift, ColorCombine[{##}, "RGB"] &]]
Finally, a comparison of the original photographs (left column) with the meanColorBalanced (middle column) and the whiteAreaBalanced photographs.
ImageAssemble[{#, meanColorBalance[#], whiteAreaBalanced[#, whitePixels]} & /@ imgList]

First approach using HistogramTransform with manual tweaking
Let's number the images i1 to i3
{i1, i2, i3} = {Import["https://i.sstatic.net/bUvXg.png"],
Import["https://i.sstatic.net/mhKLQ.png"],
Import["https://i.sstatic.net/BYku2.png"]}

Using HistogramTransform seems to be the easiest way to get close to a white balanced color image
HistogramTransform[i1, NormalDistribution[0.461, 0.207], 2]

HistogramTransform[i2, NormalDistribution[0.461, 0.207], 2]

This approach can be fine tuned
{r, g, b} = ColorSeparate[i1, "RGB"];
Manipulate[
i4= ColorCombine[{HistogramTransform[r, NormalDistribution[rm, s], 2],
HistogramTransform[g, NormalDistribution[gm, s], 2],
HistogramTransform[b, NormalDistribution[bm, s], 2]}, "RGB"],
{{rm, 0.4165}, 0.35, 0.55}, {{gm, 0.374}, 0.35, 0.55},
{{bm, 0.387}, 0.35, 0.55}, {{s, 0.213}, 0.05, 0.75}]

An additional ImageAdjust
ImageAdjust[i4, {-0.098, -0.027, 1.0975}]

gets pretty close to the "White Balanced" image.
If the upper right corner of the image is chosen to be white, the image can be further processed with
Mean@Flatten[
Map[{rc, gc, bc}*# &, ImageData@ImageCrop[i4, {32, 22}, {Left, Bottom}], {2}], 1]
Solve[% == {1, 1, 1}, {rc, gc, bc}] // Flatten
newI = Image@Map[({rc, gc, bc} /. %)*# &, ImageData[i4], {2}]

But this result is further away from the "White Balanced" reference image.
Simple white balance
As described in the Wikipedia article, the simplest white balance is to rescale the RGB channels to make white objects have white pixels. Here's a simple method where I define a white point in the original image by the 0.995 quantile of each colour channel:
image = Import["https://i.sstatic.net/bUvXg.png"];
white = Quantile[#, 0.995] & /@ Transpose[Flatten[ImageData[image], 1]]
(* {0.8, 0.737255, 0.545098} *)
ImageAssemble[{{image, ImageApply[#/white &, image]}}]

Histogram matching
In a comment to Karsten 7's answer you provided a link to some example images of the same scene taken under different conditions. To co-balance these images I propose that you isolate a region of interest and use that as a reference for HistogramTransformInterpolation.
I resized the originals to make the code faster:
files = {"http://i.imgur.com/qzFttRD.jpg", "http://i.imgur.com/fh9tAK1.jpg",
"http://i.imgur.com/b2kWM3y.jpg", "http://i.imgur.com/amNRIhh.jpg"};
images = ImageResize[Import[#], 400] & /@ files;
ImageAssemble[Partition[images, 2]]

I define a region to do the histogram matching against, using image 4 as the reference:
matchRegion = ImageTake[#, {140, 240}] &;
ref = matchRegion[images[[4]]]

Then the adjustment is like this. It attempts to match the histograms in the reference region, allowing other parts of the image (e.g. the sky) to change colour as required:
adjust[im_] := Module[{r, g, b},
{r, g, b} = HistogramTransformInterpolation[matchRegion[im], ref];
ImageApply[{r[#[[1]]], g[#[[2]]], b[#[[3]]]} &, im]]
ImageAssemble[Partition[adjust /@ images, 2]]

If I had taken a region of sky for the reference image instead, the results would all have a similar shade of blue sky but the houses would be very different.
Colour balancing
The procedure above changes the overall image brightness because it is changing the image histogram to match the reference image. Here is an attempt to change the colour balance without altering the brightness. I use a 3x3 matrix to transform the {r,g,b} triplets without changing the value of r+g+b. Technically this is not quite right, because the perceived brightness is not simply r+g+b (I believe the green component carries greater weight).
I estimate the colour correction needed by simply comparing the mean of the pixel values between two images. I've used the same sub-region as above, and this time used image 1 as the reference.
meanc[im_] := Mean[Flatten[ImageData[matchRegion[im]], 1]]
matchto[n_] := Module[{a, b, c, m},
{a, b, c} = meanc[images[[n]]]/meanc[#];
m = (1/2) {{2 a, 1 - b, 1 - c}, {1 - a, 2 b, 1 - c}, {1 - a, 1 - b, 2 c}};
ImageApply[m.# &, #]] & /@ images
ImageAssemble[Partition[matchto[1], 2]]

For comparison here is the result using image 3 as the reference:
ImageAssemble[Partition[matchto[3], 2]]












