find roots with specified coefficients

Edited to provide general analytical answer.

With,

f = -16 x^2 (t - 2 z - 2 z y) + (t + 2 z + 2 z y) ((-4 + t - 2 z) (4 + t - 2 z)
    - 4 z^2 y^2);

the Discriminant of this cubic in t is

d = Discriminant[f, t] // Factor

(* 4096 (4 + 12 x^2 + 12 x^4 + 4 x^6 - 8 z^2 + 20 x^2 z^2 + x^4 z^2 - 8 y z^2 + 
   38 x^2 y z^2 - 8 x^4 y z^2 + y^2 z^2 + 20 x^2 y^2 z^2 - 8 x^4 y^2 z^2 + 4 z^4 + 
   8 y z^4 - 2 x^2 y z^4 + 2 y^2 z^4 + 2 x^2 y^2 z^4 - 2 y^3 z^4 + 8 x^2 y^3 z^4 + 
   4 x^2 y^4 z^4 + y^2 z^6 + 2 y^3 z^6 + y^4 z^6) *)

FindInstance suggests that d is non-negative everywhere,

FindInstance[d < 0, {x, y, z}, Reals]
(* {} *)

as does random sampling of points {x, y, z}. On this basis, it seems safe to conclude the roots of f, too lengthy to reproduce here,

rt = t /. Solve[f == 0, t]

are real, despite the fact that I appears in the rt[[2]] and rt[[3]], just as observed in the question. As a further test, Im[rt[[2]]] for several slices in, for instance, x.

