I have a IMU and GPS sensors on my PCB which contains a microcontroller. An IMU consists of a gyroscope, accelerometer and a magnetometer (the magnetometer is calibrated to compensate for the effects of hard and soft iron). I am using Madgwick's filter to estimate the rotation matrix needed to rotate the read out accelerometer values onto the North-East-Down (NED) coordinate system. After subtracting the gravity vector, I feed the obtained acceleration components to the Kalman Filter. Note that the Madgwick filter estimates the gyroscope bias.

The states of my Kalman filter are:

$$ \vec x = \begin{bmatrix} p_x & p_y & p_z & v_x & v_y & v_z \end{bmatrix}^T \tag 1$$

and the control vector is:

$$ \vec u = \begin{bmatrix} a_x & a_y & a_z & a_x & a_y & a_z \end{bmatrix}^T \tag 2$$

I constructed the state space model:

$$ \dot{\vec x} = \begin{bmatrix} 1 & 0 & 0 & T & 0 & 0 \\ 
                            0 & 1 & 0 & 0 & T & 0 \\
                            0 & 0 & 1 & 0 & 0 & T \\
                            0 & 0 & 0 & 1 & 0 & 0 \\
                            0 & 0 & 0 & 0 & 1 & 0 \\
                            0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix} 
  \vec x + \begin{bmatrix}  T^2 & 0 & 0 & 0 & 0 & 0 \\ 
                            0 & T^2 & 0 & 0 & 0 & 0 \\
                            0 & 0 & T^2 & 0 & 0 & 0 \\
                            0 & 0 & 0 & T & 0 & 0 \\
                            0 & 0 & 0 & 0 & T & 0 \\
                            0 & 0 & 0 & 0 & 0 & T \end{bmatrix} \vec u = A \vec x + B \vec u  \tag 3$$

I would like to expand my state space model to incorporate additional states that would estimate the accelerometer bias. However, I don't know how to do that. My intuition tells me I need compare my acceleration data to something to obtain the bias components, but I don't know what. As I said, I have a GPS, so I am assuming there is something there that I can use, but I don't know what.

So my question is, how can the accelerometer bias terms be estimated using the Kalman Filter?