Okay, that makes things a bit simpler. Since the points are already spaced and ordered along z, I think what I would recommend is using functions that are parameterized by the z coordinate, rather than having a separate parameter t. This will reduce the number of degrees of freedom in the problem.
In other words, the solution will be a pair of functions x(z) and y(z) that give the curve's x and y locations for a given z.
If we're using polynomials (which seems like a good place to start), the first thing to do is to pick a degree. Lower degrees will give a smoother curve, while higher ones will fit the data more closely...you'll have to experiment to find what works best. Degree 3 (cubic) is probably a good place to start.
Let's let 'n' be the degree. I'll use a_i and b_i to refer to the coefficients of the x and y polynomials, respectively, where 0 \<= i \<=
n. So the x and y polynomials we are trying to solve for are:
x(z) = a_0 + z * a_1 + z\\^2 * a_2 + ... + z\\^n * a_n
y(z) = b_0 + z * b_1 + z\\^2 * b_2 + ... + z\\^n * b_n
where the a_i and b_i are unknown. Let 'm' be the number of data points and x_j, y_j, z_j be their coordinates, where 1 \<= j \<= m. Then the constraints we are trying to satisfy are:
x(z_1) = x_1, y(z_1) = y_1
x(z_2) = x_2, y(z_2) = y_2
x(z_m) = x_m, y(z_m) = y_m
This forms two linear systems of m equations in n unknowns, one for x and one for y. Since m > n the system is overdetermined and doesn't have an exact soluton, but we can still use least squares to find the approximate solution with the smallest mean-squared-error.