The vertical derivative of potential field is widely used in processing and interpreting of potential field data. The methods of calculating each order vertical derivative of potential field can mainly be divided into two kinds: one is calcu- lated in the frequency domain, the other is in the space domain. Because of the problem of poor stability for these methods, this paper describes the ISVD algorithm for calculating each order vertical derivative of potential field. The result of process- ing and analyzing bouguer gravity anomalies of single geological model and the vertically integrated geological models causing by ISVD algorithm prove that each order vertical derivative of potential field by ISVD algorithm calculating, especially the high order vertical derivative, has better stability than it calculating in the frequency domain or space domain, and it can be used to effectively identify shallow small-scale geological boundaries. In order to verify the effect of ISVD algorithm applying to the actual data, each order vertical derivative of bouguer gravity anomaly data of a block by ISVD algorithm calculating i- dentify shallow small-scale geological boundaries and determine the fault system combining with the geological background information. The results prove that effect of determining shallow small-scale geological boundaries basing on each vertical derivative of potential field calculated by ISVD algorithm is better.