According to the weight character of quaternary self-orthogonal codes, the relations between two dimensional optimal self-orthogonal codes and their weight distribution are discussed. By introducing two definitions of defying vector and projective weight vector of quaternary linear codes and using the matrix constructed with simplex codes, the relations between two dimensional optimal self-orthogonal codes and the non negative integral solutions of equation systems are setup. And, the existence problem of two dimensional optimal self-orthogonal codes is changed into the problem of determining the non negative integral solutions of equation systems.For the given code length, firstly, the minimum distance of two dimensional optimal self-orthogonal codes is determined by Griesmer bound, then, the generator matrices and the weight polynomials of all optimal self-orthogonal codes of this length are determined through solving integral equation systems. According to generator matrices obtained, equivalent relations among these optimal self-orthogonal codes are discussed by using elementary row transformations of matrices, coordinates permutations of vectors and conjugate transformations of elements. Finally the classification of two dimensional optimal self-orthogonal codes is solved, and the generator matrices of these non equivalent two dimensional optimal self-orthogonal codes and their weight polynomials are also presented. polynomials are also presented.