An algorithm for the least squares matching of overlapping 3D surfaces is presented. It estimates the transformation parameters between two or more fully 3D surfaces, using the Generalized Gauss-Markoff model, minimizing the sum of squares of the Euclidean distances between the surfaces. This formulation gives the opportunity of matching arbitrarily oriented 3D surfaces simultaneously, without using explicit tie points. Besides the mathematical model and execution aspects we give further extension of the basic model. The first extension is the simultaneous matching of sub-surface patches, which are selected in cooperative surface areas. It provides a computationally effective solution, since it matches only relevant multi-subpatches rather than the whole overlapping areas. The second extension is the matching of surface geometry and its attribute information, e.g. reflectance, color, temperature, etc., under a combined estimation model. We give practical examples for the demonstration of the basic method and the extensions.