I post here another way of doing this. Since the algorithm below builds a block matrix I suspect it is slower for large matrices than the previous answer (this was also indicated on my machine, when I tried it on some examples). My coding skills are however very limited, so the reason it is slower might be my fault.
getIntersectionbasis[l1_, l2_] := Module[{n, c1, c2, c3},
n = Dimensions[l1][[2]];
c1 = ArrayFlatten[{{l1, l1}, {l2, 0*l2}}];
c2 = RowReduce[c1];
c3 = Select[Take[#, n] == Table[0, {j, 1, n}] &][c2];
Map[Take[#, -n] &, c3]
]
The algorithm works like this: Build matrices $A$ and $B$ out of the row vectors in $\ell_1$ and $\ell_2$, and put them into a block matrix
$$C_1=\begin{pmatrix} A & A\\ B & 0\\ \end{pmatrix}.$$
Next, use row reduction on $C_1$ to obtain $C_2$. When that is done, construct $C_3$ from the rows in $C_2$ containing only zeros in their first $n$ entries (i.e. in the block where $B$ was). Finally, omitting the first $n$ zeros of each row in $C_3$ will leave you with a basis of the intersection of $\ell_1$ and $\ell_2$.