You are correct that once we have the mean feature activations over a set of images, we normalize the network sequentially, layer by layer. There is a subtlety involved, though. You can't rescale layer weights independently of the previous layers.
Let $W_i^l$ and $b_i^l$ be the weights and bias of the $i$-th convolutional filter in layer $l$. The kernel $W_i^l$ has a 3D shape with dimensions $h \times w \times c$ (height, width, channels_in), but for easier notation down the road let's reshape it to $p \times c$, where $p = h \times w$.
$F_{ij}^l \equiv max(0,\ W_i^l \bullet P_j^{l-1} + b_i^l)$ is the activation of the $i$-th filter in layer $l$ at the $j$-th position in the activation map. Here $\bullet$ designates the convolution operation (or Frobenius inner product, or multiply-add; I adopted the symbol from Baba's answer) and $P_j^{l-1}$ is the window of $h \times w \times c = p \times c$ activations in layer $l-1$'s output that the filter convolves with at the position under consideration.
Let $$\mu_i^l \equiv \mathop{\mathbb{E}}_{X, j}F_{ij}^l = \frac{1}{NM^l} \sum_X \sum_{j=1}^{M^l} F_{ij}^l = \frac{1}{NM^l} \sum_X \sum_{j=1}^{M^l} max(0,\ W_i^l \bullet P_j^{l-1} + b_i^l)$$ be the mean activation of the $i$-th filter in layer $l$ over all $N$ images in the dataset $X$ and all $M^l$ positions in the filter's activation map. This is obviously a non-negative number, and it is actually positive for all filters in the VGG networks (when mean activations are collected over a decently sized data set).
Now suppose we "normalize" the activations by dividing weights and biases by $\mu_i^l$. This would make the mean of the activation equal to 1, if the incoming activations were the same as the original non-normalized activations. That is, $\mathbb{E}_{X, j} max(0,\ \frac{W_i^l}{\mu_i^l} \bullet P_j^{l-1} + \frac{b_i^l}{\mu_i^l}) = 1$, but only if the previous layer's activations $P_j^{l-1}$ are the same as in the original non-normalized network -- the network we calculated $\mu_i^l$ in. This is only true for the first conv layer in the normalized network, the layer that convolves with the input image. For other layers this will not only result in wrong scale, but it can actually reverse the sign of the convolution and, consequently, zero out activations after passing through the ReLU. In other words, it changes the network's output.
To fix this, we need to restore incoming activations: but we can't alter the incoming values themselves, we have to undo the normalization of the previous layer using the weights of the current layer. Note that a weight in a filter only interacts with a single channel in the previous layer. So we rescale all the weights in $W_i^l$ which interact with the $k$-th channel in layer ${l-1}$ by multiplying them by $\mu_k^{l-1}$. This cancels out the normalization of the previous layer.
To formalize, let
$D^{l-1} \equiv \begin{bmatrix} \mu_1^{l-1} & 0 & \dots & 0 \\ 0 & \mu_2^{l-1} & \dots & 0 \\\vdots & & \ddots & \\ 0 & \dots & 0 & \mu_c^{l-1} \end{bmatrix} $ be the diagonal $c \times c$ matrix constructed using all $c$ mean activations from layer $l-1$.
Then, $\mathbb{E}_{X, j} max(0,\ \frac{W_i^l{D^{l-1}}}{\mu_i^l} \bullet P_j^{l-1} + \frac{b_i^l}{\mu_i^l}) = 1$. (And this is the reason we reshaped weights to 2D, so that we can multiply matrices instead of tensors, for the sake of clarity.)
Also note that max and average pooling layers will not interfere with this scheme, because they don't alter scale.
The above probably looks more complex than it is in actual code. I pushed a GitHub repo with a short Keras implementation: https://github.com/corleypc/vgg-normalize. Looking at the sample code will likely elucidate things further.