Many researchers calculated strain rate of significant differences from the same GPS measurement data. In this paper, we use the Kriging method in geostatistics to GPS velocity field. Interpolating the scattered GPS velocity data of Qinghai-Tibetan plateau and its adjacent areas to grid point values by Kriging, we calculate the strain rates from these nodal values in each grid cell similar to derivative of shape functions (essentially Lagrange interpolation function) in finite element algorithm, and obtain the stable distribution of strain rate field in Qinghai-Tibetan plateau. The results show that the main part of Qinghai-Tibetan plateau is in the state of compression in north-south direction, and extension in west-east direction. On the contrary, in the eastern part of Tibet, the strain rate is compressive in west-east and extensional in north-south direction. The orientations of principal strain rates are consistent with those of the P axis and T axis in focal mechanism. The high values of maximum compressive principal strain rates are located in the Himalayan main boundary thrust zone (MBT) and the surrounded regions. The maximum extensive principal strain rates are higher than those of the compressive ones in the main part of the interior of Qinghai-Tibetan plateau. Also, the surface dilation strain rate shows that it is in the state of surface compression in Himalayan and its surrounded areas, and in the state of surface extension in the interior of Qinghai-Tibetan plateau. The distribution of maximum shear strain rate clearly displays the outlines of some main active fault zones. The result of the strain rate in this study suggests that the contemporary tectonic strain of Tibet inherits the long term geological deformation.