Denote $(x_1, \ldots, x_k)= Y$ and $(x_{k+1}, \ldots, x_n)= X$. Then the region is given by $-\epsilon \leq |X|^2-|Y|^2 \leq \epsilon$ and $|Y|^2+|X|^2 \leq 1$.
Here is an almost-solution, which unfortunately is not as straigthforwrd as I would like due to some technical issues, but the general idea should hopefully be clear and may be helpful.
We proceed in two parts. First part, is to say that the region in the problem is homeomorphic to the region $-\epsilon \leq |X|^2-|Y|^2 \leq \epsilon$ - without the other inequality - after compactifing it at infinity; this compactified region we shall call $R$.
The next step is to provide the needed homeomorphism from the new region $R$ to products of balls as desired (described below). Finally, one needs to see that the composition of the two homeomorphisms does the right thing on the relevant boundary.
For the homeomorpism from $R$ to the product of balls, the idea is to provide a diffeomorphism to $\mathbb{B}^1 \times \mathbb{B}^1$ for the case of $k=1$ $n=2$ and then use the $O(k) \times O(n-k)$ symmetry to extend it to the general case.
Let's start with the $k=1$ and $n=2$ case first. In this case $R$ is obtained by adding 4 "corner" points at infinity.
We have $ -\epsilon \leq x^2-y^2 \leq \epsilon$ which we want to map to $\mathbb{B}^1$ in two "good" ways, so that we get a map to $\mathbb{B}^1 \times \mathbb{B}^1$. We are going to give one map, and the other will be obtained by switching $x$ and $y$.
Here by "good" we mean, for one, that the preimage of $0$ for the first map is the subset $y=0$, and the boundary $x^2-y^2 = -\epsilon$ is mapped to $\pm 1$. This suggests writing $y=0$ as $-y^2=0$ and interpolating between the two by $tx^2-y^2=-\epsilon t$, so that $t= \frac{y^2}{x^2+\epsilon}$ (here one would do well to draw a picture of the family of interpolating curves between $y=0$ and $x^2-y^2=-\epsilon$). This $t$ is between $0$ and $1$ and is a good definition for $y \geq 0$, but since $t$ is invariant under $y$ going to $-y$ it does not work well for $y<0$. To fix this we define $g_1(x,y)$ to be $t$ if $y\geq 0$ and $-t$ if $y<0$. Note that this is continuous and is $O(1)=\mathbb{Z}_2$ equivariant. We define $f_1$ as the restriction of $g_1$ to the region we started with, that is one with additional condition $x^2+y^2\leq 1$.
Further, we define $f_2(x,y)=f_1(y,x)$. We claim that $F=(f_1, f_2)$ is a homeomorphism to $\mathbb{B}^1 \times \mathbb{B}^1$. Indeed, recovering $(x,y)$ from $(t,s)$ requires solving a linear system for $x^2$ and $y^2$ which has determinant $\pm ts-1$ which is non-0 outside the 4 corners where $|t|=|s|=1$, and the corners get mapped to the 4 compactification points at infinity. One checks the continuity of this inverse map separately near the 4 corners and everywhere else.
Now, for the general case.
The compactification $R$ that is needed is the addition of the "corner stratum" $\mathbb{S}^{k-1} \times \mathbb{S}^{n-1}$, that is the $|X|^2=N+\epsilon$, $|Y|^2=N$ for large $N$ (more precisely, one would define $R$ as the original region disjoint union $\mathbb{S}^{k-1} \times \mathbb{S}^{n-1}$ and put the correct topology that makes $\mathbb{S}^{k-1} \times \mathbb{S}^{n-1}$ lie at infinity).
Let $H= (\frac{Y}{|Y|} f_1(|X|, |Y|), \frac{X}{|X|} f_2(|X|, |Y|))$. Note that this is a map to $\mathbb{B}^k \times \mathbb{B}^n$, that restricted to $|X|^2- |Y|^2=-\epsilon$ it maps to $\mathbb{S}^{k-1} \times \mathbb{B}^n$ (and that if $k=1$ and $n=2$ it gives back $F$). It is a homeomorphism, for essentially the same reasons $F$ was.
Note: All of this together does give the desired homomorphism. However, if I were to rewrite it from scratch, I would probably just restrict the map $H$ to the original domain $-\epsilon \leq |X|^2-|Y|^2 \leq \epsilon$ and $|Y|^2+|X|^2 \leq 1$, then claim that its a homeomorphism to its image inside $\mathbb{B}^k \times \mathbb{B}^n$ and then claim that the image is homeomorphic to all of $\mathbb{B}^k \times \mathbb{B}^n$.