When doing filtering, the Kinematic model comes into play only during the prediction step. That is, you predict that the next state is one which is reachable given your current kinematic constraints. I'll reiterate: this is not a hard constraint on the filter output. It simply makes unreachable states unlikely, with only the granularity available to the 2-D Gaussian distribution.
Then, the sensors are used to find which of the predicted states are likely given given the new information (GPS, etc).
If you must apply a kinematic model as a constraint (i.e., we cannot consider any state outputs which violate the kinematic model), then I know of only one option off hand: Use a particle filter to sample all nearby states, then apply a zero weight to all states which are not reachable given your Kinematics. Note, a particle filter is a general term with lots of different meanings.
A particle can represent all positions, $x,y,\theta$, e.g., a grid over the environment.
A particle can represent all nearby positions. In this case you sample the kinematic model many times with some noise, and propagate those forward. This implicitly honors the kinematic constraint.
A particle can represent possible paths (i.e., fastslam).
Then, you can re-weight each particle according to it's probability given the GPS measurements.
However, I recommend that you use the INS in place of the wheel odometry during the prediction step. It will likely be a reasonable proxy, especially over the short term. In fact, this is a very common practice in the literature.
To make this clearer, you predict your next pose $x^+$, based on your current pose $x^-$, and your INS measurement which gives you a $\delta x$ (because you integrate the accellerations). So your EKF_PREDICT function goes from:
$x_{t+1|t} = f_{kinematics}(x_{t}, u_t)$, with odometry $u_t$ instead of the control vector. To something like:
$x_{t+1|t} = x_{t}+\delta x$.
What have we lost? We have lost the ability to condition on the fact that the robot can only take certain types of movements, i.e., is kinematically constrained. To compensate for this, what I've suggested is you then condition the estimate based on the liklihood given the kinematic model. That means, during the EKF_UPDATE step.
To do this, you need to find the conditional probability:
$P(\delta x | x_{t}, u_t)$
