I don't know the history, but Wikipedia dates some of the conepts to work by Lorentz(1853-1928).
The key is this:
To characterize a two-port network requires that we relate the terminal quantities V1, V2, I1, and I2, out of which two are independent.
You can pick any two of those variables as independent and then the other two as dependent variables. Because the network is linear, once you know the relationship between one set of independent and dependent variables, you can do matrix transformations to find the relationship between some other choice of variables.
Depending which you pick, there's a different name for the matrix that relates them.
If you pick \$I_1\$ and \$I_2\$ as independent, then you get the "Z", or impedance, matrix and the equation
$${\bf V}={\bf Z I},$$
shorthand for
$$V_1 = Z_{11}I_1 + Z_{12}I_2$$
$$V_2 = Z_{21}I_1 + Z_{22}I_2$$
If you pick \$V_1\$ and \$V_2\$ as independent, then you get the "Y" or admittance matrix, and the equation
$${\bf I}={\bf Y V}.$$
Note that if we premultiply the two sides of this equation with the inverse of \$\bf Y\$ we have
$${\bf Y}^{-1}{\bf I}={\bf V},$$
which tells us that \${\bf Y}^{-1}\$ is just another name for \$\bf Z\$.
If you pick one current and one voltage as independent then that's considered a "hybrid" representation and the matrix is named "\$\bf H\$".
You can even pick linear combinations of the I and V variables as the independent and dependent variables and get things like the S-matrix ("scattering matrix") and S-parameters used in RF work.
Usually we pick a representation that makes our circuit easy to understand. Or we might initially analyzer a problem using a representation that makes it easy to understand the function of some subnetwork, and then use well-known transformations between the different matrices to get a form that's easier to solve, or easier to evaluate how it will combine with other subnetworks in the full network.
For example, we often represent a common-emitter BJT circuit using the H matrix because we know the collector current is very nearly proportional to the base current, but that it varies slightly due to the collector voltage (this is exactly the circuit given in your example, and why the BJT current gain \$\beta\$ is sometimes called \$h_{fe}\$).
Edit
I guess you are trying to focus on this part of your question?
But when I think of some one writing the following:
V1 = hi×I1 + hr×V2
my mind gets locked.
You can see it in the circuit model (b) that you posted.
The output current is the sum of the current through two parallel branches. One is a CCCS with gain \$h_f\$, so the current through that branch is \$h_f I_1\$. The other is a conductor with conductance \$h_o\$, so the current through that branch is \$h_o V_2\$. By KCL, the total current into port two (\$I_2\$) is therefore
$$I_2 = h_f I_1 + h_o V_2.$$
Like Neil points out in the comments, this means that \$h_f\$ must be a unitless quantity and \$h_o\$ must have units of conductance.
Edit 2
You found this definition of the superposition principle:
This is simply wrong, or at least it's a very incomplete version of the superposition principle.
First, not all resistors are linear, and a circuit containing a nonlinear resistor can't be solved by superposition.
Second, other elements than those specified can be included in a linear circuit, such as current sources, linear capacitors and inductors, and linear controlled sources (VCVS, VCCS, CCVS, or CCCS).
Probably this definition was put in a book at a point where voltage sources and linear resistors were the only components that had been introduced. Later, the author should have given a more complete version to allow for other component types.