In fact a weaker version of Aleksandrov's theorem is true for subharmonic functions. Since $\Delta u$ is a Radon measure, the following result follows from Proposition 4.4 in [1].

>**Theorem.** If $u:\Omega\to\mathbb{R}$ is subharmonic and locally integrable, then
for almost all $x\in\Omega$ there exist a vector $Du(x)$ and a matrix $D^2u(x)$ such that
for all $1\leq p<n/(n-2)$ ($1\leq p<\infty$, if $n=2$)
$$
\left(\frac{1}{|B(x,r)|}\int_{B(x,r)}|u(y)-u(x)-u(x)(y-x)-\frac{1}{2}(y-x)^TD^2u(x)(y-x)|^p\right)^{1/p}\\=o(r^2).
$$


Although it is not stated in [AG] that the matrix $D^2u(x)$ is symmetric, you can always assume that since it appears in a quadratic form and in a quadratic form you can always replace a matrix by a symmetric one since
$$
\langle Ax,x\rangle=\big\langle\frac{1}{2}(A+A^T)x,x\big\rangle.
$$

As a corollary, the authors prove Imomkulov's theorem (see the answer of user111), but they do it independently (21 years later!) since they do not quote the work of Imomkulov.

If you want to compare it with the statement of Aleksandrov's theorem, see https://mathoverflow.net/q/352421/121665.

**[AG] G. Alberti, S. Bianchini, C. G. Stefano,** On the $L^p$-differentiability of certain classes of functions. *Rev. Mat. Iberoam.* 30 (2014), no. 1, 349–367.