Using
$$ K(ik) =
\frac{1}{\sqrt{1+k^2}}K\left(\sqrt{\frac{k^2}{k^2+1}}\right) $$
and a substitution $t^2 = \frac{k^2}{1+k^2}$, rewrite the integral as
$$ \int_0^\infty K(i k)^3\,dk = \int_0^1 K(t)^3\,dt. $$
There is a paper "Moments of elliptic integrals and critical L-values"
by Rogers, Wan and Zucker (http://arxiv.org/abs/1303.2259; also one of
the authors' earlier papers: http://arxiv.org/abs/1101.1132), and the
authors, by relating this integral to an L-series of a modular form (their theorems 1 and 2),
show that
$$ \int_0^1 K(k)^3\,dk = \frac{3}{5}K(1/\sqrt{2})^4 =
\frac{3\Gamma(\frac14)^8}{1280\pi^2}, $$
using $K(1/\sqrt{2}) = \frac14 \pi^{-1/2}\Gamma(\frac14)^2$.
$$I:=\int_{0}^{\infty}\frac{\ln{(x)}}{\sqrt{x}\,\sqrt{x+1}\,\sqrt{2x+1}}\mathrm{d}x.$$
After first multiplying and dividing the integrand by 2, substitute $x=\frac{t}{2}$:
$$I=\int_{0}^{\infty}\frac{2\ln{(x)}}{\sqrt{2x}\,\sqrt{2x+2}\,\sqrt{2x+1}}\mathrm{d}x=\int_{0}^{\infty}\frac{\ln{\left(\frac{t}{2}\right)}}{\sqrt{t}\,\sqrt{t+2}\,\sqrt{t+1}}\mathrm{d}t.$$
Next, substituting $t=\frac{1}{u}$ yields:
$$\begin{align}
I
&=-\int_{0}^{\infty}\frac{\ln{(2u)}}{\sqrt{u}\sqrt{u+1}\sqrt{2u+1}}\mathrm{d}u\\
&=-\int_{0}^{\infty}\frac{\ln{(2)}}{\sqrt{u}\sqrt{u+1}\sqrt{2u+1}}\mathrm{d}u-\int_{0}^{\infty}\frac{\ln{(u)}}{\sqrt{u}\sqrt{u+1}\sqrt{2u+1}}\mathrm{d}u\\
&=-\int_{0}^{\infty}\frac{\ln{(2)}}{\sqrt{u}\sqrt{u+1}\sqrt{2u+1}}\mathrm{d}u-I\\
\implies I&=-\frac{\ln{(2)}}{2}\int_{0}^{\infty}\frac{\mathrm{d}x}{\sqrt{x}\sqrt{x+1}\sqrt{2x+1}}.
\end{align}$$
Making the sequence of substitutions $x=\frac{u-1}{2}$, then $u=\frac{1}{t}$, and finally $t=\sqrt{w}$, puts this integral into the form of a beta function:
$$\begin{align}
\int_{0}^{\infty}\frac{\mathrm{d}x}{\sqrt{x}\sqrt{x+1}\sqrt{2x+1}}
&=\int_{1}^{\infty}\frac{\mathrm{d}u}{\sqrt{u-1}\sqrt{u+1}\sqrt{u}}\\
&=\int_{1}^{\infty}\frac{\mathrm{d}u}{\sqrt{u^2-1}\sqrt{u}}\\
&=\int_{1}^{0}\frac{t^{3/2}}{\sqrt{1-t^2}}\frac{(-1)}{t^2}\mathrm{d}t\\
&=\int_{0}^{1}\frac{\mathrm{d}t}{\sqrt{t}\,\sqrt{1-t^2}}\\
&=\frac12\int_{0}^{1}\frac{\mathrm{d}w}{w^{3/4}\,\sqrt{1-w}}\\
&=\frac12\operatorname{B}{\left(\frac14,\frac12\right)}\\
&=\frac12\frac{\Gamma{\left(\frac12\right)}\Gamma{\left(\frac14\right)}}{\Gamma{\left(\frac34\right)}}\\
&=\frac{\pi^{3/2}}{2^{1/2}\Gamma^2{\left(\frac34\right)}}
\end{align}$$
Hence,
$$I=-\frac{\ln{(2)}}{2}\frac{\pi^{3/2}}{2^{1/2}\Gamma^2{\left(\frac34\right)}}=-\frac{\pi^{3/2}\,\ln{(2)}}{2^{3/2}\,\Gamma^2{\left(\frac34\right)}}.~~~\blacksquare$$
Possible Alternative: You could also derive the answer from the complete elliptic integral of the first kind instead of from the beta function by making the substitution $t=z^2$ instead of $t=\sqrt{w}$.
$$\begin{align}
\int_{0}^{\infty}\frac{\mathrm{d}x}{\sqrt{x}\sqrt{x+1}\sqrt{2x+1}}
&=\int_{1}^{\infty}\frac{\mathrm{d}u}{\sqrt{u-1}\sqrt{u+1}\sqrt{u}}\\
&=\int_{1}^{\infty}\frac{\mathrm{d}u}{\sqrt{u^2-1}\sqrt{u}}\\
&=\int_{1}^{0}\frac{t^{3/2}}{\sqrt{1-t^2}}\frac{(-1)}{t^2}\mathrm{d}t\\
&=\int_{0}^{1}\frac{\mathrm{d}t}{\sqrt{t}\,\sqrt{1-t^2}}\\
&=2\int_{0}^{1}\frac{\mathrm{d}z}{\sqrt{1-z^4}}\\
&=2\,K{(-1)}\\
&=\frac{\Gamma^2{\left(\frac14\right)}}{2\sqrt{2\pi}}\\
&=\frac{\pi^{3/2}}{2^{1/2}\Gamma^2{\left(\frac34\right)}}.
\end{align}$$
Best Answer
$\newcommand{\+}{^{\dagger}}% \newcommand{\angles}[1]{\left\langle #1 \right\rangle}% \newcommand{\braces}[1]{\left\lbrace #1 \right\rbrace}% \newcommand{\bracks}[1]{\left\lbrack #1 \right\rbrack}% \newcommand{\ceil}[1]{\,\left\lceil #1 \right\rceil\,}% \newcommand{\dd}{{\rm d}}% \newcommand{\ds}[1]{\displaystyle{#1}}% \newcommand{\equalby}[1]{{#1 \atop {= \atop \vphantom{\huge A}}}}% \newcommand{\expo}[1]{\,{\rm e}^{#1}\,}% \newcommand{\fermi}{\,{\rm f}}% \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}% \newcommand{\half}{{1 \over 2}}% \newcommand{\ic}{{\rm i}}% \newcommand{\iff}{\Longleftrightarrow} \newcommand{\imp}{\Longrightarrow}% \newcommand{\isdiv}{\,\left.\right\vert\,}% \newcommand{\ket}[1]{\left\vert #1\right\rangle}% \newcommand{\ol}[1]{\overline{#1}}% \newcommand{\pars}[1]{\left( #1 \right)}% \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\pp}{{\cal P}}% \newcommand{\root}[2][]{\,\sqrt[#1]{\,#2\,}\,}% \newcommand{\sech}{\,{\rm sech}}% \newcommand{\sgn}{\,{\rm sgn}}% \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}} \newcommand{\ul}[1]{\underline{#1}}% \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}$
Since $\ds{{\rm erf}\pars{x} \stackrel{{\rm def.}}{=}{2 \over \root{\pi}}\int_{0}^{x}\expo{-y^{2}}\,\dd y}$: \begin{align} {\cal I}\pars{\mu}&=\int_{0}^{\infty}x^{\mu}\expo{-x^{2}} {2 \over \root{\pi}}\int_{0}^{\infty}\Theta\pars{x - y}\expo{-y^{2}}\,\dd y\,\dd x \\[3mm]&= {2 \over \root{\pi}}\int_{0}^{\pi/2}\dd\theta\,\cos^{\mu}\pars{\theta} \Theta\pars{\cos\pars{\theta} - \sin\pars{\theta}} \overbrace{\int_{0}^{\infty}\dd r\,r^{\mu + 1}\expo{-r^{2}}} ^{\Gamma\pars{1 + \mu/2}/2} \\[3mm]&={1 \over \root{\pi}}\,\Gamma\pars{1 + {\mu \over 2}} \int_{0}^{\pi/4}\dd\theta\,\cos^{\mu}\pars{\theta} \end{align}
\begin{align} I&=\lim_{\mu \to 2}\totald{{\cal I}\pars{\mu}}{\mu}= {1 \over 2\root{\pi}}\,\overbrace{\Psi\pars{2}}^{\ds{1 - \gamma}}\ \overbrace{\int_{0}^{\pi/4}\dd\theta\,\cos^{2}\pars{\theta}}^{\ds{\pars{\pi + 2}/8}} \\[3mm]&+ {1 \over \root{\pi}}\,\overbrace{\Gamma\pars{2}}^{\ds{1}} \overbrace{% \int_{0}^{\pi/4}\dd\theta\,\cos^{2}\pars{\theta}\ln\pars{\cos\pars{\theta}}} ^{\ds{\braces{4G + \pi - 2\bracks{1 + \ln\pars{2}} - \pi\ln\pars{4}}/16}} \end{align} $\Gamma\pars{z}$ and $\Psi\pars{z}$ are the $\it Gamma$ and $\it Digamma$ functions, respectively. $\gamma$ and $G$ are the $\it Euler-Mascheroni$ and Catalan constants, respectively.
$$ \begin{array}{|l|}\hline \mbox{}\\ \quad{\displaystyle\int_{0}^{\infty}x^{2}\expo{-x^{2}} \,\mathrm{erf}\pars{x}\ln\pars{x}\,\dd x}\quad \\[2mm] = \quad{{\displaystyle\quad\pars{\pi + 2}\pars{1 - \gamma} + 4G + \pi - 2\bracks{1 + \ln\pars{2}} - \pi\ln\pars{4}\quad} \over {\displaystyle 16\root{\pi}}}\quad \\ \mbox{}\\ \hline \end{array} \approx 0.0436462 $$