The answer from S tomio has the essential ingredient, that is to identify the coordinate axes from the symmetry properties.
However, one has to be careful about the conventions. Note already that in the 4th edition of Boyd (2020) one has a $3m$ susceptibility table (i.e. a "$d_{il}$ matrix") that conflicts with yours:
This $d_{il}$ matrix turns out to be compatible with the tensor element relations in the table posted by S tomio, which also appears in Boyd 4th edition, which asserts that the mirror plane is perpendicular to the $\hat{x}$ direction (i.e. is parallel to the $yz$ plane).
To orient ourselves, let's see why, for instance, the $d_{11} = d_{111}$ element vanishes.
This is proportional to the $x$-component of the polarization when the electric field $\mathbf{E}$ is also aligned in the $x$ direction (i.e. $\mathbf{E} = E \hat{x}$).
In equation form this reads
$$
d_{111} = \hat{x} \cdot P(\hat{x},\hat{x}) \propto \hat{x} \cdot P(E\hat{x},E\hat{x}) = \hat{x} \cdot P(\mathbf{E},\mathbf{E})
$$
However, if the crystal structure is unchanged under the reflection $R: x \to -x, y \to y, z \to z$, then the effect of reflecting the electric fields on the resulting polarization vector should be the same as reflecting the polarization, i.e.
$$
P(R \mathbf{E}, R \mathbf{E}) = R\left(P(\mathbf{E},\mathbf{E})\right)
$$
Now under the reflection $R$ we have $R\hat{x} = -\hat{x}$, and by the same principle the $x$-component of the polarization is also inverted, i.e.
$$
\hat{x} \cdot \left(RP\right) = \left( R\hat{x} \right) \cdot P = \left(-\hat{x}\right)\cdot P = - \left( \hat{x} \cdot P \right)
$$
where on the first step we exploit the identity $\mathbf{w} \cdot \left( O\mathbf{v}\right) = \left(O^T \mathbf{w}\right)\cdot \mathbf{v}$ for orthogonal transformations $O$ and that reflections are orthogonal transformations that additionally satisy $R = R^T$.
Putting this together we get
$$
d_{111} = \hat{x} \cdot P(\hat{x},\hat{x}) = \hat{x} \cdot P(-\hat{x},-\hat{x}) = \hat{x} \cdot P(R\hat{x},R\hat{x}) = \hat{x} \cdot R \left(P(\hat{x},\hat{x})\right) = - \hat{x} \cdot P(\hat{x},\hat{x}) = -d_{111}
$$
so that $d_{111} = d_{11}$ necessarily vanishes.
By the same argument, any tensor element $d_{ijk}$ containing just one "1" index, e.g. $d_{122} = d_{12}$ or $d_{123} = d_{14}$, also necessarily vanishes.
You can check that this checks out on the table from Boyd that I reproduce above.
Another convention is that the high symmetry axis should align with the $z$-axis.
In the $3m$ case this is a 3-fold axis.
All of the symmetry operations of the $3m$ point group are either rotations about this axis or reflections about planes containing this axis, so the $z$-coordinate is unmodified by any point group operation.
Therefore there can be no symmetry relations between tensor elements $d_{ijk}$ that modify any "3" indices.
You can check that this is the case in the table reproduced in S tomio's post.
We for instance have $xzx = yzy$, where the 3-fold rotational symmetry of the point group can "mix" the $x$ and $y$ components, but leaves the middle $z$ index unchanged.
Similarly there are no symmetry relations relating the $zzz$ tensor component to any other component, since all the symmetry operations have no effect on the component $\hat{z} \cdot P(\hat{z},\hat{z})$. This is reflected in the Boyd's $d_{il}$ matrix for the $3m$ group I reproduce above, where the $d_{33} = d_{333}$ component appears only once.
So we have established a relationship between the $xyz$ coordinate system that the $d_{il}$ are referenced to and the point group symmetry of the crystal by looking at the constraints among the $d_{il}$ induced by the crystal symmetries.
Namely, the 3-fold symmetry axis coincides with the $z$-axis, and one of the three equivalent reflection planes coincides with the $yz$ plane.
The other two reflection planes can be generated by $\pm120^\circ$ rotations of the $yz$ plane.
From inspection of the crystal structure using the visualizer at:
https://next-gen.materialsproject.org/materials/mp-5020/#crystal_structure ,
I am nearly certain that the 3-fold symmetry axis intersects a barium ion and is parallel to the c-axis of the conventional unit cell, i.e. the one with a length of 7.01 Angstroms unequal to the length 5.70 Angstroms of the other two.
The reflection plane contains the symmetry axis and (again I am nearly sure) a neighboring titanium ion. This description is sufficient for establishing a coordinate system from the positions of the ions, which I believe answers your question.
You however have an issue I mentioned earlier, which is that the $d_{il}$ matrix you present does not conform to the table presented by S tomio or the $d_{il}$ matrix from Boyd.
I believe the issue is that the $d_{il}$ matrix you use is either not actually the right table for the $3m$ point group, or the $d_{il}$ matrix you use is defined differently than the convention used in Boyd 4th edition.
You will need to check this.
The other possibility I considered was that it uses the non-standard convention of putting the $zx$ plane as a reflection plane.
This would explain why for your $d_{il}$ matrix the $d_{12} = d_{122}$ element is nonzero and the $d_{21} = d_{211}$ element vanishes, but this would not explain then why the $d_{22} = d_{222}$ does not vanish.
