We consider the partition function $Z_{\ell}(\vec x,0\vert \vec y,t)$ of $\ell$ non-intersecting continuous directed polymers of length $t$ in dimension $1+1$, in a white noise environment, starting from positions $\vec x$ and terminating at positions $\vec y$. When $\ell=1$, it is well known that for fixed $x$, the field $\log Z_1(x,0\vert y,t)$ solves the Kardar-Parisi-Zhang equation and admits the Brownian motion as a stationary measure. In particular, as $t$ goes to infinity, $Z_1(x,0\vert y,t)/Z_1(x,0\vert 0,t) $ converges to the exponential of a Brownian motion $B(y)$. In this article, we show an analogue of this result for any $\ell$. We show that $Z_{\ell}(\vec x,0\vert \vec y,t)/Z_{\ell}(\vec x,0\vert \vec 0,t) $ converges as $t$ goes to infinity to an explicit functional $Z_{\ell}^{\rm stat}(\vec y)$ of $\ell$ independent Brownian motions. This functional $Z_{\ell}^{\rm stat}(\vec y)$ admits a simple description as the partition sum for $\ell$ non-intersecting semi-discrete polymers on $\ell$ lines. We discuss applications to the endpoints and midpoints distribution for long non-crossing polymers and derive explicit formula in the case of two polymers. To obtain these results, we show that the stationary measure of the O'Connell-Warren multilayer stochastic heat equation is given by a collection of independent Brownian motions. This in turn is shown via analogous results in a discrete setup for the so-called log-gamma polymer and exploit the connection between non-intersecting log-gamma polymers and the geometric RSK correspondence found in \cite{corwin2014tropical}.