There are multiple (equivalent) ways of deriving the canonical form. One of them—from the original Gorini et al. paper on the Lindblad form—is to start from some orthonormal basis $\{F_j\}_j$ of $\mathbb C^{n\times n}$ (with respect to the usual Hilbert-Schmidt inner product) and to rewrite the Lindbladian as
$$
\mathcal L=-i[H,\cdot]+\sum_{j,k}c_{jk}\Big(F_j(\cdot)F_k^\dagger -\frac12(\cdot)F_k^\dagger F_j-\frac12F_k^\dagger F_j(\cdot) \Big)\,.
$$
The resulting matrix $C:=(c_{jk})_{j,k}\in\mathbb C^{n^2\times n^2}$—which is also called Kossakowski matrix—is positive semi-definite because $\mathcal L$ is a Lindbladian so $C$ can be diagonalized as $C=\sum_j\mu_j|g_j\rangle\langle g_j|$ for some $\mu_j>0$. If I recall correctly, the orthonormal system $\{g_j\}_j\subset\mathbb C^{n^2}$ tells us how to build a new orthonormal basis $\{L_j\}_j$ from the old $\{F_j\}_j$ one (i.e. via linear combinations according to the entries of $g_j$) which thus form the (orthonormal) jump operators $L_j$.
The second way to obtain the canonical form—which is probably more instructive as it is slightly more explicit—is via conditional complete positivity (cf. also this answer). The idea here is to diagonalize (some variant of) the Choi matrix, much like how one finds Kraus operators for a given completely positive map. Now starting from the Lindbladian (or, more generally, a linear, Hermitian-preserving, and conditionally completely positive map) $\mathcal L$ we consider the matrix
$$
W:=n({\bf1}-|\Omega\rangle\langle\Omega|)\big( ({\rm id}\otimes\mathcal L)(|\Omega\rangle\langle\Omega|) \big)({\bf1}-|\Omega\rangle\langle\Omega|)
$$
where $|\Omega\rangle:=\frac1{\sqrt n}\sum_j|j\rangle\otimes|j\rangle$ is the maximally entangled state. Because $\mathcal L$ is of Lindblad form (more precisely: because $\mathcal L$ is conditionally completely positive) the matrix $W$ is positive semi-definite meaning we can again diagonalize it: $W=\sum_j\mu_j|g_j\rangle\langle g_j|$. Then the constructive proof that conditional complete positivity implies the form $K(\cdot)+(\cdot)K^\dagger +\Phi$ shows that the orthonormal operators $L_j:={\rm vec}^{-1}(g_j)\in\mathbb C^{n\times n}$ are the Lindblad/jump operators we were looking for; here ${\rm vec}$ is the usual vectorization and ${\rm vec}^{-1}$ is the inverse of that, i.e. it re-arranges a big vector back into a square matrix. Altogether, these $L_j$ and the eigenvalues $\mu_j$ of $W$ satisfy
$$
\mathcal L=-i[H,\cdot]+\sum_j\mu_j\Big(L_j(\cdot)L_j^\dagger -\frac12(\cdot)L_j^\dagger L_j-\frac12 L_j^\dagger L_j(\cdot)\Big)\,,\tag3
$$
where $H:=\frac i2(K-K^\dagger )$ is the corresponding Hamiltonian.
Let me illustrate this procedure by applying it to your example.
First a disclaimed: I'll drop the time-dependence because it does not change anything and just makes the expressions more lengthy.
Now we start from
\begin{align*}
\dot{\rho} &= \mathcal L(\rho)=(\tfrac32\gamma + \tilde{\gamma}) (2 \sigma_x \rho \sigma_x + 2 \sigma_y \rho \sigma_y - 4\rho) \nonumber \\&- \gamma (2 \sigma_- \rho \sigma_+ - \sigma_+ \sigma_- \rho - \rho \sigma_+ \sigma_-) \nonumber \\&- \gamma (2 \sigma_+ \rho \sigma_- - \sigma_- \sigma_+ \rho - \rho \sigma_- \sigma_+)
\end{align*}
(refer to my comment for where the factor $\frac32$ is coming from).
Next we compute the Choi matrix
\begin{align*}
2({\rm id}\otimes\mathcal L)(|\Omega\rangle\langle\Omega|)=\begin{pmatrix}
-4 (\gamma +\tilde\gamma ) & 0 & 0 & -4 (\gamma +\tilde\gamma ) \\
0 & 4 (\gamma +\tilde\gamma ) & 0 & 0 \\
0 & 0 & 4 (\gamma +\tilde\gamma ) & 0 \\
-4 (\gamma +\tilde\gamma ) & 0 & 0 & -4 (\gamma +\tilde\gamma )
\end{pmatrix}
\end{align*}
as well as the $W$-matrix:
\begin{align*}
W&=2({\bf1}-|\Omega\rangle\langle\Omega|)\big( ({\rm id}\otimes\mathcal L)(|\Omega\rangle\langle\Omega|) \big)({\bf1}-|\Omega\rangle\langle\Omega|)\\
&=
\begin{pmatrix}
\frac{1}{2} & 0 & 0 & -\frac{1}{2} \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
-\frac{1}{2} & 0 & 0 & \frac{1}{2}
\end{pmatrix}
\begin{pmatrix}
-4 (\gamma +\tilde\gamma ) & 0 & 0 & -4 (\gamma +\tilde\gamma ) \\
0 & 4 (\gamma +\tilde\gamma ) & 0 & 0 \\
0 & 0 & 4 (\gamma +\tilde\gamma ) & 0 \\
-4 (\gamma +\tilde\gamma ) & 0 & 0 & -4 (\gamma +\tilde\gamma )
\end{pmatrix}
\begin{pmatrix}
\frac{1}{2} & 0 & 0 & -\frac{1}{2} \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
-\frac{1}{2} & 0 & 0 & \frac{1}{2}
\end{pmatrix}\\
&=\begin{pmatrix}
0 & 0 & 0 & 0 \\
0 & 4 (\gamma +\tilde\gamma ) & 0 & 0 \\
0 & 0 & 4 (\gamma +\tilde\gamma ) & 0 \\
0 & 0 & 0 & 0
\end{pmatrix}\,.
\end{align*}
Luckily $W$ is already diagonal: $W=4 (\gamma +\tilde\gamma )|01\rangle\langle 01|+4 (\gamma +\tilde\gamma )|10\rangle\langle 10|$. Therefore $\mathcal L$ satisfies Eq.(3) with jump operators
\begin{align*}
L_1&={\rm vec}^{-1}|01\rangle={\rm vec}^{-1}\begin{pmatrix}0\\1\\0\\0\end{pmatrix}=\begin{pmatrix}0&0\\1&0\end{pmatrix}=\sigma_-\\
L_2&={\rm vec}^{-1}|10\rangle={\rm vec}^{-1}\begin{pmatrix}0\\0\\1\\0\end{pmatrix}=\begin{pmatrix}0&1\\0&0\end{pmatrix}=\sigma_+
\end{align*}
and rates $\gamma_1=\gamma_2=4 (\gamma +\tilde\gamma )$ (also because there was no Hamiltonian part to begin with, so $H=0$).
The discrepancy to Eq.(2) comes from the fact that the Lindblad operators are only unique up to unitary transformations (cf. the end of this section or
Eq.(3.72) in the book of Breuer & Petruccione)
which in this case is
\begin{align*}
\frac1{\sqrt2}\sigma_x&=\frac1{\sqrt2}\sigma_++\frac1{\sqrt2}\sigma_+-\\
\frac1{\sqrt2}\sigma_y&=-\frac i{\sqrt2}\sigma_++\frac i{\sqrt2}\sigma_+-\tag4
\end{align*}
(this is an allowed transformation because
$$
\begin{pmatrix}
\frac1{\sqrt2}&\frac1{\sqrt2}\\
-\frac i{\sqrt2}&\frac i{\sqrt2}
\end{pmatrix}
$$
is unitary);
also the factor $4(\gamma+\tilde\gamma)$ becomes $2(\gamma+\tilde\gamma)$ because of the normalization of $\sigma_x,\sigma_y$ from Eq.(4).
Altogether this shows that the version from Eq.(2) with jump operators $\sigma_x,\sigma_y$ is equivalent to Eq.(3) meaning both are canonical forms of $\mathcal L$.
