The natural boundary integral equation (NBIE) is developed to calculate potential derivatives for potential problems with anisotropic media. Firstly, the governing equation of the two-dimensional anisotropic potential problem is transformed into standard Laplace equation by a coordinate transformation method. Then a potential derivative boundary integral equation named as NBIE is extended to solve the anisotropic potential problem. The most important virtue of the NBIE is that the singularity of the integral kernel function is reduced by one order in comparison with the conventional potential derivative boundary integral equation(CDBIE). Therefore the new potential derivative boundary integral equation only contains strongly singular integrals rather than hyper-singular integrals. Thus the NBIE can calculate more accurate potential derivative results for both boundary nodes and interior points. Moreover, in combination with the analytical integral regularization algorithm of nearly singular integrals, the NBIE can obtain more accurate potential derivatives of interior points very close to the boundary than the CDBIE. Numerical examples on heat conduction in anisotropic media demonstrate the accuracy and efficiency of the NBIE.