Just to simplify the notation, I use that $u(s)$ and $f(t,s)$ vanish for $s<0$ or $t<0$, so I can remove the integration bounds and all integrals run from $-\infty$ to $\infty$. I might then as well take a Fourier transform instead of a Laplace transform, ${\cal F}(\omega)=\int e^{i\omega t}F(t)dt$. The desired relation between the transforms ${\cal F}(\omega)$ of $F(t)$ and the transforms ${\cal F}(\omega,\omega')$ of $f(s,t)$ and ${\cal U}(\omega)$ of $u(t)$ is
$${\cal F}(\omega)=(2\pi)^{-1}\int {\cal F}(\omega,\omega'){\cal U}(\omega'){\cal U}(\omega-\omega')d\omega'.$$
You started out with a double convolution and upon transformation one convolution is left.
Derivation: $${\cal F}(\omega)=\int e^{i\omega t}f(t-s,s-k)u(s)u(k)dkdsdt =$$
$$\int e^{i\omega \tau}e^{i\omega s}f(\tau,s-k)u(s)u(k)dkdsd\tau =$$
$$\int e^{i\omega \tau}e^{i\omega\sigma}e^{i\omega k}f(\tau,\sigma)u(\sigma+k)u(k)dkd\sigma d\tau =$$
$$(2\pi)^{-1}\int e^{i\omega\tau}e^{i\omega\sigma}e^{i\omega k}f(\tau,\sigma)u(\sigma+k){\cal U}(\omega')e^{-i\omega'k}d\omega'dkd\sigma d\tau =$$
$$(2\pi)^{-1}\int e^{i\omega \tau}e^{i\omega \sigma}f(\tau,\sigma){\cal U}(\omega')u(\sigma+k)e^{i(\omega-\omega')k}d\omega'dkd\sigma d\tau =$$
$$(2\pi)^{-1}\int e^{i\omega \tau}e^{i\omega \sigma}f(\tau,\sigma){\cal U}(\omega'){\cal U}(\omega-\omega')e^{-i(\omega-\omega')\sigma}d\omega'd\sigma d\tau =$$
$$(2\pi)^{-1}\int {\cal F}(\omega,\omega'){\cal U}(\omega'){\cal U}(\omega-\omega')d\omega'.$$