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. $\phantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa}$ $\Box$