Phase unwrapping (PU) represents a crucial bottleneck restricting the practical application of InSAR technique. According to the relationship between the same relative elevation and each interferometric phase differential, a new robust branch and bound pure integer programming phase unwrapping (BB-PIP-PU) algorithm for multi-baseline InSAR based on operations research is proposed. This algorithm constructs an N dimensional space with the ambiguity number as the axis. It takes the intercept of the ith ambiguity number axis as the objective function, and the directed ray intersecting the N-1 dimensional plane as the constraint condition. Firstly, it takes the optimal solution of ith ambiguity number axis as the solution starting point. Then, it increments by one pixel successively to verify whether the solution set composed of the remaining ambiguity numbers meets the integer condition, until the condition is established. Finally, the optimal integer solution of the N ambiguity numbers is determined. To improve noise robustness, a rectangular window is innovatively designed to replace ambiguity number at the center mutation point with ambiguity number with the highest frequency in the rectangular window. The significant advantages of this algorithm are that it has better unwrapping ability in the abrupt topographic change region, compared with the branch-cut, LS and MCF algorithms. Moreover, it weakens the baseline requirement of interferometric pairs, compared with the CRT algorithm. Finally, experimental results on both simulated and real data illustrate the good performance of our approach in terms of the effectiveness, universality, reliability, and noise robustness ability.