Rademacher theorem

While I am still looking for references, let me sketch a proof. This will be a brief sketch only.

Theorem. If $f:\mathbb{R}^n\to\mathbb{R}^m$ is Lipschitz, differentiable at $x_o$ and $\operatorname{rank} Df(x_o)=k$, then in any neighborhood of $x_o$ the set of points satisfying $\operatorname{rank} Df\geq k$ has positive measure.

Proof. Choose $k$-dimensional subspaces $V$ and $W$ passing through $x_o$ and $f(x_o)$ respectively so that the restriction of the derivative $Df(x_o):V\to W$ is an isomorphism. $Df(x_o)$ maps a sphere $\mathbb{S}^{k-1}(\epsilon)\subset V$ of radius $\epsilon$ to the boundary $\partial\mathbb{E}(\epsilon)$ of a non-degenerate ellipsoid $\mathbb{E}(\epsilon)\subset W$. $\mathbb{S}^{k-1}(\epsilon)$ is the boundary of a ball $\mathbb{B}^k(\epsilon)\subset V$.

Let $\pi:\mathbb{R}^m\to W$ be the orthogonal projection. If $\epsilon$ is small enough, the definition of derivative implies that $(\pi\circ f)(\mathbb{S}^{k-1}(\epsilon))$ is very close to $Df(x_o)(\mathbb{S}^{k-1}(\epsilon))=\partial\mathbb{E}(\epsilon)$ and it follows that (for small $\epsilon$), $\mathbb{E}(\epsilon/2)\subset(\pi\circ f)(\mathbb{B}^k(\epsilon))$. (In fact to prove this inclusion one needs to use an argument based on degree or the fact that there is no continuous retraction of the ball onto its boundary; this is however standard and I am not going to provide details.)

Consider a torus: $\mathbb{T}^n_\delta\subset\mathbb{R}^n$ being the $\delta<\epsilon$ neighborhood of $\mathbb{S}^{k-1}(\epsilon)$. There is a natural parametrization of the torus by the product $\mathbb{S}^{k-1}(\epsilon)\times\mathbb{B}^{n-k+1}(\delta)$. For $z\in \mathbb{B}^{n-k+1}(\delta)$ let $\mathbb{S}^{k-1}_z$ be the corresponding sphere inside $\mathbb{T}^n_\delta$ and let $\mathbb{B}^k_z$ be the $k$-dimensional ball in the subspace parallel to $V$ whose boundary is $\mathbb{S}^{k-1}_z$.

Now it follows from the continuity of $f$ and $\mathbb{E}(\epsilon/2)\subset(\pi\circ f)(\mathbb{B}^k(\epsilon))$, that for $\delta$ small enough $\mathbb{E}(\epsilon/4)\subset(\pi\circ f)(\mathbb{B}^k_z)$ for all $z\in \mathbb{B}^{n-k+1}(\delta)$. Therefore the image $(\pi\circ f)(\mathbb{B}^k_z)$ has positive $k$-dimensional measure and it follows from the area formula that rank of the derivative of $\ pi\circ f$ has to be $k$ on a set of positive measure in $\mathbb{B}^k_z$. Now it follows from the Fubini theorem that the rank of derivative of $\pi\circ f$ is $k$ on a set of positive measure in $\mathbb{R}^n$ (in a small neighborhood of $x_o$). Since $\operatorname{rank}(Df)\geq \operatorname{rank}(D(\pi\circ f))$, the result follows.
$\Box$

While I don't have a reference, let me add another perspective. The following result is proved e.g. in the book Functions of Bounded Variation and Free Discontinuity Problems by Ambrosio, Fusco and Pallara (see the proof of Theorem 2.16).

Theorem. Let $n\ge k$. If $f_j\to f_\infty$ is a converging sequence of Lipschitz maps from $\mathbb R^n$ to $\mathbb R^k$ (pointwise or equivalently locally uniformly), then $$ \det(\partial_{1}f_j,\dots,\partial_kf_j)\rightharpoonup\det(\partial_{1}f_\infty,\dots,\partial_kf_\infty) $$ in the weak*-$L^\infty$ topology.

Your observation can be seen as a nice corollary of the above theorem. Indeed, assuming wlog that $m=k$ and $L(x):=\sum_{i=1}^k\partial_i f(x_0)x_i$ is an invertible linear map, you can apply the theorem to $f_j:=r_j^{-1}f(x_0+r_j\cdot)$, for any sequence of radii $r_j\downarrow 0$, getting that $$ \mathcal{L}^n(\{x\in B_1:\det(\partial_{1}f_j(x),\dots,\partial_kf_j(x))\neq 0\}>0 $$ eventually (since $f_\infty=L$), which trivially implies your claim.