I am implementing precise orbit determination for LEO satellites by EKF using GNSS carrier phase measurements, and I have a question about how to estimate the phase ambiguity.

First, I implemented the following based on this post.

Create the modifier for each continuous phase measurement for each GNSS satellite.

Set a different name for each modifier and add the same modifier for each continuous measurement.

(Is there any mistake in usage?)

However, in this method, the size of the state vector would be larger by the amount of modifiers created. Therefore, If the estimation period is long and the number of measurement data increases, or if cycle slips occur frequently, the computation time will slow down due to the increased size of the state vector.

So, is there another way to implement this that avoids this problem? For example, creating only one modifier for each GNSS satellite, and resetting the estimate and covariance when a slip occurs, etc.

This is a known problem that should really apply only to batch least square estimation. As on both sides of a cycle slip the ambiguities are independent, we must estimate them as different parameters, so the size of the state vector effectively increases. Their is even a worst case: clock offsets, which should really be modeled as stochastic parameters (the quadratic drift being only a smoothing model that is not suitable for high accuracy).

As you tell you are using a Kalman filter, you should be able to keep the same parameter. You could use some tricks to bump the covariance for the parameter up when their is a cycle slip. Finding the proper value is probably tricky.

You can look at CovarianceMatrixProvider and set up you own implementation that will use the state date to return either a regular covariance or a bumped up covariance at cycle slip occurrences.