  We study $\mathbb{Z}_3$ symmetry-protected topological (SPT) phases in
one-dimensional spin systems with $Z_3 \times Z_3$ symmetry. We construct
ground-state wave functions of the matrix product form for nontrivial
$\mathbb{Z}_3$ phases and their parent Hamiltonian from a cocycle of the group
cohomology $H^2(Z_3\times Z_3,U(1))$. The Hamiltonian is an SU(3) version of
the Affleck-Kennedy-Lieb-Tasaki (AKLT) model, consisting of bilinear and
biquadratic terms of su(3) generators in the adjoint representation. A
generalization to the SU($N$) case, the SU($N$) AKLT Hamiltonian, is also
presented which realizes nontrivial $\mathbb{Z}_N$ SPT phases. We use the
infinite-size variant of the density matrix renormalization group (iDMRG)
method to determine the ground-state phase diagram of the SU(3)
bilinear-biquadratic model as a function of the parameter $\theta$ controlling
the ratio of the bilinear and biquadratic coupling constants. The nontrivial
$\mathbb{Z}_3$ SPT phase is found for a range of the parameter $\theta$
including the point of vanishing biquadratic term ($\theta=0$) as well as the
SU(3) AKLT point [$\theta=\arctan(2/9)$]. A continuous phase transition to the
SU(3) dimer phase takes place at $\theta \approx -0.027\pi$, with a central
charge $c\approx3.2$. For SU(3) symmetric cases we define string order
parameters for the $\mathbb{Z}_3$ SPT phases in a similar way to the
conventional Haldane phase. We propose simple spin models that effectively
realize the SU(3) and SU(4) AKLT models.
