Abstract: Parameter tuning of Proportional Integral (PI) controllers was a key step in the automatic control of canals, directly affecting system stability, control accuracy, and regulation smoothness. However, traditional PI tuning methods suffered from high computational costs, heavy reliance on engineering experience, and difficulty in balancing multiple performance criteria, severely restricting practical applications. Therefore, this study proposed an efficient and reliable PI parameter tuning method to solve the above problems and improve the overall control performance of canal systems. A two-stage Bayesian optimization method was designed for this purpose. In the first stage, the system settling time was taken as the core evaluation index. A stability discriminant function was established based on a Gaussian process surrogate model with an anisotropic composite kernel, combining the radial basis function and the Matérn kernel to capture the nonlinear behavior of settling time across the parameter space. This function characterized the stability of the PI control system under different parameter combinations. A Bayesian optimization algorithm guided adaptive parameter sampling to screen out the feasible region that ensured stable canal system. The discrimination threshold was determined by the optimal settling time obtained from the first-stage optimization multiplied by a safety factor, with prediction uncertainty explicitly incorporated to avoid misclassification of marginally stable points. In the second stage, a comprehensive performance index was constructed within the feasible region by the weighted combination of the normalized Integral Absolute Error(I_IAE) and the normalized Total Variation(T_TV) of the control action. A dynamic acquisition strategy, combining Expected Improvement and the Lower Confidence Bound criteria with a linearly decreasing exploration weight, was adopted to guide the Bayesian optimization process, achieving coordinated improvement of control accuracy and regulation smoothness. The proposed method was validated using the first long reach of the gentle slope test canal recommended by the American Society of Civil Engineers (ASCE). Particle Swarm Optimization (PSO) and standard Bayesian optimization were used as benchmark methods. All simulations were performed using a one-dimensional unsteady flow model based on the Saint-Venant equations. Results showed that the proposed method quickly and accurately determined the stability boundary of PI control parameters, requiring only about 50 simulation evaluations to obtain the stable feasible region. Within the feasible region, 100% of the sampled parameter combinations satisfied the stability requirement, while outside the region only 2% did so, demonstrating the high precision of the identified boundary. Compared with PSO, the proposed method reduced settling time by 34.5%, decreased I_IAE by 15.8% and T_TV by 9.9%, and cut total evaluations by 86.7%, significantly improving tuning efficiency. Compared with standard Bayesian optimization, the proposed method reduced T_TV by 32.6%, effectively suppressing control oscillations during canal operation and improving system smoothness, although it led to a slight increase (about 12%) in the Integral Absolute Error, indicating a successful trade-off between control accuracy and smoothness. In the multi-objective optimization process, when the weight coefficient of the comprehensive performance index was in the range of 0.3 to 0.7, the method could effectively balance control accuracy and regulation smoothness, meeting practical needs of canals. By ensuring stable canal operation, the proposed two-stage Bayesian optimization method effectively overcame the high computational cost and experience-dependent tuning of traditional approaches. It significantly improved tuning efficiency and overall system control performance, providing a reliable theoretical basis and technical support for canal automation.