Parameter identification determines the essential system parameters required to build real-world dynamical systems by fusing crucial physical relationships and experimental data. However, the data-driven approach faces many difficulties, such as discontinuous or inconsistent time trajectories and noisy measurements. The ill-posedness of the inverse problem comes from the chaotic divergence of the forward dynamics. Motivated by the challenges, we shift from the Lagrangian particle perspective to the state space flow field's Eulerian description. Instead of using pure time trajectories as the inference data, we treat statistics accumulated from the Direct Numerical Simulation (DNS) as the observable. The continuous analog of the latter is closely related to the physical invariant measure, a stationary distributional solution to the continuity equation. The connection motivates us to build a regularized forward model in the form of a PDE and reformulate the original parameter identification problem as a data-fitting, PDE-constrained optimization problem. A finite-volume upwind scheme and the so-called teleportation regularization are used to discretize and regularize the forward problem. We present theoretical regularity analysis for evaluating gradients of optimal transport costs and introduce two different formulations for efficient gradient calculation. Numerical results using the quadratic Wasserstein metric from optimal transport demonstrate the robustness of the novel approach for chaotic system parameter identification.