Let me first give you a heuristic "reason", why the regulator in the class number formula looks different from the regulator in the Birch and Swinnerton-Dyer conjecture. It is often more convenient (and more canonical) to combine the regulator and the torsion term in the Birch and Swinnerton-Dyer conjecture: one chooses a free subgroup $A$ of the Mordell-Weil group of $E(K)$ with finite index and looks at the quantity $R(A)/[E(K):A]^2$, where $R(A)$ is the absolute value of the determinant of the Néron-Tate height pairing on a basis of $A$. As you can easily check, the square in the denominator insures that this quantity is independent of the choice of $A$. Now, in the class number formula, the torsion term is replaced by the number of roots of unity in $K$ and that term is not squared. So for such a canonical formulation to be possible, it is reasonable to expect that not the regulator of the number field should be defined through a symmetric pairing on the units, but its square. Then you could make the same definition as in the elliptic curves case, and it would be independent of the choice of finite index subgroup.
Now, that we have established this, there are several ways to bring the two situations closer together.
You could do the naïve thing: take the matrix $M=(\log|u_i|_{v_j})$, where $u_i$ runs over a basis of the free part of the units (or more generally $S$-units for any set of places $S$ which includes the Archimedean ones), and where $v_j$ runs over all but one Archimedean place (or more generally all but one place in $S$), v_0, say. The absolute values have to be suitably normalised (see e.g. Tate's book on Stark's conjecture). Now take the symmetric matrix $MM^{tr}$ and define this to be the matrix of a new pairing on the units. In other words, you would define your symmetric pairing as
$$
(u_1,u_2) = \sum_{v\in S\backslash\{v_0\}} \log|u_1|_v\log|u_2|_v.
$$
Then, it's clear that you have a symmetric pairing and that the determinant of that pairing with respect to any basis on the free part of the units is $R(K)^2$. As I mentioned above, the square was expected.
Depending on what you want to do, this pairing might not be the best one to consider. For example if now $F/K$ is a finite extension and you consider the analogous pairing on the $S$-units of $F$ and restrict it to $K$, then it's not clear how to compare it to the pairing on $K$. In the elliptic curves case by contrast, the former is $[F:K]$ times the latter. To fix this, we can make the following definition:
$$
(u_1,u_2) = \sum_{v\in S} \frac{1}{e_vf_v}\log|u_1|_v\log|u_2|_v,
$$
where $e$ and $f$ are the absolute ramification index and residue field degree respectively. With this definition, the compatibility upon restriction to subfields is the same as in the elliptic curves case. On the other hand, the relationship with the actual regulator is slightly less obvious. It is however quite easy to show (and I have done it in http://arxiv.org/abs/0904.2416, Lemma 2.12, in case you are interested in the details) that the determinant of this pairing is given by
$$
\frac{\sum e_vf_v}{\prod e_vf_v}R(K)^2,
$$
with both the sum and the product again ranging over the places in $S$.
So I guess, the moral is that one shouldn't seek analogies between the BSD-formula and the class number formula but rather between the BSD and the square of the class number formula (note that also sha, which is supposed to be the elliptic curve analogue of the class number, has square order whenever it is finite). There is also a corresponding heuristic on the analytic side.
The method explained in Husemöller's book on elliptic curves is as follows:
Take a general quartic $v^2=f_4(u)=a_ou^4+a_1u^3+a_2u^2+a_3u+a_4$, and let
$$u=\frac{ax+b}{cx+d}\qquad v=\frac{ad-bc}{(cx+d)^2} y$$
Substituting you get:
$$v^2=\frac{(ad-bc)^2}{(cx+d)^4}y^2=f_4\bigg(\frac{ax+b}{cx+d}\bigg)$$
which implies
$$(ad-bc)^2y^2=f_4\bigg(\frac{ax+b}{cx+d}\bigg)(cx+d)^4=\sum_{i=0}^4a_i(ax+b)^{4-i}(cx+d)^i=$$
$$=c^4f_4\bigg(\frac{a}{c}\bigg)x^4+f_3(x)$$
where $f_3(x)$ is a cubic polynomial whose coefficient of $x^3$ is $c^3f'_4(a/c)$. For $a/c$ a simple root of $f_4$ and $ad-bc=1$, this leaves the cubic equation $y^2=f_3(x)$.
From here you can use the tools you mention in the question to take care of the cubic.
Best Answer
I proved that if $E/\mathbf{Q}$ is given using by a minimal Weiestrass equation, then
$ \#E(Z) \le C^{\text{rank} E(Q) + n(j) + 1} $
where $n(j)$ is the number of distinct primes dividing the denominator of the $j$-invariant of $E$ and $C$ is an absolute constant. This is in J. Reine Angew. Math. 378 (1987), 60-100.
Mark Hindry and I proved that if you assume the abc conjecture, then you can remove the n(j) in the above estimate. This is in Invent. Math. 93 (1988), 419-450. It is a conjecture due to Lang.
The papers contain more general results for (quasi)-S-integral points over number fields.