Your probability changes with $t$ but as Michael said, you don't know how. linearly or not ? It looks like a model selection problem where your probablity $p$ :
$p=\Phi(g(t,\theta))$ may depend on a highly non linear $g(t,\theta)$ function. $\Phi$ is just a bounding function that guarantees between 0 and 1 probabilities.
A simple exploratory approach would be to try several probits for $\Phi$ with different non linear $g()$ and to perform a $g()$ model selection based on standard Information Criterias.
To answer your re-eddited question:
As you said using probit would imply numerical solutions only but you may use a logistic function instead :
Logistic function: $P[\theta(t+1)] = \frac{1}{1+\exp{(\theta(t)+\epsilon)}}$
Linearized by : $ \log{\frac{P}{1-P}} = \theta(t)+\epsilon $
I'm not sure how this can work under Kalman filter approach, but still believe that a non linear specification like $\theta(t+1)=a t^3 +bt^2+ct + d $ or many others without a random term will do the job. You already have randomnes in the bernouilli event (Markov Chain) and you are adding an additional source of it due to $\epsilon$. I suppose you agree that that parsimony is very important. Unless your main objective is to apply a given method (HMM and Kalman Filter) and not to give the simplest valid solution to your problem.