I don't know much about superlattices, or really anything related to topological physics. I know, however, how to mindlessly apply Fourier transform to this kind of tight-binding Hamiltonian. I see, however, why you might be confused with all these different operators.
For now, just forget about the "2D" $\mathbf{k}$ sum at the beginning. We are only performing a Fourier transform along $z$. Also, you can basically forget about the Pauli matrices being operators and treat them as numbers during the Fourier transform step. Let me rewrite the Hamiltonian in a more convenient way:
$$H = \sum_{\boldsymbol{k}}\sum_{n} \Big[ v_f \tau_z(k_y \sigma_x - k_x \sigma_y) + V\tau_z + \Delta_T\tau_x \Big] c_{\boldsymbol{k},n}^{\dagger} c_{\boldsymbol{k},n}^{\phantom{.}} + \Delta_N \tau_{+} c_{\boldsymbol{k},n}^{\dagger} c_{\boldsymbol{k},n+1}^{\phantom{.}} + \Delta_N \tau_{-} c_{\boldsymbol{k},n}^{\dagger} c_{\boldsymbol{k},n-1}^{\phantom{.}} $$
(I have grouped the terms together depending on the relative value of $i$ and $j$, and replaced $i$ by $n$ so that we don't confuse it with the complex number $i$ later).
So we have three different terms: terms like $c_n^{\dagger} c_n^{\phantom{.}}$, terms like $c_n^{\dagger} c_{n+1}^{\phantom{.}}$ and terms like $c_n^{\dagger} c_{n-1}^{\phantom{.}}$, with Pauli matrices and 2D $k$ - dependent factors in front, which we will treat as numbers for now. The Fourier transform technique consists in performing the following transformations:
\begin{align}
c_{k_z}^{\dagger} &= \frac{1}{\sqrt{L}} \sum_{n} e^{+ik_znd} c_n^{\dagger}\\
c_{k_z}^{\phantom{.}} &= \frac{1}{\sqrt{L}} \sum_{n} e^{-ik_znd} c_n^{\phantom{.}},
\end{align}
with the corresponding inverse transformations:
\begin{align}
c_n^{\dagger} &= \frac{1}{\sqrt{L}} \sum_{k_z} e^{-ik_znd} c_{k_z}^{\dagger} \\
c_n^{\phantom{.}} &= \frac{1}{\sqrt{L}} \sum_{k_z} e^{+ik_znd} c_{k_z}^{\phantom{.}},
\end{align}
where $L$ is the total size of your system in the $z$ direction, and $d$ is the period of the superlattice. You can check that this choice of normalization makes the $c_{k_z}^{\dagger}$'s and $c_{k_z}^{\phantom{.}}$'s "true" fermionic operators as they verify the anticommutation relation:
$$\left\{ c_{k_z}^{\phantom{.}}, c_{k'_z}^{\dagger} \right\} = \delta_{k_z, k'_z} $$
Here because the system is finite of size $L$, the $k_z$'s can only take values which are multiples of $\frac{2 \pi}{L}$, but the same would hold for an infinite sized system ($L \to \infty$), you would just need to be more careful about the normalization.
The next step is to actually perform the substitution. Because the factors in front of each family of terms do not depend on $n$, we can simply look at the following sums and multiply them by whatever is in front of each term in the Hamiltonian:
\begin{equation}
\begin{split}
\sum_{n} c_n^{\dagger} c_n^{\phantom{.}} &= \frac{1}{L} \sum_{k_z, k'_z} \sum_{n} e^{-i(k_z-k'_z)nd} c_{k_z}^{\dagger} c_{k'_z}^{\phantom{.}} &= \sum_{k_z, k'_z} \delta_{k_z, k'_z} c_{k_z}^{\dagger} c_{k'_z}^{\phantom{.}} &= \sum_{k_z} c_{k_z}^{\dagger} c_{k_z}^{\phantom{.}} \\
\sum_{n} c_n^{\dagger} c_{n+1}^{\phantom{.}} &= \frac{1}{L} \sum_{k_z, k'_z} \sum_{n} e^{-i(k_z-k'_z)nd} e^{ik'_zd} c_{k_z}^{\dagger} c_{k'_z}^{\phantom{.}} &= \sum_{k_z, k'_z} \delta_{k_z, k'_z} e^{ik'_zd} c_{k_z}^{\dagger} c_{k'_z}^{\phantom{.}} &= \sum_{k_z} e^{ik_zd} c_{k_z}^{\dagger} c_{k_z}^{\phantom{.}} \\
\sum_{n} c_n^{\dagger} c_{n-1}^{\phantom{.}} &= \frac{1}{L} \sum_{k_z, k'_z} \sum_{n} e^{-i(k_z-k'_z)nd} e^{-ik'_zd} c_{k_z}^{\dagger} c_{k'_z}^{\phantom{.}} &= \sum_{k_z, k'_z} \delta_{k_z, k'_z} e^{-ik'_zd} c_{k_z}^{\dagger} c_{k'_z}^{\phantom{.}} &= \sum_{k_z} e^{-ik_zd} c_{k_z}^{\dagger} c_{k_z}^{\phantom{.}} \\
\end{split}
\end{equation}
All things, considered, this yields the following Hamiltonian:
$$H = \sum_{\overrightarrow{k}} \Big[ v_f \tau_z(k_y \sigma_x - k_x \sigma_y) + V\tau_z + \Delta_T\tau_x + \Delta_N e^{-ik_zd} \tau_{+} + \Delta_N e^{+ik_zd} \tau_{-} \Big] c_{\overrightarrow{k}}^{\dagger} c_{\overrightarrow{k}}^{\phantom{k}}$$
which is of the form:
$$H = \sum_{\overrightarrow{k}} H'\left(\overrightarrow{k}\right) c_{\overrightarrow{k}}^{\dagger} c_{\overrightarrow{k}}^{\phantom{k}},$$
with $\overrightarrow{k} = (\mathbf{k}, k_z) = (k_x, k_y, k_z)$ a 3D wavevector. The last step which you can try to do by yourself because it is not related to Fourier transform is to remember that $H'\left(\overrightarrow{k}\right)$ is actually an operator, which needs to be diagonalized. It can be seen as a $4 \times 4$ matrix acting on the product of the two spin spaces associated respectively with $\sigma$ and $\tau$.