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.
I'd like to write a better response, but I must be brief.
For now, let me offer some places to read. Long story short, it is predicted that there's a relationship between special values of $p$-adic $L$-functions and syntomic regulators (which are the analogue of Beilinson's regulators in the $p$-adic world).
The beautiful paper of Manfred Kolster and Thong Nguyen Quong Do is, I think, a very readable resource.
The best results I know in this direction are Besser's papers here and here, which use rigid syntomic cohomology.
Besser's overview talk at the conference in Loen (notes available here) was a real joy.
Best Answer
Conjecturally, the answer is yes, but the amount of work required is not trivial at all. The general set-up is roughly as follows: the special values of $L$-functions (in your case, for Tate motives) are predicted by the Tamagawa Number Conjecture, and by the Equivariant Tamagawa Number Conjecture (ETNC) when one wishes to incorporate the action of some Galois group (as you do). The ETNC links the leading term at 0 of the $L$-function to the determinants of some cohomological complexes. However, in order to do this coherently for any locally constant character (again, as you want to), one needs rather strong hypotheses on the complexes involved: namely, they need to be semisimple at $\rho$ if one wishes to interpolate at $\rho$. Here, the bad news start: showing that a complex is semisimple at $\rho$ amounts to Leopoldt's conjecture for Tate motives so is in general very hard. But things should be fine in your favourite case. In this way, you can construct (leading terms of) $p$-adic $L$-functions for Tate motives.
If you are ready to admit all conjectures, all of this has been known for quite some time, see for instance B.Perrin-Riou Fonctions $L$ $p$-adiques des représentations $p$-adiques section 4.2 Conjecture CP(M). For a more recent and more flexible formulation, I recommend D.Burns and O.Venjakob On the leading terms of Zeta Isomorphism... section 3.2.2 and 3.2.3.