Why not to just use  a formula of taylor series to calculate the coefficients behind x?

You can basically do a readdata function as proposed by mmcdara. It is a good function to work with that kind of files

data := readdata(strDir, float, 2) # strDir is your file directory. 
#Read data detects columns automatically. You could have 1 2 3 .. columns, as you choose.

time:= data[1 .. -1, 1]; # the data is a list so to work accordingly
time:= convert(time, Vector); # could be useful to translate a list to vector
strain:=  data[1 .. -1, 2]*10^19;
strain:= convert(strain, Vector);

dataplot(time, strain, opts);# opts - whatever options you want

Your parameter a is basically can go outside of the integration.  You can basically integrate your function without a and then and then multyply be a

> f := (x, y) -> 1/2*sinh(y - x^2) + 1/2*tanh(x - y^3):
> sol := evalf(int(int(f(x, y), x = -6 .. 5), y = -5 .. 5)):
> g := a*sol
   g = -1.352318125*10^16*a