Plot3D[Table[Im[rt[[2]]] // Chop, {x, -5, 5}], {y, -5, 5}, {z, -5, 5}, 
    PlotRange -> {-10^-8, 10^-8}]

is zero everywhere, again substantiating that all rt are real.

The general solution, based on 3.8.2 of Abramowitz and Stegun, is derived as follows.

a = CoefficientList[f, t]
(* {-32 z + 32 x^2 z - 32 y z + 32 x^2 y z + 8 z^3 + 8 y z^3 - 
    8 y^2 z^3 - 8 y^3 z^3, -16 - 16 x^2 - 4 z^2 - 8 y z^2 - 4 y^2 z^2, -2 z + 2 y z, 1} *)

q = Simplify[a[[2]]/3 - a[[3]]^2/9]
(* -(16/9) (3 + 3 x^2 + (1 + y + y^2) z^2) *)
r = Simplify[(a[[2]] a[[3]] - 3 a[[1]])/6 - a[[3]]^3/27]
(* 32/27 z (-9 x^2 (1 + 2 y) + (2 + y) (9 + (-1 - y + 2 y^2) z^2)) *)
Factor[q^3 + r^2];

The last expression is proportional to the negative of d and, therefore, is non-positive. We can then write the key expression,

s = ComplexExpand[(r + I Sqrt[-q^3 - r^2])^(1/3), TargetFunctions -> {Re, Im}]

which, after a bit of additional manipulation becomes

(* (c1 + c2^2)^(1/6) Cos[1/3 ArcTan[c2, Sqrt[c1]]] + 
   I (c1 + c2^2)^(1/6) Sin[1/3 ArcTan[c2, Sqrt[c1]]] *)

where

rul = {c1 -> (4096/729 (3 + 3 x^2 + (1 + y + y^2) z^2)^3 - 1024/729 z^2 (-9 x^2 (1 + 2 y) 
           + (2 + y) (9 + (-1 - y + 2 y^2) z^2))^2), 
       c2 -> 32/27 z (-9 x^2 (1 + 2 y) + (2 + y) (9 + (-1 - y + 2 y^2) z^2))}

The roots of f are, then,

re = {2 s[[1]] - a[[3]]/3, -s[[1]] + I s[[2]] Sqrt[3] - a[[3]]/3, 
    -s[[1]] - I s[[2]] Sqrt[3] - a[[3]]/3} // Simplify
(* {2/3 (z - y z + 3 (c1 + c2^2)^(1/6) Cos[1/3 ArcTan[c2, Sqrt[c1]]]), 
    1/3 (-2 (-1 + y) z - 3 (c1 + c2^2)^(1/6) Cos[1/3 ArcTan[c2, Sqrt[c1]]] - 
        3 Sqrt[3] (c1 + c2^2)^(1/6) Sin[1/3 ArcTan[c2, Sqrt[c1]]]), 
    -(2/3) (-1 + y) z - (c1 + c2^2)^(1/6) Cos[1/3 ArcTan[c2, Sqrt[c1]]] + 
        Sqrt[3] (c1 + c2^2)^(1/6) Sin[1/3 ArcTan[c2, Sqrt[c1]]]} *)

which explicitly answers the question posed by the OP. That this is equivalent to the solutions rt, two of which do contain I explicitly, can be seen from, for instance,

re /. rul /. {x -> 1, y -> 0.5, z -> 1}
(* {6.81903, -6.03763, 0.218601} *)
rt /. {x -> 1, y -> 0.5, z -> 1}
(* {6.81903 - 4.44089*10^-16 I, -6.03763 - 6.66134*10^-16 I, 0.218601 + 1.33227*10^-15 I}*)

or, if exact values are desired,

re /. rul /. {x -> 1, y -> 1/2, z -> 1}
(* {1/3 (1 + 4 Sqrt[31] Cos[1/3 ArcTan[(9 Sqrt[367])/8]]), 
    1/3 (1 - 2 Sqrt[31] Cos[1/3 ArcTan[(9 Sqrt[367])/8]] - 
        2 Sqrt[93] Sin[1/3 ArcTan[(9 Sqrt[367])/8]]), 
    1/3 (1 - 2 Sqrt[31] Cos[1/3 ArcTan[(9 Sqrt[367])/8]] + 
        2 Sqrt[93] Sin[1/3 ArcTan[(9 Sqrt[367])/8]])} *)

ComplexExpand[Simplify[rt /. {x -> 1, y -> 1/2, z -> 1}]] // FullSimplify

produces the same answer.

roots = t /. {Roots[-16 x^2 (t - 2 z - 2 z y) + (t + 2 z + 
            2 z y) ((-4 + t - 2 z) (4 + t - 2 z) - 4 z^2 y^2) == 0, t] // 
      ToRules} // Simplify;

roots // LeafCount

(*  861  *)

To get a representation without I use ComplexExpand

roots2 = ComplexExpand[roots, TargetFunctions -> {Re, Im}] // Simplify;

roots2 // LeafCount

(*  25973  *)

The first representation using I is much less complicated.

For your specific example,

ex1 = roots /. {x -> 1, y -> 1/2, z -> 1} // FullSimplify

(*  {Root[9 - 41 #1 - #1^2 + #1^3 &, 3], 
 Root[9 - 41 #1 - #1^2 + #1^3 &, 1], 
 Root[9 - 41 #1 - #1^2 + #1^3 &, 2]}  *)

ans1 = ex1 // N

(*  {6.81903, -6.03763, 0.218601}  *)

Converting the Root expressions

ex1R = ex1 // ToRadicals // ComplexExpand // Simplify

(*  {1/3 (1 + 4 Sqrt[31] Cos[1/3 ArcTan[(9 Sqrt[367])/8]]), 
 1/3 (1 - 2 Sqrt[31] Cos[1/3 ArcTan[(9 Sqrt[367])/8]] - 
    2 Sqrt[93] Sin[1/3 ArcTan[(9 Sqrt[367])/8]]), 
 1/3 (1 - 2 Sqrt[31] Cos[1/3 ArcTan[(9 Sqrt[367])/8]] + 
    2 Sqrt[93] Sin[1/3 ArcTan[(9 Sqrt[367])/8]])}  *)

For the second form

ex2 = roots2 /. {x -> 1, y -> 1/2, z -> 1} // FullSimplify

(*  {(1/3)*(1 + 4*Sqrt[31]*
          Cos[(1/3)*ArcTan[(9*Sqrt[367])/
                  8]]), (1/3)*
     (1 - 2*Sqrt[31]*
          Cos[(1/3)*ArcTan[(9*Sqrt[367])/
                  8]] - 2*Sqrt[93]*
          Sin[(1/3)*ArcTan[(9*Sqrt[367])/
                  8]]), (1/3)*
     (1 - 2*Sqrt[31]*
          Cos[(1/3)*ArcTan[(9*Sqrt[367])/
                  8]] + 2*Sqrt[93]*
          Sin[(1/3)*ArcTan[(9*Sqrt[367])/
                  8]])}  *)

This is identical to earlier result

ex2 === ex1R

(*  True  *)

Demonstrating that the original Root expression is equivalent

(ex2 // N) == ans1

(*  True  *)