A new approach for debris removal sequence optimization with multiple spacecraft is developed and tested. The motion planning of the multiple rendezvous involves complex Traveling Salesman Problem (TSP) combinatorial decision, in which both the continuous parameters of the orbital maneuvers and integer indices of the debris are determined. Expected value and variance of trajectory conditioned on preceding visit, not the state of target itself, are propagated and analyzed. With the metric, likelihood of expected trajectory can be given, and the mixed integer problem is reformulated into continuous parameter estimation with maximum likelihood. A non-population gradient based optimizer is then used to search the target sequence. The method of TSP optimization is compared to the continuous parameter estimator of Extended Kalman Filter. Using this method, the classic static TSP-LIB benchmarks with up to 100 nodes are tested. We then extend the method to the TSP with space dynamics. Effectiveness and efficiency of the proposed method are tested by a sequential visit of near-earth debris considering J2 effects with a single and multiple spacecrafts.