The key is to do the Taylor expansion properly (I change the notation slightly, using initially $\delta$ instead $d$. I also highlight in $\color{red}{red}$ the terms in the Taylor expansion that differ to yours):

First, let's define $\delta V$:

$$\delta V := V(S_1(t_0)+\delta S_1, S_2(t_0)+\delta S_2, t_0+\delta t))- V(S_1(t_0), S_2(t_0), t_0))$$

Now let's do a Taylor expansion around zero (i.e. $S_1(t_0), S_2(t_0), t_0)$ ):

$$\delta V=\frac{\partial V}{\partial t} \delta t + \frac{\partial V}{\partial S_1} \delta S_1 + \frac{\partial V}{\partial S_2} \delta S_2 + \frac{1}{2}\frac{\partial^2 V}{\partial S_1^2} \color{red}{\delta S_1^2} + \frac{\partial^2 V}{\partial S_1 \partial S_2} \color{red}{\delta S_1 \delta S_2} + \frac{1}{2} \frac{\partial^2 V}{\partial S_2^2} \color{red}{\delta S_2^2}$$

We ignore all higher order terms (i.e. $\delta t \delta S_1$, $\delta t \delta S_2$ , $\delta t ^2$ etc.), assuming they go to zero as $\delta \to 0$). 

Now evaluating the terms:

$$\delta S_1^2=(a_1\delta t + b_1 \delta X_1)^2=a_1^2\delta t^2+2a_1b_1\delta t \delta X_1+b_1^2\delta X_1^2$$

Now we can argue as $\delta t \to 0$ that: 

- $a_1^2\delta t^2 \to 0$
- $2a_1b_1\delta t \delta X_1 = 2a_1b_1\delta t \sqrt{\delta t} X_1 = 2a_1b_1\delta t^{\frac{3}{2}}X_1\to 0$ (using the scaling property of Brownian motion, i.e. for any finite $\delta t > 0$, we have that $\delta X_1 := X_1(\delta t) \stackrel{d}{=}\sqrt{\delta t}X_1(1)$)
- $b_1^2\delta X_1^2 \to b_1^2 dt$ (using the scaling property again and computing: $\mathbb{E}[b_1^2\delta X_1^2] = b_1^2\delta t\mathbb{E}[X_1^2] = b_1^2\delta t $, $V(\delta b_1^2\delta X_1^2) = b_1^2\delta t^2 \mathbb{E}[X_1^2] \to 0$ because \delta t^2 \to 0 )

So using the above $b_1^2\delta S_1^2 \to b_1^2dt$

Using the same machinery, we get that $b_2^2 \delta t S_2^2 \to b_2^2 dt$. And we can use the same machinery to compute $\delta S_1 \delta S_2$:

$$\delta S_1 \delta S_2 = a_1a_2 \delta t^2 + a_1b_2 \delta t \delta X_1+b_1a_2 \delta t \delta X_1+ b_1 b_2 \delta X_1 \delta X_2$$

Using the results already shown, everything goes to zero except for the last term:

$$b_1b_2\delta X_1 \delta X_2= b_1b_2 X_1(\delta t) X_2 (\delta t) = b_1b_2 \sqrt{\delta}X_1(1) \sqrt{\delta} X_2 (1) = b1b2 \delta t X_1X_2 $$

We know that the variances are (by definition of Brownian motion) $V(X_1(t))=\mathbb{E}[X_1(t)^2]= t$ and $V(X_2(t))=\mathbb{E}[X_1(t)^2]= t$ and the correlation between $X_1(t)$ and $X_2(t)$ is $\rho$, which means that:

$$Cov(X_1(t), X_2(t))=\rho \sqrt{V(X_1)V(X_2)}= \rho t$$ 

We also know that:

$$Cov(X_1(t), X_2(t)) = \mathbb{E}[X_1(t) X_2(t)]-\mathbb{E}[X_1(t)]\mathbb{E}[X_1(t)] = \mathbb{E}[X_1(t) X_2(t)]=t\mathbb{E}[X_1 X_2] $$

Using the last two equations, we can say that:

$$\mathbb{E}[X_1 X_2]=\rho $$

Therefore, the expected value of the cross term $\delta S_1 \delta S_2$ is $b_1b_2\delta t \rho$ which goes to $b_1b_2 \rho dt$ as $\delta t \to 0$. The variance of this term goes to $0$ (using $V( b1b2 \delta t X_1X_2)=b_1^2b_2^2\delta t^2 V(X_1X_1)$ which goes to zero due to $\delta t^2 \to 0$), so this term converges to $b_1b_2 \rho dt$.

So the formula is (after grouping all terms):

$$dV =\left( \frac{\partial V}{\partial t} +  \frac{\partial V}{\partial S_1}a_1+ \frac{\partial V}{\partial S_2}a_2+\frac{1}{2} \frac{\partial^2 V}{\partial S_1^2}b_1^2+\frac{\partial^2 V}{\partial S_1\partial S_2}\rho b_1b_2+\frac{1}{2} \frac{\partial^2 V}{\partial S_2^2}b_2^2\right) dt + b_1\frac{\partial V}{\partial S_1}dX_1 +b_2\frac{\partial V}{\partial S_2}dX_2 $$